CN109470135A - CSAMT data inactivity bearing calibration - Google Patents
CSAMT data inactivity bearing calibration Download PDFInfo
- Publication number
- CN109470135A CN109470135A CN201811337149.3A CN201811337149A CN109470135A CN 109470135 A CN109470135 A CN 109470135A CN 201811337149 A CN201811337149 A CN 201811337149A CN 109470135 A CN109470135 A CN 109470135A
- Authority
- CN
- China
- Prior art keywords
- point
- signal
- imf component
- value
- apparent resistivity
- 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.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
- G01B7/26—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring depth
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The present invention relates to a kind of CSAMT data inactivity bearing calibrations, method includes the following steps: fitting controlled-source audiomagnetotelluric original signal according to controlled-source audiomagnetotelluric discrete point apparent resistivity data, each IMF component of empirical mode decomposition is carried out to it；Static relays are identified according to the cross-correlation coefficient between the apparent resistivity value of whole frequency points of observation point and reference point, if there is static relays, then original signal is reconstructed according to each IMF component, the controlled-source audiomagnetotelluric reconstruction signal s'(t after obtaining static shift correction).The present invention can retain the information of anomalous body in initial data, keep its undistorted while suppressing static relays.
Description
Technical field
The present invention relates to a kind of static calibration methods for CSAMT data.
Background technique
Controlled-source audiomagnetotellurics method (Controlled Source Audio-frequency
Magnetotellurics, abbreviation CSAMT Carnot model sounding curve, therefore also known as controlled-source audiomagnetotellurics sounding
Method, it is excellent to have that work efficiency is high, depth of exploration range is big, spatial resolution is high, the influence of topography is small, high resistant electric screening action is small etc.
Point.CSAMT is in signal acquisition process, and earth's surface or near surface often have the conductive heterogeneous body in part, when electric current flows through unevenness
It when the surface of even body, just will form charge accumulation effects, electric field caused to be distorted, generate static relays.In " apparent resistivity-
The apparent resistivity curve of adjacent each measuring point is moved up and down along apparent resistivity axis in frequency " log-log coordinate system, is reflected in view electricity
On resistance rate pseudosection map, it may appear that the little shape dense band that rises steeply of lateral extent covers underground as caused by ore body or object
It is abnormal, the judgement of practical geological condition is caused to mislead.
Can be divided into the bearing calibration of static relays: (1) curvilinear translation method does static correction；(2) numerical analysis method does quiet
Correction；(3) joint inversion method does static correction；(4) static correction is done using phase conversion information；(5) quiet school is done using magnetic field data
Just；(6) direct two dimension and multidimensional inversion method；(7) wavelet filteration method etc..The quality of curvilinear translation method effect is very big
Treatment people is depended in degree to the human-subject test of data, judgement and process experience to static relays play very crucial work
With being easy to miss by those is not caused by static relays extremely as static relays to " correction "；For certain terrain backgrounds
The static relays that simple static relays, the especially influence of topography generate, can be identified using the method for numerical simulation
And correction is made, and these are often unknown under practical field condition, so can only study to use as reference；Joint is anti-
The method of drilling is equivalent to repetition and has done two kinds of geophysical methods, at high cost, not very practical；Wavelet transformation is a kind of change window change
It changes, static relays caused by earth's surface inhomogeneities is mathematically intuitively distinguished extremely with macrotectonics, there is higher point
Resolution and correction accuracy, however wavelet method there are the problem of be exactly wavelet basis function selection, Selection of Wavelet Basis it is different with
And the difference of Decomposition order all can cause decomposition result larger difference occur, and then will affect subsequent static shift correction effect.
It is American National aerospace that Hilbert-Huang, which converts (Hilbert-huang Transform, abbreviation HHT) method,
Non-stationary, nonlinear properties are had clearly physics by the new signal processing method of one kind that N.E.Huang of office et al. is proposed
Meaning can obtain one time of amplitude, the one frequency distribution feature of signal, and have adaptivity.It is being pressed for conventional method
The problems of when static relays processed, the invention proposes the static shift corrections of the CSAMT data based on HHT.
Summary of the invention
The technical problem to be solved in the present invention is to provide a kind of CSAMT data inactivity bearing calibration, base is not present in this method
The problem of function selects can retain macrotectonics while suppressing static relays while realizing static shift correction and believe extremely
Breath.
In order to solve the above-mentioned technical problem, CSAMT data inactivity bearing calibration of the invention the following steps are included:
Step 1: acquisition controlled-source audiomagnetotelluric discrete point apparent resistivity data, fits controllable source audio the earth
Electromagnetism original signal s (t) carries out empirical mode decomposition (Empirical Mode Decomposition, abbreviation EMD) to it and obtains
To n intrinsic mode function c1(t)、c2(t)、...ci(t)、...cn(t), and using each intrinsic mode function as IMF component；
Step 2: according to the cross-correlation coefficient r between the apparent resistivity value of whole frequency points of observation point and reference pointxyIdentification
Static relays are then reconstructed original signal s (t) according to step 3 if there is static relays；
Step 3: original signal s (t) is reconstructed according to each IMF component that step 1 obtains, after obtaining static shift correction
Controlled-source audiomagnetotelluric reconstruction signal s'(t)；
The method that original signal s (t) is reconstructed is as follows:
(1) using the adaptive threshold based on white Gaussian noise as reconstruct threshold value, if i-th of IMF component ci(t) weight
Structure threshold value is Ti:
Wherein, c indicates threshold coefficient, value range 0.5-0.7；M indicates the signal length chosen, i.e., each IMF component
Discrete point number, EiIndicate i-th of IMF component ci(t) energy intensity, its calculation formula is:
Wherein,It is first IMF component c1(t) variance；
(2) h is used(ik)(t) i-th of IMF component c is indicatedi(t) in time tkThe discrete value at place, k=1,2 ... M, then i-th
A IMF component ci(t) corresponding reconstruct sub-signal is
Then controlled-source audiomagnetotelluric reconstruction signal is s'(t):
The preferred c=0.7 of the threshold coefficient.
In the step one, empirical mode decomposition is carried out to the apparent resistivity data of CSAMT and obtains multiple intrinsic mode
The method of function component is as follows:
(1), controlled-source audiomagnetotelluric original signal s (t) all maximum point and minimum point are found out；
(2), it connects all maximum point and minimum point to form upper packet respectively using cubic spline interpolation
Winding thread S1, lower envelope line S2: coenvelope line S is subtracted with each discrete point apparent resistivity value ρ of original signal s (t)1With lower envelope line S2
Corresponding points mean value, obtain multiple discrete difference points and it be fitted to obtain first difference matched curve h11(t)；By
One difference matched curve h11(t) all maximum points and minimum point connect to form difference signal coenvelope line H11、
Lower envelope line H12；Investigate first difference matched curve h11(t) whether meet two conditions of IMF, that is: extreme point number
(including maximum and minimum) and zero crossing number it is equal or it is most difference 1；The corresponding difference signal coenvelope in arbitrary point
Line H11On discrete point value and difference signal lower envelope line H12On discrete point value average value be 0；By first if meeting
A difference matched curve h11(t) it is used as the 1st intrinsic mode function c1(t), step 3 is gone to；Otherwise quasi- with first difference
Close curve h11(t) substitution original signal s (t) repeats above-mentioned screening process, meets IMF until filtering out from original signal s (t)
Two conditions difference matched curve and as the 1st intrinsic mode function c1(t)；
(3) by remaining signal r1(t)=s (t)-c1(t) is obtained according to step (1), (two) as new signal
Two intrinsic mode function c2(t)；And so on, finally obtain n intrinsic mode function c1(t)、c2(t)、...ci(t)、
...cn(t) and a nondecomposable residual components rn(t):
N intrinsic mode function c1(t)、c2(t)、...ci(t)、...cn(t) it is used as IMF component.
In the step 2, identification static relays method is as follows: the apparent resistivity value of whole frequency points of observation point is regarded as
It is one group of data series xi, apparent resistivity value corresponding to whole frequency points with reference also to point is also one group of data series yi, by this
Two groups of data carry out relevant matches, seek cross-correlation coefficient r between the twoxy,
Wherein N is frequency point quantity,For the average value of whole frequency point apparent resistivity values of observation point,The whole of reference point
The average value of frequency point apparent resistivity value；
If rxy>=0.85 is thought that there are static relays.
The present invention can also include the following steps:
Step 4: each IMF component c obtained to step 11(t)、c2(t)、...ci(t)、...cn(t) make Hilbert change
Get its instantaneous frequency and Hilbert spectrum and marginal spectrum in return；For any IMF component ci(t), instantaneous frequency ωi(τ),
Hilbert spectrum H (ω, τ) of original signal s (t) and the marginal spectrum h (ω) of original signal s (t) are as follows:
In formula, yi(τ) is IMF component ci(t) signal after Hilbert transform, τ are signals after Hilbert transform
Time variable, θi(τ) is phase,
H(ωi, τ) and=H (ωi(τ),τ)≡ai(τ)
Wherein H (ωi, τ) and it is IMF component ci(t) spectral line in corresponding Hilbert spectrum, aiVariable Amplitude when (τ) is；
The beneficial effects of the present invention are: (1) traditional static bearing calibration (such as: wavelet method) all there is basic function selection
Difficult problem, and this problem is not present in method proposed by the present invention, avoids this restrictive condition.(2) present invention proposes
CSAMT static shift correction can while suppressing static relays, retain initial data in anomalous body information, lose it not
Very.(3) the Hilbert spectrum and marginal spectrum h (ω) and the HHT of initial data spectrum and HHT marginal spectrum obtained according to step 4 carries out
It compares, can intuitively reflect the correction situation of static relays.The difference for showing signal day part can be understood in time-frequency spectrum
Frequency, it can be seen that dominant frequency range present in signal from marginal spectrum.
Detailed description of the invention:
Fig. 1: CSAMT data inactivity bearing calibration flow chart of the invention.
Fig. 2: EMD decomposes the process schematic of controlled source audio-frequency magnetotelluric data screening IMF component.
Fig. 3: the HHT spectrum of controlled source audio-frequency magnetotelluric data.
Fig. 4: the HHT marginal spectrum of controlled source audio-frequency magnetotelluric data.
Fig. 5: homogeneous half space forward model.
Fig. 6 (a), Fig. 6 (b), Fig. 6 (c) are respectively that homogeneous half space forward model is handled without static shift correction, based on small
The processing of wave method static shift correction and the apparent resistivity isoline figure that static shift correction is carried out using the present invention.
Fig. 7: the homogeneous half space forward model containing cap rock.
Fig. 8 (a), Fig. 8 (b), Fig. 8 (c) be respectively the homogeneous half space forward model containing cap rock without static shift correction at
Reason, the apparent resistivity isoline figure that static shift correction is carried out based on the processing of wavelet method static shift correction and using the present invention.
Specific embodiment:
Invention is further described in detail with embodiment with reference to the accompanying drawing.
As shown in Figure 1, CSAMT data inactivity bearing calibration of the invention the following steps are included:
Step 1: acquisition controlled-source audiomagnetotelluric discrete point apparent resistivity data, fits controllable source audio the earth
Electromagnetism original signal s (t) carries out empirical mode decomposition (Empirical Mode Decomposition, abbreviation EMD) to it and obtains
To n intrinsic mode function c1(t)、c2(t)、...ci(t)、...cn(t), and using each intrinsic mode function as IMF component；
As shown in Fig. 2, specific decomposition step is as follows:
(1), controlled-source audiomagnetotelluric original signal s (t) all maximum point and minimum point are found out；
(2), it connects all maximum point and minimum point to form upper packet respectively using cubic spline interpolation
Winding thread S1, lower envelope line S2；Envelope S is subtracted with each discrete point apparent resistivity value ρ1, lower envelope line S2Mean value (the figure of corresponding points
Middle m (t) is mean value matched curve) it obtains multiple discrete difference points and it is fitted to obtain first difference matched curve h11
(t)；Using cubic spline interpolation respectively by first difference matched curve h11(t) all maximum points and minimum point connect
It picks up to form difference signal coenvelope line H11, lower envelope line H12；Investigate first difference matched curve h11(t) whether meet
Two conditions of IMF, that is: 1. extreme point numbers (including maximum and minimum) and zero crossing number it is equal or at most difference 1
It is a；2. the corresponding difference signal coenvelope line H in arbitrary point11On point and difference signal lower envelope line H12On point average value be
0；By first difference matched curve h if meeting11(t) it is used as the 1st intrinsic mode function c1(t), step 2 is gone to；It is no
Then with first difference matched curve h11(t) substitution original signal s (t) repeats above-mentioned screening process, until from original signal s
(t) the difference matched curve for two conditions for meeting IMF is filtered out in and as the 1st intrinsic mode function c1(t)；c1
(t) radio-frequency component in original signal s (t) is represented；
(3) by remaining signal r1(t)=s (t)-c1(t) as new signal, continue to sieve according to step (1), (two)
Choosing, obtains second intrinsic mode function c2(t)；And so on, finally obtain n intrinsic mode function c1(t)、c2
(t)、...ci(t)、...cn(t) and a nondecomposable residual components rn(t): residual components rn(t) mean value of representation signal
Or trend term；Original signal s (t) so has just been broken down into n IMF component (i.e. n intrinsic mode function c1(t)、c2
(t)、...ci(t)、...cnAnd a trend term (i.e. residual components r (t))nThe sum of (t)):
The n intrinsic mode function c that empirical mode decomposition obtains1(t)、c2(t)、...ci(t)、...cn(t) according to frequency
Sequence arrangement from high to low.
Step 2: according to the cross-correlation coefficient r between the apparent resistivity value of whole frequency points of observation point and reference pointxyIdentification
Static relays, the specific method is as follows:
Firstly, if controlled source audio-frequency magnetotelluric data there are static relays, can be there are quiet in isogram
Occur significantly rise steeply band at the observation point of state effect (as shown in Fig. 6 (a) and Fig. 8 (a)).In order to confirm to rise steeply
Be implicitly present in static relays at the observation point of band, it would be desirable to calculate the observation point and around it observation point cross-correlation coefficient.
The characteristics of according to static relays, in log-log coordinate, the observation point curve influenced by static state and the reference not influenced by static state
Point curve form is constant, and in combination the characteristics of lower electrical consecutive variations, the apparent resistivity value of whole frequency points of observation point is regarded as
It is one group of data series xi, apparent resistivity value corresponding to whole frequency points with reference also to point is also one group of data series yi, by this
Two groups of data carry out relevant matches, seek cross-correlation coefficient r between the twoxy, it is believed that if tracing pattern is same or similar,
And their cross-correlation coefficient is big, illustrates that this is the data-bias as caused by static relays, should be corrected, whereas if
Related coefficient is small, then determines that this is by the truthful data of abnormal caused reflection underground electrical property.It is general choose close on observation point and
There is the apparent resistivity data of obvious numerical value difference, or the regional background apparent resistivity Value Data obtained by other means is made
For reference point apparent resistivity data.Cross-correlation coefficient rxyCalculation formula is as follows:
Wherein N is frequency point quantity,For the average value of whole frequency point apparent resistivity values of observation point,The whole of reference point
The average value of frequency point apparent resistivity value；
If rxy>=0.85 is thought there are static relays, and original signal s (t) is reconstructed according to step 3；
Step 3: original signal s (t) is reconstructed according to each IMF component that step 1 obtains, after obtaining static shift correction
Controlled-source audiomagnetotelluric reconstruction signal s'(t)；
The method that original signal s (t) is reconstructed is as follows:
(1) when reconstructing threshold value selection, if reconstruct threshold value is too small, possible static relays cannot be suppressed completely；If reconstruct
Threshold value is too big, can also mute the macrotectonics exception information of original signal s (t), it is contemplated that the characteristics of EMD decomposed signal, weight
Structure threshold value uses the adaptive threshold based on white Gaussian noise.If i-th of IMF component ci(t) reconstruct threshold value is Ti:
Wherein, c indicates threshold coefficient, and value range is 0.5-0.7 (analyzing by many experiments, preferably c=0.7)；M table
Show the signal length of selection, i.e., the discrete point number that each IMF component is chosen, EiIndicate i-th of IMF component ci(t) energy is strong
Degree, its calculation formula is:
Wherein,It is first IMF component c1(t) variance；
(2) h is used(ik)(t) i-th of IMF component c is indicatedi(t) in time tkThe discrete value at place, k=1,2 ... M, then i-th
A IMF component ci(t) corresponding reconstruct sub-signal is
Then controlled-source audiomagnetotelluric reconstruction signal is s'(t):
Step 4: HHT time-frequency spectrum is also known as HHT energy spectrum (abbreviation HHT spectrum, as shown in Figure 3), can show that controllable source sound
Frequency Magnetotelluric signal energy with time-frequency the specific regularity of distribution.It is equivalent to Fourier spectrum, but than Fourier spectrum with higher
Frequency resolution.Hilbert marginal spectrum (as shown in Figure 4) is by obtaining to Hilbert spectral integral, is a kind of signal spectrum
Estimation method provides amplitude (energy) distribution total to frequency, represents amplitude accumulation of the entire signal on time span
Effect breaches Fourier transformation only to the effective shortcoming of stationary signal.
Each IMF component c that step 1 is obtained1(t)、c2(t)、...ci(t)、...cn(t) it converts to obtain as Hilbert
Its instantaneous frequency and Hilbert spectrum and marginal spectrum；To IMF component ci(t) mathematic(al) representation for making Hilbert transformation is as follows:
In formula, yi(τ) is IMF component ci(t) signal after Hilbert transform, τ are signals after Hilbert transform
Time variable, yiThe analytic signal of (τ) is as follows:
zi(τ)=ci(t)+jyi(τ)
J is the imaginary part unit of analytic signal.
By Eulerian equation: ejθ=cos θ+jsin θ and yiKnown to the analytic signal expression formula of (τ):
It enablesThen it follows that
zi(τ)=ci(t)+jyi(τ)=a (τ) ejθ(τ),
Wherein ai(τ) is the amplitude of time-varying, and
Similarly, it can derive that the expression formula of phase theta (τ) is as follows:
Equally, instantaneous frequency ωi(τ) may be defined as:
When Variable Amplitude aiThe time-frequency distributions of (τ) are defined as:
H(ωi, τ) and=H (ωi(τ), τ) ≡ ai(τ)
H (ω, τ) is a wherein spectral line for Hilbert spectrum.
Finally, important Hilbert spectrum accumulation to together to get the Hilbert spectrum of original signal:
Hilbert marginal spectrum h (ω) can be equally expressed as followsin:
Wherein, h (ω) is the superposition of amplitude distribution in unit frequency.
Specific embodiment:
Fig. 5 be establish homogeneous half space model, phantom thicknesses be 900 meters, amount to 14 frequency points, respectively 1 to
8192Hz, survey line distance AB are 2000 meters, amount to 41 measuring points, and measuring point spacing is 50 meters.In the homogeneous half space, setting
One electric conductor can namely cause the static body of static relays.For the static posture near measuring point 20, top buried depth is 0.5
Rice, size 50m*50m, resistivity size are 500 Ω m.An exception is also provided with immediately below the static body
Body, the top buried depth of the anomalous body are 600 meters, size 200m*300m, and resistivity value is 1000 Ω m, the model
Background resistivity is 200 Ω m.
Fig. 6 (a) is the apparent resistivity isoline figure that model shown in Fig. 5 is handled without static shift correction, as shown in the figure,
It, can be true by the cross-correlation coefficient for calculating at these measuring points in measuring point 20 to there is the isopleth of shape of rising steeply between measuring point 22
Determine the presence of static relays.Fig. 6 (b) is the apparent resistivity isoline figure that static shift correction is carried out based on wavelet method, can from figure
To find out, static relays have obtained good inhibition, however, distortion phenomenon but has occurred in the position of anomalous body, with institute in model
The abnormal posture of setting, which is set, is deviated.Fig. 6 (c) is the apparent resistivity isoline figure that static shift correction is carried out based on the present invention, from
In figure, it is apparent that while static relays have obtained effective inhibition, the information of anomalous body is also accurately embodied
Out.
Fig. 7 is the homogeneous half space model that a surface is covered with cap rock, and model depth of cover is 50 meters, uniform below cap rock
Half space is with a thickness of 900 meters, total 14 frequency points, and respectively 1 to 8192Hz, survey line distance is 2000 meters, amount to 41 measuring points,
Measuring point spacing is 50 meters.It is provided with two static bodies in the model, is located near measuring point 10 and measuring point 16.Static body
Top buried depth is 50 meters, and size is 50m*50m, and resistivity size distinguishes 10 Ω m and 500 Ω m.In measuring point 17
To an anomalous body is also provided between measuring point 25, the top buried depth of the anomalous body is 650 meters, size 200m*
300m, resistivity value are 100 Ω m, and the background resistivity of the model is 50 Ω m.
Fig. 8 (a) is the apparent resistivity isoline figure that model shown in Fig. 7 is handled without static shift correction, as shown in the figure,
There is the isopleth for the shape that rises steeply near measuring point 10, can determine static relays by calculating the cross-correlation coefficient at the measuring point
In the presence of.And near measuring point 16, cap rock and the high resistance body of lower section have been connected in together, can not carry out effective differentiation of information, seriously
The data for hindering next step are explained.Fig. 8 (b) is the apparent resistivity isoline figure that static shift correction is carried out based on wavelet method, from
As can be seen that static relays have obtained good inhibition, however, the information of subsurface anomaly body is similarly suppressed in figure.
Fig. 8 (c) is the apparent resistivity isoline figure that static shift correction is carried out based on the present invention, from the graph, it is apparent that imitating in static state
While should having obtained effective inhibition, the information of anomalous body also accurately embodies out.
Claims (5)
1. a kind of CSAMT data inactivity bearing calibration, it is characterised in that the following steps are included:
Step 1: acquisition controlled-source audiomagnetotelluric discrete point apparent resistivity data, fits controlled-source audiomagnetotelluric
Original signal s (t) carries out empirical mode decomposition (Empirical Mode Decomposition, abbreviation EMD) to it and obtains n
A intrinsic mode function c1(t)、c2(t)、...ci(t)、...cn(t), and using each intrinsic mode function as IMF component；
Step 2: according to the cross-correlation coefficient r between the apparent resistivity value of whole frequency points of observation point and reference pointxyIdentification is static
Effect is then reconstructed original signal s (t) according to step 3 if there is static relays；
Step 3: original signal s (t) is reconstructed according to each IMF component that step 1 obtains, after obtaining static shift correction can
Control source audio magnetotelluric reconstruction signal s'(t)；
The method that original signal s (t) is reconstructed is as follows:
(1) using the adaptive threshold based on white Gaussian noise as reconstruct threshold value, if i-th of IMF component ci(t) reconstruct threshold
Value is Ti:
Wherein, c indicates threshold coefficient, value range 0.5-0.7；M indicate choose signal length, i.e., each IMF component from
Scatterplot number, EiIndicate i-th of IMF component ci(t) energy intensity, its calculation formula is:
Wherein,It is first IMF component c1(t) variance；
(2) h is used(ik)(t) i-th of IMF component c is indicatedi(t) in time tkThe discrete value at place, k=1,2 ... M, then i-th of IMF
Component ci(t) corresponding reconstruct sub-signal is
Then controlled-source audiomagnetotelluric reconstruction signal is s'(t):
2. CSAMT data inactivity bearing calibration according to claim 1, it is characterised in that the threshold coefficient c=
0.7。
3. CSAMT data inactivity bearing calibration according to claim 1, it is characterised in that right in the step one
The method that the apparent resistivity data progress empirical mode decomposition of CSAMT obtains multiple intrinsic mode function components is as follows:
(1), controlled-source audiomagnetotelluric original signal s (t) all maximum point and minimum point are found out；
(2), it connects all maximum point and minimum point to form coenvelope line respectively using cubic spline interpolation
S1, lower envelope line S2: coenvelope line S is subtracted with each discrete point apparent resistivity value ρ of original signal s (t)1With lower envelope line S2Pair
Mean value should be put, multiple discrete difference points are obtained and it is fitted to obtain first difference matched curve h11(t)；By first
Difference matched curve h11(t) all maximum points and minimum point connect to form difference signal coenvelope line H11, lower packet
Winding thread H12；Investigate first difference matched curve h11(t) whether meet two conditions of IMF, that is: extreme point number (including pole
Big value and minimum) and zero crossing number it is equal or most differ 1；The corresponding difference signal coenvelope line H in arbitrary point11On
Discrete point value and difference signal lower envelope line H12On discrete point value average value be 0；First difference is intended if meeting
Close curve h11(t) it is used as the 1st intrinsic mode function c1(t), step 3 is gone to；Otherwise with first difference matched curve h11
(t) substitution original signal s (t) repeats above-mentioned screening process, until filtering out two items for meeting IMF from original signal s (t)
The difference matched curve of part and as the 1st intrinsic mode function c1(t)；
(3) by remaining signal r1(t)=s (t)-c1(t) second is obtained according to step (1), (two) as new signal
Intrinsic mode function c2(t)；And so on, finally obtain n intrinsic mode function c1(t)、c2(t)、...ci(t)、...cn(t)
With a nondecomposable residual components rn(t):
N intrinsic mode function c1(t)、c2(t)、...ci(t)、...cn(t) it is used as IMF component.
4. CSAMT data inactivity bearing calibration according to claim 1, it is characterised in that in the step 2, identify quiet
State effect method is as follows: regarding the apparent resistivity value of whole frequency points of observation point as one group of data series xi, with reference also to point
Whole frequency points corresponding to apparent resistivity value be also one group of data series yi, this two groups of data are subjected to relevant matches, are sought
Cross-correlation coefficient r between the twoxy,
Wherein N is frequency point quantity,For the average value of whole frequency point apparent resistivity values of observation point,Whole frequency points of reference point
The average value of apparent resistivity value；
If rxy>=0.85 is thought that there are static relays.
5. CSAMT data inactivity bearing calibration according to claim 1, it is characterised in that further include following step:
Step 4: each IMF component c obtained to step 11(t)、c2(t)、...ci(t)、...cn(t) it is converted as Hilbert
To its instantaneous frequency and Hilbert spectrum and marginal spectrum；For any IMF component ci(t), instantaneous frequency ωi(τ), it is original
Hilbert spectrum H (ω, τ) of signal s (t) and the marginal spectrum h (ω) of original signal s (t) are as follows:
In formula, yi(τ) is IMF component ci(t) signal after Hilbert transform, when τ is signal after Hilbert transform
Between variable, θi(τ) is phase,
H(ωi, τ) and=H (ωi(τ),τ)≡ai(τ)
Wherein H (ωi, τ) and it is IMF component ci(t) spectral line in corresponding Hilbert spectrum, aiVariable Amplitude when (τ) is；
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811337149.3A CN109470135A (en) | 2018-11-12 | 2018-11-12 | CSAMT data inactivity bearing calibration |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811337149.3A CN109470135A (en) | 2018-11-12 | 2018-11-12 | CSAMT data inactivity bearing calibration |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109470135A true CN109470135A (en) | 2019-03-15 |
Family
ID=65671653
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811337149.3A Pending CN109470135A (en) | 2018-11-12 | 2018-11-12 | CSAMT data inactivity bearing calibration |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109470135A (en) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110348568A (en) * | 2019-07-16 | 2019-10-18 | 山东科技大学 | A kind of deep Mined-Out Areas method suitable for strong electromagnetic area |
| CN111352163A (en) * | 2020-03-03 | 2020-06-30 | 吉林大学 | Magnetotelluric depth sounding-based static effect correction method and system |
| CN111965712A (en) * | 2020-10-21 | 2020-11-20 | 国网江西省电力有限公司电力科学研究院 | Method for correcting static effect of controllable source audio magnetotelluric method |
| CN112099100A (en) * | 2020-08-25 | 2020-12-18 | 吉林大学 | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering |
| CN112379449A (en) * | 2020-10-30 | 2021-02-19 | 中国石油天然气集团有限公司 | Method and device for processing electromagnetic data of controllable source |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH03105279A (en) * | 1989-09-19 | 1991-05-02 | Central Res Inst Of Electric Power Ind | Receiving device for csamt method |
| CN105911595A (en) * | 2016-02-02 | 2016-08-31 | 中国科学院地质与地球物理研究所 | Method and apparatus for obtaining controllable source audio-frequency magnetotelluric (CSAMT) apparent phase information |
| CN106646666A (en) * | 2017-01-16 | 2017-05-10 | 中南大学 | Static effect correcting method based on plane wave electromagnetic sounding |
| CN107329183A (en) * | 2017-07-14 | 2017-11-07 | 中国地质科学院地球物理地球化学勘查研究所 | A kind of controlled-source audiomagnetotellurics sounding collecting method and device |
| CN108169800A (en) * | 2017-12-27 | 2018-06-15 | 江苏省有色金属华东地质勘查局地球化学勘查与海洋地质调查研究院 | Controlled-source audiomagnetotellurics method apparent resistivity near-field calibrating method |
-
2018
- 2018-11-12 CN CN201811337149.3A patent/CN109470135A/en active Pending
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH03105279A (en) * | 1989-09-19 | 1991-05-02 | Central Res Inst Of Electric Power Ind | Receiving device for csamt method |
| CN105911595A (en) * | 2016-02-02 | 2016-08-31 | 中国科学院地质与地球物理研究所 | Method and apparatus for obtaining controllable source audio-frequency magnetotelluric (CSAMT) apparent phase information |
| CN106646666A (en) * | 2017-01-16 | 2017-05-10 | 中南大学 | Static effect correcting method based on plane wave electromagnetic sounding |
| CN107329183A (en) * | 2017-07-14 | 2017-11-07 | 中国地质科学院地球物理地球化学勘查研究所 | A kind of controlled-source audiomagnetotellurics sounding collecting method and device |
| CN108169800A (en) * | 2017-12-27 | 2018-06-15 | 江苏省有色金属华东地质勘查局地球化学勘查与海洋地质调查研究院 | Controlled-source audiomagnetotellurics method apparent resistivity near-field calibrating method |
Non-Patent Citations (5)
| Title |
|---|
| BEN K.STERNBERG等: "Correction for the static shift in magnetotelluticusing transient electromagnetic soundings", 《GEOPHYSICS》 * |
| 于生宝等: "基于小波变换模极大值法和阈值法的CSAMT静态校正", 《地球物理学报》 * |
| 伍亮等: "大地电磁测深法中静态效应及其反演", 《地球物理学进展》 * |
| 罗延钟等: "可控源音频大地电磁法的静态效应校正", 《物探与化探》 * |
| 郑建波等: "C#与Matlab混合编程的CSAMT静态校正软件设计", 《实验室研究与探索》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110348568A (en) * | 2019-07-16 | 2019-10-18 | 山东科技大学 | A kind of deep Mined-Out Areas method suitable for strong electromagnetic area |
| CN111352163A (en) * | 2020-03-03 | 2020-06-30 | 吉林大学 | Magnetotelluric depth sounding-based static effect correction method and system |
| CN111352163B (en) * | 2020-03-03 | 2021-04-02 | 吉林大学 | Magnetotelluric depth sounding-based static effect correction method and system |
| CN112099100A (en) * | 2020-08-25 | 2020-12-18 | 吉林大学 | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering |
| CN111965712A (en) * | 2020-10-21 | 2020-11-20 | 国网江西省电力有限公司电力科学研究院 | Method for correcting static effect of controllable source audio magnetotelluric method |
| CN111965712B (en) * | 2020-10-21 | 2021-03-02 | 国网江西省电力有限公司电力科学研究院 | Method for correcting static effect of controllable source audio magnetotelluric method |
| CN112379449A (en) * | 2020-10-30 | 2021-02-19 | 中国石油天然气集团有限公司 | Method and device for processing electromagnetic data of controllable source |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109470135A (en) | CSAMT data inactivity bearing calibration | |
| Wang et al. | Time-frequency analysis of seismic data using synchrosqueezing transform | |
| US8311745B2 (en) | Method for spatial filtering of electromagnetic survey data | |
| Costabel et al. | Despiking of magnetic resonance signals in time and wavelet domains | |
| CN109100814B (en) | A kind of audio magnetotelluric method signal antinoise method based on noise classification | |
| CN106814402B (en) | Transient electromagnetic signal Prestack Noise Suppression Methods | |
| Wei et al. | Robust estimation of the discrete spectrum of relaxations for electromagnetic induction responses | |
| Liu et al. | A modified empirical mode decomposition method for multiperiod time-series detrending and the application in full-waveform induced polarization data | |
| CN109188542B (en) | Far reference magnetotelluric impedance calculation method for wave zone correlation detection | |
| CN107526103B (en) | The acquiring method of Processing Seismic Data and its threshold and useful signal frequency | |
| CN108680782A (en) | The voltage flicker parameter detection method decomposed based on extreme point symmetric pattern | |
| CN110135022B (en) | Generalized skin depth calculation method based on polarized medium model | |
| LIU et al. | Robust estimation method of sea magnetotelluric impedance based on correlative coefficient | |
| Volkovitsky et al. | Airborne EM systems variety: what is the difference? | |
| CN111665569B (en) | Dual-mode frequency domain passive source electric field method | |
| Nelson et al. | Theoretical response of a time-domain, airborne, electromagnetic system | |
| Xu et al. | Combining the Wiener Filter With Calibration Device to Improve the Accuracy of the Helicopter Transient Electromagnetic System | |
| CN114924328B (en) | Urban artificial source electromagnetic exploration method and system with vertical magnetic field reference channel | |
| Song et al. | Emi-based diagnosis to grounding grids by combining ensemble empirical mode decomposition and ICA | |
| Zach et al. | Data preprocessing and starting model preparation for 3D inversion of marine CSEM surveys | |
| Lavrukhin et al. | Search for the optimal method for obtaining estimations for an automated magnetotelluric signal processing system | |
| Bo et al. | New correction method for controlled-source electromagnetics source effects | |
| CN111308561B (en) | Method for removing strong noise of electromagnetic signal | |
| Cuevas et al. | 1D layered resolution analysis of two marine EM exploration methods | |
| Rupf et al. | New approaches for automated data processing of annually laminated sediments |
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 | ||
| AD01 | Patent right deemed abandoned | ||
| AD01 | Patent right deemed abandoned |
Effective date of abandoning: 20200911 |