CN109740453A - A kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation - Google Patents
A kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation Download PDFInfo
- Publication number
- CN109740453A CN109740453A CN201811554885.4A CN201811554885A CN109740453A CN 109740453 A CN109740453 A CN 109740453A CN 201811554885 A CN201811554885 A CN 201811554885A CN 109740453 A CN109740453 A CN 109740453A
- Authority
- CN
- China
- Prior art keywords
- magnetic field
- data
- track
- component
- magnetostatic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
Earthquake monitoring field of the present invention is a kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation, comprising: reads the magnetic field data of satellite, and chooses valid data according to flag bit；The tranquiler magnetostatic track of geomagnetic activity is chosen according to geomagnetic index；The magnetic field three component seismic data of magnetostatic track is pre-processed by difference and wavelet transform, removal amplitude with the slow magnetostatic field part of variation, obtains the situation of change of magnetic field data greatly；Dimension-reduction treatment is carried out to pretreated magnetic field three component seismic data using principal component analysis, main feature in stick signal while removing redundancy；Continuous wavelet transform is carried out to the magnetic field data after dimensionality reduction, orbital energy intensity is defined by wavelet coefficient, and carry out abnormal track using it and extract.The present invention compensates for exclusive use, and wherein a certain component sufficiently cannot cannot reflect the deficiency that seismic precursor influences magnetic field data using data useful information and its anomaly extracting result comprehensively.
Description
Technical field
Earthquake monitoring field of the present invention, particularly the satellite magnetic field data seismic precursor to be a kind of based on wavelet transformation
Abnormal extraction method.
Background technique
Seismic precursor refers to the abnormal phenomenon occurred before earthquake occurs, the earth magnetism occurred including focus and neighbouring substance,
The geophysical anomalies such as ground electricity, gravity.Its observation method may include ground observation and moonscope two types, from global model
Ground observation Station distribution is more sparse from the point of view of enclosing, and to compensate for ground observation observation scope limited for moonscope, and needs to spread
If the deficiency of a large amount of stations.Lithosphere-atmosphere-Ionosphere coupling mechanism proposition simultaneously shows that seismic precursor generates different
Often ionosphere can be propagated to by lithosphere, atmosphere, and it is had an impact.Therefore, earthquake is extracted using satellite data
Precursory anomaly is feasible.Magnetic field can be abnormal variation during earthquake preparation, and the magnetic field data of Satellite observation includes
Scalar resultant field and three-component vector data, vector data reflect the situation of change in different directions magnetic field, and containing more can be with
The information of abnormal conditions before reflecting earthquake.
CN106021710A discloses satellite orbit abnormality recognition method before a kind of shake based on atmospheric ionization layer parameter, packet
Include unruly-value rejecting and the standardization of primary ionization parameter；Under different earthquake earthquake magnitude, it is close that probability is divided using different threshold values
Function is spent, satellite exception track identification model before shaking is established；The track towards multi-parameter is obtained using integrated evaluating method to identify
The contents such as model.For the features such as atmosphere ionization supplemental characteristic amount is huge, information dimension is high, gives a kind of seismic precursor and defend
The recognition methods of star orbital road exception, can be used for satellite orbit abnormality detection and discrimination of earthquake precursory.
CN107356969A discloses a kind of seismic precursor analysis method based on satellite thermal infrared data and GIS, uses
Thermal Infrared Data adds the mode of GIS to carry out earthquake prediction analysis, can improve accuracy rate, at the same for the data source of analysis compared with
Extensively, there is preferable research effect.All pairs of earthquakes can be occurred in the way of GIS for the research direction for having expanded seismic precursor
There may be the quantization factors of influence to be all included in spatial operation, is not limited solely to geology coefficient and architectonic
Analysis.
Mehdi Akhoondzadeh etc. analyzes the Ecuador on April 16th, 2016 using intermediate value and quartile method
Shake, carries out the magnetic field three component seismic data of the Swarm satellite group Alpha star on April 30,1 day to 2016 November in 2015 respectively
Processing, has extracted corresponding precursory anomaly.Angelo De Santis etc. is with analyzing the Nepal on April 25th, 2015
Shake daily adds up the east orientation magnetic-field component in BEFORE AND AFTER EARTHQUAKE one month extremely, and is fitted discovery by S curve and deposits
The phenomenon that abnormal number quicklys increase, and using acceleration maximum as precursory anomaly day.Above-mentioned technology is respectively to magnetic field three
Component data individually carries out handling or using only one of component progress treatment research, there are more redundancies in the data used
Information or the data used are not comprehensive, fail that the useful information contained in three component seismic data is sufficiently efficiently utilized to carry out earthquake
The extraction of precursory anomaly, abnormal results are only obtained by the data of a certain component, fail sufficiently comprehensive reflection seismic precursor pair
The influence in magnetic field.
Summary of the invention
Before a kind of satellite magnetic field data earthquake based on wavelet transformation is provided
Million abnormal extraction methods carry out dimension-reduction treatment to magnetic field three component seismic data by principal component analytical method, reduce redundant data simultaneously
The main feature for retaining three component seismic data, more useful informations are transformed into dimensionality reduction data, and is obtained by continuous wavelet transform
It is distributed to frequency domain energy at that time, further defines and seek orbital energy intensity, effectively extract the abnormal track in dimensionality reduction data,
After excluding non-shake influence factor, Earthquake Precursor Anomalies are obtained.The present invention can make full use of useful in the three component seismic data of magnetic field
Information, the comprehensive influence for reflecting seismic precursor to magnetic field data compensate for above-mentioned technology and are used alone wherein that a certain component is not
Sufficiently it cannot can reflect that seismic precursor influences not magnetic field data comprehensively using data useful information and its anomaly extracting result
Foot.
The invention is realized in this way
A kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation, this method comprises:
A, the magnetic field data of satellite is read, and valid data are chosen according to flag bit；
B, the tranquiler magnetostatic track of geomagnetic activity is chosen according to geomagnetic index；
C, the magnetic field three component seismic data of magnetostatic track is pre-processed by difference and wavelet transform, removes amplitude
The big and slow magnetostatic field part of variation, obtains the situation of change of magnetic field data；
D, dimension-reduction treatment is carried out to pretreated magnetic field three component seismic data using principal component analysis, removes redundancy
Main feature in stick signal simultaneously；
E, continuous wavelet transform is carried out to the magnetic field data after dimensionality reduction, orbital energy intensity is defined by wavelet coefficient, and
Abnormal track is carried out using it to extract；
F, excluding non-seismic factor influences caused exception；
G, the extraction result of Earthquake Precursor Anomalies is exported.
Further, step a includes: the magnetic field three-component B for reading the magnetic field data in moonscope resultx、ByAnd Bz,
And invalid data is removed according to the flag bit of data.
Further, step b according to geomagnetic index choose magnetostatic track include: by geomagnetic index ap and Dst removal by
Geomagnetic activity interferes biggish track, the condition that magnetostatic track meets are as follows: current time ap < 12, | Dst |≤20；It is first 23 small
When ap≤32, | Dst |≤30.
Further, step c includes first carrying out first-order difference processing respectively by track to the magnetic field three component seismic data of selection
Changes of magnetic field amount is obtained, then wavelet transform is carried out to data after difference processing respectively, removes low-frequency component therein.
Further, wavelet transform obtains the formula of every layer of radio-frequency component and low-frequency component are as follows:
Wherein, j is Decomposition order, and L indicates that low-frequency component, H indicate that radio-frequency component, g are low-pass filter, and h is high pass filter
Wave device, dB are magnetic field differential data, and the residual error of every track magnetic field three-component differential data and low-frequency component is pre-processed results.
Further, step d, described being handled using principal component analysis satellite magnetic field Data Dimensionality Reduction include:
The magnetic field three component seismic data of every track after pretreatment is expressed as according to time series:
Xi=[xi(1),xi(2)...,xi(m)], i=1,2 ..., n
Wherein m is track length, and n is data dimension, obtains matrixThe expression formula of Y are as follows:
The covariance matrix C of calculating matrix YYThe element γ of (n × n)uv, calculation formula are as follows:
Wherein, xiuAnd xivThe i-th row u column element and the i-th row v column element in respectively matrix Y；WithIt is respectively
The mean value of u column and v column element；
Calculate the characteristic value and feature vector of covariance matrix:
CY=R Λ RTWherein, Λ (n × n) is the characteristic value diagonal matrix arranged from big to small, and R (1 × n) is character pair
The list of feature values is shown as λ by the feature vector of value1,λ2,...,λn(λ1> λ2> ... > λn)；
The initial data of matrix Y is subjected to linear projection using matrix R and obtains principal component Φ,
Φ=RY=[Φ1,Φ2,...,Φn]T
Wherein, Φ1,Φ2,...,ΦnIt is the 1st to the n principal component (1 × m), the contribution rate of preceding k principal component is right by its
Characteristic value is answered to calculate:
The preceding k principal component that choosing makes contribution rate to 60% or more represents the main feature of magnetic field data, and as drop
Tie up the result of processing.
Further, step e carries out continuous wavelet transform to the magnetic field data after dimensionality reduction, defines track by wavelet coefficient
Energy intensity, and carry out abnormal track extraction using it and include:
Signal is transformed into time-frequency domain by time domain by continuous wavelet transform；
The energy intensity at signal each moment is calculated in time-frequency domain；
Data are averaged using sliding rectangular window, to ensure the aberrant continuation regular hour, it is different to exclude mutability
Often；
Reflect its anomalous variation situation using the maximum average energy intensity on this track.
Further, further includes: calculate the mean μ of all orbital energy intensityZAnd standard deviation sigmaZ, and threshold value is set as ThZ=
μZ+kk×σZ, when the energy intensity of track is greater than the threshold value of setting, which is considered as abnormal track.
Compared with prior art, the present invention beneficial effect is:
The present invention is based on the satellite magnetic field data Earthquake Precursor Anomalies extracting methods of wavelet transformation, obtain magnetic field by difference
Variable quantity removes the big magnetostatic field variation tendency of periodic amplitude using wavelet transform and obtains changes of magnetic field situation；Benefit
Dimension-reduction treatment is carried out to three-component magnetic field data with principal component analysis, extracts the main feature of three component seismic data, it will be more useful
Information is retained in dimensionality reduction data；Continuous wavelet transform finally is carried out to data and obtains frequency domain energy distribution at that time, it is further fixed
Justice simultaneously seeks orbital energy intensity, and the anomalous variation of prominent signal can obtain good anomaly extracting effect, and guarantee to extract
The abnormal duration excludes mutation sexual abnormality.The present invention can make full use of in the three-component of magnetic field by the process
Information effectively extracts the Earthquake Precursor Anomalies variation with the duration.
Detailed description of the invention
Fig. 1 is the satellite magnetic field data Earthquake Precursor Anomalies extracting method flow chart based on wavelet transformation；
Fig. 2 is Ecuador earthquake epicenter, influence area and survey region schematic diagram；
Fig. 3 is satellite magnetic field three component seismic data primitive curve, and wherein A is Bx component, and B is By component, and C is Bz component；
Fig. 4 a is satellite magnetic field three-component differential data and low-frequency component curve, wherein A is Bx component, and B is By component, C
For Bz component；
Fig. 4 b is satellite magnetic field three-component pre-processed results curve, wherein A is Bx component, and B is By component, and C is Bz points
Amount；
Fig. 5 is satellite magnetic field three component seismic data first principal component curve；
Fig. 6 a is satellite magnetic field dimensionality reduction data wavelet conversion coefficient result figure；
Fig. 6 b is satellite magnetic field dimensionality reduction data average energy intensity variation curve；
Fig. 7 is satellite magnetic field dimensionality reduction data-track energy intensity change curve；
Fig. 8 a is on 2 29th, 2,016 3 Earthquake Precursor Anomalies result figures of track；
Fig. 8 b is 1 Earthquake Precursor Anomalies result figure of track on March 14th, 2016.
Specific embodiment
In order to make the objectives, technical solutions, and advantages of the present invention clearer, with reference to embodiments, to the present invention
It is further elaborated.It should be appreciated that the specific embodiments described herein are merely illustrative of the present invention, it is not used to
Limit the present invention.
Extraction the present invention is based on wavelet transformation to satellite magnetic field data progress Earthquake Precursor Anomalies, reads satellite first
Magnetic field data, and valid data are chosen according to flag bit；The tranquiler track of geomagnetic activity is chosen further according to geomagnetic index；Pass through
Difference and wavelet transform pre-process the magnetic field data of magnetostatic track, and removal amplitude is big and changes slow magnetostatic field
Part obtains the situation of change of magnetic field data；Dimensionality reduction is carried out to pretreated magnetic field three component seismic data using principal component analysis
Processing, most important feature in stick signal while removing redundancy, makes in signal comprising more useful informations；To dimensionality reduction
Magnetic field data afterwards carries out continuous wavelet transform, defines orbital energy intensity by wavelet coefficient, prominent anomalous variation guarantees to mention
The abnormal duration is taken, excludes mutation sexual abnormality, and carry out abnormal track using it and extract；Non-seismic factor influence is excluded to lead
The exception of cause；Export the extraction result of Earthquake Precursor Anomalies.The present invention can be effectively using wavelet transformation to defending after dimensionality reduction
Star magnetic field data carries out Earthquake Precursor Anomalies and extracts.
It is shown in Figure 1, a kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation, including with
Lower step:
A, the magnetic field data of satellite is read, and valid data are chosen according to flag bit；
B, the tranquiler magnetostatic track of geomagnetic activity is chosen according to geomagnetic index；
C, the magnetic field three component seismic data of magnetostatic track is pre-processed by difference and wavelet transform, removes amplitude
The big and slow magnetostatic field part of variation, obtains the situation of change of magnetic field data；
D, dimension-reduction treatment is carried out to pretreated magnetic field three component seismic data using principal component analysis, removes redundancy
Main feature in stick signal simultaneously；
E, continuous wavelet transform is carried out to the magnetic field data after dimensionality reduction, orbital energy intensity is defined by wavelet coefficient, and
Abnormal track is carried out using it to extract；
F, excluding non-seismic factor influences caused exception
G, the extraction result of Earthquake Precursor Anomalies is exported；
Step a reads the magnetic field data of satellite, and chooses valid data according to flag bit, is to read moonscope result
In magnetic field three-component Bx、ByAnd Bz, it refers respectively to: north orientation, east orientation and vertical component, and removed according to the flag bit of data
Effective magnetic field data is used only in invalid data, and flag bit is the parameter for characterizing magnetic field data quality condition, can pass through mark
Whether position judgement goes wrong when the instrument is measured, and removal is existing method due to the wrong caused invalid data of apparatus measures,
Details are not described herein again.Since the time of measurement data and space can change with satellite flight, therefore according to track to data
Carry out processing and anomaly extracting.
Step b chooses magnetostatic track according to geomagnetic index, is dry by geomagnetic activity by geomagnetic index ap and Dst removal
Biggish track is disturbed, the exception extracted is avoided to be caused by geomagnetic activity.Geomagnetic index reflects the geomagnetic activity in a period of time
Intensity, magnetostatic track need the condition met are as follows: current time ap < 12, | Dst |≤20；Preceding 23 hour ap≤32, | Dst |
≤30。
Step c is the magnetic field three-component first to selection using difference and wavelet transform satellite magnetic field data prediction
Data carry out first-order difference by track respectively and handle to obtain changes of magnetic field amount, then carry out discrete wavelet transformer to differential data respectively
It changes, removes low-frequency component therein, the situation of change of prominent data.Satellite observation magnetic field be earth's magnetic field, wherein comprising main field,
Four part of lithospheric magnetic field, variation magnetic field and induced magnetic field.Earth's magnetic field amplitude is 104NT magnitude, and become caused by seismic precursor
Change is mainly manifested in lithospheric magnetic field, and only several nT are even more small, very faint for overall amplitude, directly to original
Beginning measurement magnetic field is operated, and exception can not be extracted.Therefore the variable quantity of magnetic field data is obtained by first-order difference first, due to
There are cyclically-varying, the result after difference still exists relative to the biggish low frequency of anomalous variation in earth's magnetic field of measurement itself
Ingredient.Differentiated data are decomposed using wavelet transform, the period is eliminated by the low-frequency component in removal data
Property variation influence, finally obtain the situation of change in magnetic field, result can be used for the subsequent operations such as principal component analysis.
Wavelet transform obtains the formula of every layer of radio-frequency component and low-frequency component are as follows:
Wherein, j is Decomposition order, and L indicates that low-frequency component, H indicate that radio-frequency component, g are low-pass filter, and h is high pass filter
Wave device, dB are magnetic field differential data.The residual error of every track magnetic field three-component differential data and low-frequency component is pre-processed results.
Step d is handled using principal component analysis satellite magnetic field Data Dimensionality Reduction, is the principal component by seeking three component seismic data
Dimension-reduction treatment is carried out to data, and retains the main feature of multidimensional data as far as possible.The changes of magnetic field that Earthquake Precursor Anomalies generate
There should be reflection on the different components in magnetic field, dimension-reduction treatment is carried out to magnetic field three component seismic data using principal component analysis, it can
It is simplified the calculating process of subsequent processing to remove redundancy, while retaining most important feature in three-component magnetic field.
Dimensionality reduction result contains the main feature of three component datas simultaneously, analyzes it processing than individually carrying out to some component
Processing result more comprehensively, can preferably react information of earthquake.It specifically includes:
The magnetic field three component seismic data of every track after pretreatment is expressed as according to time series:
Xi=[xi(1),xi(2)...,xi(m)], i=1,2 ..., n (2)
Wherein m is track length, and n is data dimension.Matrix then can be obtainedThe expression formula of Y
Are as follows:
The covariance matrix C of calculating matrix YYThe element γ of (n × n)uv, calculation formula are as follows:
Wherein, xiuAnd xivThe i-th row u column element and the i-th row v column element in respectively matrix Y；WithIt is respectively
The mean value of u column and v column element.
Calculate the characteristic value and feature vector of covariance matrix:
CY=R Λ RT (5)
Wherein, Λ (n × n) is the characteristic value diagonal matrix arranged from big to small, and R (1 × n) is the feature of corresponding eigenvalue
The list of feature values is shown as λ by vector1,λ2,...,λn(λ1> λ2> ... > λn)。
Principal component Φ is can be obtained into initial data progress linear projection using matrix R.
Φ=RY=[Φ1,Φ2,...,Φn]T (6)
Wherein, Φ1,Φ2,...,ΦnIt is the 1st to the n principal component (1 × m), the contribution rate of preceding k principal component is right by its
Characteristic value is answered to calculate:
The preceding k principal component that choosing makes contribution rate to 60% or more represents the main feature of magnetic field data, and as drop
Tie up the result of processing.
Step e carries out continuous wavelet transform to satellite magnetic field dimensionality reduction data, and seeks track using obtained wavelet coefficient
Energy intensity extracts abnormal track.The result that dimension-reduction treatment obtains contains three-component data information, but wherein
Exception be not obvious, directly by time domain progress anomaly extracting have difficulties.Therefore by continuous wavelet transform by signal by when
Domain transforms to time-frequency domain, calculates the energy intensity at signal each moment, its prominent anomalous variation in time-frequency domain, and use sliding square
The equalization that shape window carries out data excludes mutation sexual abnormality, finally utilizes this rail to ensure the aberrant continuation regular hour
Maximum average energy intensity on road reflects its anomalous variation situation, and ideal effect can be obtained by carrying out anomaly extracting at this time
Fruit.Pass through continuous wavelet transform and calculate orbital energy intensity, very effective can extract the abnormal rail in dimensionality reduction data
Road.
The formula of continuous wavelet transform are as follows:
Wherein a is scale factor, and b is shift factor, and Ψ is mother wavelet function, and f (t) is pending data.Every track
The continuous wavelet transform result of magnetic field dimensionality reduction data is a matrix:
Wherein m is the length of track, and n is determined by the scale factor a chosen.
The wavelet coefficient meter of time-frequency domain can be can use from the energy of frequency-domain calculations signal according to principle of conservation of energy
The energy at energy intensity characterization magnetic field dimensionality reduction data each moment is calculated, to protrude the anomalous variation in magnetic field.It is calculated by wavelet coefficient
The formula of energy intensity are as follows:
Then obtain the energy intensity data sequence of every track are as follows: Ze=[Ze(1),Ze(2),...,Ze(m)]。
To ensure the aberrant continuation regular hour, mutation sexual abnormality is excluded, using a length of l of window, step-length is that the rectangular window of s exists
It is slided on track, seeks the mean value of each window self-energy intensity, and being rounded a maximum value for all mean values of track is this track
Energy intensity Zemax, the energy intensity that this track has is represented, the sequence of all orbital energy intensity is obtained are as follows:
Zemax=[Zemax(1),Zemax(2),...,Zemax(p)] (10)
Wherein p is magnetostatic track number.
Calculate all orbital energy intensity ZemaxMean μZAnd standard deviation sigmaZ, and threshold value is set as ThZ=μZ+kk×σZ,
Middle kk is empirical coefficient, and when the energy intensity of track is greater than the threshold value of setting, which is considered as abnormal track.
Step f, Ionospheric measurement parameter are affected by various factors, such as season, geographical location, thermal testing, ionization
Layer disturbance and terrestrial magnetic disturbance.To further confirm that extract exception is caused by earthquake, must exclude different caused by other factors influence
Often.Step b removes the track influenced by terrestrial magnetic disturbance by geomagnetic index, since the influence that non-shake factor generates belongs to entirely
Ball interference, therefore occur in earthquake effect region but also the exception outside earthquake effect region by removing on same track not only,
And it appears in earthquake effect region but extends to the exception outside influence area or only occur in outside earthquake effect region
Abnormal, excluding non-shake factor influences caused exception.
Step g, the exclusion influenced by non-shake factor, using finally obtained abnormal track as Earthquake Precursor Anomalies track
It is exported.
Below with reference to embodiment, the present invention is described in detail.
For 7.8 grades of Ecuador's earthquakes that on April 16th, 2016 occurs, with the magnetic field three of Swarm satellite group Alpha star
For component data (1Hz).According to Dobrovolsky's formula R=100.43M(earthquake magnitude that M is earthquake) selection Ecuador
The survey region of shake by the square area of side length of 4518.8km is survey region centered on earthquake centre, chooses first 90 days of shake
30 days after to shake, i.e., on January 17th, 2016 to May 16 is the time range studied.Ecuador earthquake epicenter, earthquake effect
As shown in Fig. 2, wherein five-pointed star is earthquake centre, border circular areas is earthquake effect region, square area for region and survey region
For survey region.
A, typing Swarm satellite group Alpha star on January 17th, 2016 to magnetic field three component seismic data on May 16 (1Hz), counts
It is denoted as B respectively according to by north orientation, east orientation and vertical componentx、ByAnd Bz.The track in region after study is picked out, and according in data
Flag bit remove invalid data, be used only effective magnetic field three component seismic data.Since the time of measurement data and space can
It changes with satellite flight, therefore processing and anomaly extracting is carried out to data according to track.It is on the April 1st, 2016 of track 1
Example, the three-component primitive curve in magnetic field are as shown in Figure 3, it can be seen that magnetic field initial data amplitude is very big, can not observe that its is micro-
Weak variation, exception can not be extracted by directly operating on it.
B, the geomagnetic index ap and Dst in research on utilization time range, according to magnetostatic condition: current time ap < 12, |
Dst|≤20；Preceding 23 hour ap≤32, | Dst |≤30.Removal is interfered biggish track by geomagnetic activity, avoids extracting different
Often caused by geomagnetic activity, obtain it is magnetostatic under the conditions of 226 tracks magnetic field three component seismic data.
C, first-order difference is carried out respectively by track to the magnetic field three component seismic data of selection to handle to obtain changes of magnetic field amount, then divide
It is other that wavelet transform is carried out to differential data, remove low-frequency component therein, the situation of change of prominent data.Discrete wavelet transformer
Change the formula for obtaining every layer of radio-frequency component and low-frequency component are as follows:
The characteristics of for magnetic field data itself, taking Decomposition order j is 1 to 3, and selects the low-pass filter of db4 wavelet basis
Indicate that low-frequency component, H indicate radio-frequency component as g and h, L with high-pass filter, dB is magnetic field differential data.It will be calculated
3rd layer of low-frequency component dB3,LAs removal item, every three-component differential data in track magnetic field is residual with the 3rd layer of low-frequency component
Difference is pre-processed results.
By taking on the April 1st, 2016 of track 1 as an example, differential data and low-frequency component are respectively such as the solid line and dotted line in Fig. 4 a
It is shown, it can be seen that still there are the biggish low-frequency cycle ingredient of amplitude in differential data, need further to be removed, it is discrete small
The low-frequency component that wave conversion calculates is corresponding with differential data good, will not generate after subtracting each other pseudo- abnormal.Utilize wavelet transform
Removal low-frequency component can effectively remove magnetostatic field background, the situation of change in prominent magnetic field itself.Differential data subtract low frequency at
The residual error curve divided is as shown in Figure 4 b, and step processing can be very good to obtain the situation of change in magnetic field known to result figure.
D, the magnetic field data of every track after pretreatment is expressed as in temporal sequence
Xi=[xi(1),xi(2)...,xi(m)], i=1,2 ..., n (2)
Wherein m is track length, and n is that data latitude takes 3.Matrix then can be obtainedThe expression formula of Y
Are as follows:
The covariance matrix C of calculating matrix YYThe element γ of (3 × 3)uv, calculation formula are as follows:
Wherein, xiuAnd xivThe u column and v column data of i-th row in respectively matrix Y；WithBe respectively u column and
The mean value of v column data.
Calculate the characteristic value and feature vector of covariance matrix:
CY=R Λ RT (5)
In formula, Λ (3 × 3) is the characteristic value diagonal matrix arranged from big to small, and R (1 × 3) is the feature of corresponding eigenvalue
The list of feature values is shown as λ by vector1,λ2,λ3(λ1> λ2> λ3)。
Initial data, which is carried out linear projection, with matrix R can be obtained principal component Φ.
Φ=RY=[Φ1,Φ2,Φ3]T (6)
Wherein, Φ1,Φ2,Φ3For the 1st to the 3rd principal component (1 × m), the contribution rate of preceding k principal component is corresponded to special by it
Value indicative calculates:
It is computed the contribution rate ζ of discovery first principal component1> 60%, therefore choose first principal component Φ1Represent magnetic field data
Main feature, and the result as dimension-reduction treatment.
By taking on the April 1st, 2016 of track 1 as an example, dimensionality reduction result is as shown in figure 5, with Fig. 4 b comparison as can be seen that principal component
Data after dimensionality reduction remain most important feature in signal, and the information for eliminating redundancy obtains the calculating process of subsequent processing
To simplification, be conducive to subsequent abnormal identification and extraction.
E, continuous wavelet transform is carried out to magnetic field dimensionality reduction data, signal is transformed into time-frequency domain by time domain.Continuous wavelet becomes
The formula changed are as follows:
According to the characteristic of magnetic field dimensionality reduction data, it is db4 wavelet basis function that enable scale factor a, which be 512, Ψ, and f (t) is Φ1。
The continuous wavelet transform result of every track dimensionality reduction data is a matrix:
Wherein m is the length of track, and n is determined as 512 by scale factor a.By taking on the April 1st, 2016 of track 1 as an example, magnetic field
The wavelet transform result of dimensionality reduction data is as shown in Figure 6 a, and the available orbital data is used in the Energy distribution situation of time-frequency domain
In the energy for seeking energy intensity and characterizing each moment.
The wavelet coefficient of time-frequency domain can be utilized to calculate energy from the energy of frequency-domain calculations signal according to principle of conservation of energy
The energy at intensity characterization magnetic field dimensionality reduction data each moment is measured, to protrude the anomalous variation in magnetic field.Energy is calculated by wavelet coefficient
The formula of intensity are as follows:
Then obtain the energy intensity data sequence of every track are as follows: Ze=[Ze(1),Ze(2),...,Ze(m)]。
To ensure the aberrant continuation regular hour, mutation sexual abnormality is excluded, using a length of l=20 of window, step-length is s=20's
Rectangular window slides in orbit, seeks the mean value of each window self-energy intensity.By taking on the April 1st, 2016 of track 1 as an example, it is averaged
Energy intensity curve is as shown in Figure 6 b, and the dimensionality reduction data of comparison diagram 5 can protrude abnormal change it is found that seeking the energy intensity of data
Change, seeking mean value can reflect mean intensity of the magnetic field energy within this time, eliminate the influence of Sudden Anomalies.
It is rounded the energy intensity Z that a maximum value for all mean values of track is this trackemax, represent what this track had
Energy intensity obtains the sequence of all orbital energy intensity are as follows:
Zemax=[Zemax(1),Zemax(2),...,Zemax(226)] (10)
It is 226 that wherein p, which is magnetostatic track number,.
Calculate all orbital energy intensity ZemaxMean μZAnd standard deviation sigmaZ, kk=3 is taken to set threshold value as ThZ=μZ+3×
σZ, when the characteristic value of track is greater than the threshold value of setting, which is considered as abnormal track.Anomaly extracting result such as Fig. 7 institute
Show, wherein dash-dotted horizontal line is respectively mean μ from bottom to top on the day of vertical dotted line is earthquakeZWith threshold value ThZ=μZ+3×σZ, star
Number be the abnormal track beyond threshold value.Track before shaking in figure beyond threshold value is respectively on January 25th, 2016 track 63,2 months
29 days 3,2 months tracks 3 of track and track 1 on March 14.It can be seen that from result figure and pass through continuous wavelet transform and calculate orbital energy
Intensity can protrude exception, effectively extract the abnormal track in all magnetostatic tracks.
F, Ionospheric measurement parameter are affected by various factors, such as season, geographical location, thermal testing, ionosphere disturb
Dynamic and terrestrial magnetic disturbance.To further confirm that extract exception is caused by earthquake, must exclude exception caused by other factors influence.Step
Rapid b removes the track influenced by terrestrial magnetic disturbance by geomagnetic index, due to the influence that non-shake factor generates belong to it is global
Interference, therefore by removing exception 2016 not only occurred in earthquake effect region again outside earthquake effect region on same track
Track 3, and appeared in exception that is in earthquake effect region but extending to outside influence area or only occurred in earthquake January 25
Abnormal on 2 6th, 2016 tracks 3 outside influence area, excluding non-shake factor influences caused exception.
G, the exclusion influenced by non-shake factor, finally only on 2 29th, 2016 track 3 and 2016 rails on March 14,
Road 1 is exception, and as shown in figs. 8 a and 8b, the treatment process known to extraction result through the invention can be in degaussing field data
Amplitude is big and the slow magnetostatic field of variation obtains the situation of change of magnetic field data, and makes full use of the magnetic field three component seismic data to include
Information effectively protrudes anomalous variation, guarantees to extract the abnormal duration, excludes mutation sexual abnormality, obtains good exception
Extract result.29 days 2 months tracks 3 and track 1 on March 14 are exported as Ecuador's Earthquake Precursor Anomalies track.
The foregoing is merely illustrative of the preferred embodiments of the present invention, is not intended to limit the invention, all in essence of the invention
Made any modifications, equivalent replacements, and improvements etc., should all be included in the protection scope of the present invention within mind and principle.
Claims (8)
1. a kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation, which is characterized in that this method packet
It includes:
A, the magnetic field data of satellite is read, and valid data are chosen according to flag bit；
B, the tranquiler magnetostatic track of geomagnetic activity is chosen according to geomagnetic index；
C, the magnetic field three component seismic data of magnetostatic track is pre-processed by difference and wavelet transform, removal amplitude greatly and
Change slow magnetostatic field part, obtains the situation of change of magnetic field data；
D, dimension-reduction treatment is carried out to pretreated magnetic field three component seismic data using principal component analysis, while removing redundancy
Main feature in stick signal；
E, continuous wavelet transform is carried out to the magnetic field data after dimensionality reduction, orbital energy intensity is defined by wavelet coefficient, and utilize
It carries out abnormal track and extracts；
F, excluding non-seismic factor influences caused exception；
G, the extraction result of Earthquake Precursor Anomalies is exported.
2. according to the method for claim 1, which is characterized in that step a includes: the magnetic field number read in moonscope result
According to magnetic field three-component Bx、ByAnd Bz, and invalid data is removed according to the flag bit of data.
3. according to the method for claim 1, which is characterized in that it includes: logical that step b, which chooses magnetostatic track according to geomagnetic index,
It crosses geomagnetic index ap and Dst removal and is interfered biggish track, the condition that magnetostatic track meets are as follows: current time ap by geomagnetic activity
< 12, | Dst |≤20；Preceding 23 hour ap≤32, | Dst |≤30.
4. according to the method for claim 1, which is characterized in that step c include first to the magnetic field three component seismic data of selection by
Track carries out first-order difference respectively and handles to obtain changes of magnetic field amount, then carries out discrete wavelet transformer to data after difference processing respectively
It changes, removes low-frequency component therein.
5. according to the method for claim 4, which is characterized in that wavelet transform obtain every layer of radio-frequency component and low frequency at
The formula divided are as follows:
Wherein, j is Decomposition order, and L indicates that low-frequency component, H indicate radio-frequency component, and g is low-pass filter, and h is high-pass filter,
DB is magnetic field differential data, and the residual error of every track magnetic field three-component differential data and low-frequency component is pre-processed results.
6. according to the method for claim 1, which is characterized in that step d, described utilizes principal component analysis satellite magnetic field number
Include: according to dimension-reduction treatment
The magnetic field three component seismic data of every track after pretreatment is expressed as according to time series:
Xi=[xi(1),xi(2)...,xi(m)], i=1,2 ..., n
Wherein m is track length, and n is data dimension, obtains matrixThe expression formula of Y are as follows:
The covariance matrix C of calculating matrix YYThe element γ of (n × n)uv, calculation formula are as follows:
Wherein, xiuAnd xivThe i-th row u column element and the i-th row v column element in respectively matrix Y；WithIt is u column respectively
With the mean value of v column element；
Calculate the characteristic value and feature vector of covariance matrix:
CY=R Λ RTWherein, Λ (n × n) is the characteristic value diagonal matrix arranged from big to small, and R (1 × n) is corresponding eigenvalue
The list of feature values is shown as λ by feature vector1,λ2,...,λn(λ1> λ2> ... > λn)；
The initial data of matrix Y is subjected to linear projection using matrix R and obtains principal component Φ,
Φ=RY=[Φ1,Φ2,...,Φn]T
Wherein, Φ1,Φ2,...,ΦnIt is the 1st to the n principal component (1 × m), the contribution rate of preceding k principal component is corresponded to special by it
Value indicative calculates:
The preceding k principal component that choosing makes contribution rate to 60% or more represents the main feature of magnetic field data, and as dimensionality reduction at
The result of reason.
7. according to the method for claim 1, which is characterized in that step e carries out continuous wavelet to the magnetic field data after dimensionality reduction
Transformation defines orbital energy intensity by wavelet coefficient, and carries out abnormal track extraction using it and include:
Signal is transformed into time-frequency domain by time domain by continuous wavelet transform；
The energy intensity at signal each moment is calculated in time-frequency domain；
Data are averaged using sliding rectangular window, to ensure the aberrant continuation regular hour, exclude mutation sexual abnormality；
Reflect its anomalous variation situation using the maximum average energy intensity on this track.
8. according to the method for claim 7, which is characterized in that further include: calculate the mean μ of all orbital energy intensityZWith
Standard deviation sigmaZ, and threshold value is set as ThZ=μZ+kk×σZ, when the energy intensity of track is greater than the threshold value of setting, which is recognized
To be abnormal track.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811554885.4A CN109740453B (en) | 2018-12-19 | 2018-12-19 | Satellite magnetic field data earthquake precursor anomaly extraction method based on wavelet transformation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811554885.4A CN109740453B (en) | 2018-12-19 | 2018-12-19 | Satellite magnetic field data earthquake precursor anomaly extraction method based on wavelet transformation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109740453A true CN109740453A (en) | 2019-05-10 |
| CN109740453B CN109740453B (en) | 2022-03-29 |
Family
ID=66360608
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811554885.4A Active CN109740453B (en) | 2018-12-19 | 2018-12-19 | Satellite magnetic field data earthquake precursor anomaly extraction method based on wavelet transformation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109740453B (en) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110673206A (en) * | 2019-08-26 | 2020-01-10 | 吉林大学 | Satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization |
| CN112327371A (en) * | 2020-11-06 | 2021-02-05 | 吉林大学 | Satellite magnetic field data time-varying background field establishment method based on variational modal decomposition |
| CN112401906A (en) * | 2020-11-10 | 2021-02-26 | 河北省科学院应用数学研究所 | Automatic electroencephalogram artifact removing method based on amplitude |
| CN113435259A (en) * | 2021-06-07 | 2021-09-24 | 吉林大学 | Tensor decomposition-based satellite magnetic field data fusion seismic anomaly extraction method |
| CN113568057A (en) * | 2021-07-30 | 2021-10-29 | 应急管理部国家自然灾害防治研究院 | Zhang Heng I satellite induction type magnetometer data processing method and system |
| CN117665949A (en) * | 2023-12-04 | 2024-03-08 | 中国科学院地质与地球物理研究所 | Fluxgate aeromagnetic measurement method and system |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EA200801036A1 (en) * | 2008-05-06 | 2009-12-30 | Леонид Николаевич Дода | METHOD FOR PREDICTING EARTHQUAKES |
| CN106021710A (en) * | 2016-05-13 | 2016-10-12 | 南京航空航天大学 | Seismic precursor satellite orbit anomaly identification method based on atmosphere ionosphere parameter |
| CN106918836A (en) * | 2017-03-31 | 2017-07-04 | 吉林大学 | Borehole strain data exception extraction method based on principal component analysis |
| CN107037486A (en) * | 2017-03-31 | 2017-08-11 | 中国地质大学（武汉） | The Time-frequency Spectrum Analysis method and system of earth natural pulses electromagnetic field data processing |
| CN107356969A (en) * | 2017-09-06 | 2017-11-17 | 四川易利数字城市科技有限公司 | A kind of seismic precursor analysis method based on satellite thermal infrared data and GIS |
| CN107422381A (en) * | 2017-09-18 | 2017-12-01 | 西南石油大学 | A kind of earthquake low-frequency information fluid prediction method based on EEMD ICA |
| CN109031403A (en) * | 2018-08-20 | 2018-12-18 | 吉林大学 | A kind of borehole strain data exception extraction method based on S-K feature |
-
2018
- 2018-12-19 CN CN201811554885.4A patent/CN109740453B/en active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EA200801036A1 (en) * | 2008-05-06 | 2009-12-30 | Леонид Николаевич Дода | METHOD FOR PREDICTING EARTHQUAKES |
| CN106021710A (en) * | 2016-05-13 | 2016-10-12 | 南京航空航天大学 | Seismic precursor satellite orbit anomaly identification method based on atmosphere ionosphere parameter |
| CN106918836A (en) * | 2017-03-31 | 2017-07-04 | 吉林大学 | Borehole strain data exception extraction method based on principal component analysis |
| CN107037486A (en) * | 2017-03-31 | 2017-08-11 | 中国地质大学（武汉） | The Time-frequency Spectrum Analysis method and system of earth natural pulses electromagnetic field data processing |
| CN107356969A (en) * | 2017-09-06 | 2017-11-17 | 四川易利数字城市科技有限公司 | A kind of seismic precursor analysis method based on satellite thermal infrared data and GIS |
| CN107422381A (en) * | 2017-09-18 | 2017-12-01 | 西南石油大学 | A kind of earthquake low-frequency information fluid prediction method based on EEMD ICA |
| CN109031403A (en) * | 2018-08-20 | 2018-12-18 | 吉林大学 | A kind of borehole strain data exception extraction method based on S-K feature |
Non-Patent Citations (3)
| Title |
|---|
| P.HAN等: "Evaluation of ULF electromagnetic phenomena associated with the 2000 Izu Islands earthquake swarm by wavelet transform analysis", 《NATURAL HAZARDS AND EARTH SYSTEM SCIENCES》 * |
| QIU ZEHUA等: "Abnormal strain changes observed by a borehole strainmeter at Guza Station before the Ms7 0 Lushan earthquake", 《GEODESY AND GEODYOAMICS》 * |
| 池成全等: "基于主成分分析的芦山Ms7.0地震前钻孔应变数据异常提取", 《2017中国地球科学联合学术年会论文集（十五）》 * |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110673206A (en) * | 2019-08-26 | 2020-01-10 | 吉林大学 | Satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization |
| CN110673206B (en) * | 2019-08-26 | 2020-12-29 | 吉林大学 | Satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization |
| CN112327371A (en) * | 2020-11-06 | 2021-02-05 | 吉林大学 | Satellite magnetic field data time-varying background field establishment method based on variational modal decomposition |
| CN112327371B (en) * | 2020-11-06 | 2021-07-30 | 吉林大学 | Satellite magnetic field data time-varying background field establishment method based on variational modal decomposition |
| CN112401906A (en) * | 2020-11-10 | 2021-02-26 | 河北省科学院应用数学研究所 | Automatic electroencephalogram artifact removing method based on amplitude |
| CN112401906B (en) * | 2020-11-10 | 2021-12-14 | 河北省科学院应用数学研究所 | Automatic electroencephalogram artifact removing method based on amplitude |
| CN113435259A (en) * | 2021-06-07 | 2021-09-24 | 吉林大学 | Tensor decomposition-based satellite magnetic field data fusion seismic anomaly extraction method |
| CN113435259B (en) * | 2021-06-07 | 2022-06-03 | 吉林大学 | Tensor decomposition-based satellite magnetic field data fusion earthquake anomaly extraction method |
| CN113568057A (en) * | 2021-07-30 | 2021-10-29 | 应急管理部国家自然灾害防治研究院 | Zhang Heng I satellite induction type magnetometer data processing method and system |
| CN117665949A (en) * | 2023-12-04 | 2024-03-08 | 中国科学院地质与地球物理研究所 | Fluxgate aeromagnetic measurement method and system |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109740453B (en) | 2022-03-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109740453A (en) | A kind of satellite magnetic field data Earthquake Precursor Anomalies extracting method based on wavelet transformation | |
| Yang et al. | Unsupervised 3-D random noise attenuation using deep skip autoencoder | |
| Saad et al. | Scalodeep: A highly generalized deep learning framework for real‐time earthquake detection | |
| CN106021710B (en) | Satellite orbit abnormality recognition method before shake based on atmospheric ionization layer parameter | |
| Liebrand et al. | Antarctic ice sheet and oceanographic response to eccentricity forcing during the early Miocene | |
| Li et al. | Identifying P-phase arrivals with noise: An improved Kurtosis method based on DWT and STA/LTA | |
| CN110673206B (en) | Satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization | |
| Stolz et al. | Superconducting sensors and methods in geophysical applications | |
| Li et al. | Novel wavelet threshold denoising method to highlight the first break of noisy microseismic recordings | |
| Akhoondzadeh | Support vector machines for TEC seismo-ionospheric anomalies detection | |
| CN110068857B (en) | Swarm double-star magnetic field data earthquake precursor anomaly extraction method based on principal component analysis | |
| Hafez et al. | Seismic noise study for accurate P-wave arrival detection via MODWT | |
| Jiang et al. | Comparison of the earthquake detection abilities of PhaseNet and EQTransformer with the Yangbi and Maduo earthquakes | |
| Francis et al. | Nonlinear prediction of the ionospheric parameter ƒ o F 2 on hourly, daily, and monthly timescales | |
| US4449099A (en) | Recording decimated magnetotelluric measurements with coherence | |
| Xue et al. | LSTM-autoencoder network for the detection of seismic electric signals | |
| CN113158830A (en) | Residual gravity abnormal field separation method | |
| Langel | The use of low altitude satellite data bases for modeling of core and crustal fields and the separation of external and internal fields | |
| Lan et al. | An analysis of mechanical constraints when using superconducting gravimeters for far-field pre-seismic anomaly detection | |
| CN112327371B (en) | Satellite magnetic field data time-varying background field establishment method based on variational modal decomposition | |
| Du et al. | Magnetic anomaly detection based on singular spectrum analysis and orthonormal basis functions | |
| CN113435259A (en) | Tensor decomposition-based satellite magnetic field data fusion seismic anomaly extraction method | |
| Vorobev et al. | Development and application of problem-oriented digital twins for magnetic observatories and variation stations | |
| Love et al. | Extreme-event geoelectric hazard maps | |
| Yu et al. | A Novel Method based on Proximate Wavelet Coefficient Recovery for Magnetic Resonance Sounding Signal Denoising in Complex Interference Environments |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |