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
- signal
- points
- imf component
- data
- difference
- 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
- 230000003068 static effect Effects 0.000 claims abstract description 74
- 238000000034 method Methods 0.000 claims abstract description 49
- 238000012937 correction Methods 0.000 claims abstract description 33
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 10
- 238000001228 spectrum Methods 0.000 claims description 32
- 230000008569 process Effects 0.000 claims description 7
- 230000009466 transformation Effects 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 5
- 230000000007 visual effect Effects 0.000 claims description 5
- 230000002547 anomalous effect Effects 0.000 abstract description 3
- 230000002159 abnormal effect Effects 0.000 description 9
- 230000000694 effects Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000005856 abnormality Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000013519 translation Methods 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000001149 cognitive effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
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 invention relates to a static correction method aiming at CSAMT data.
Background
Controllable Source Audio frequency magnetotelluric method (CSAMT Carny resistivity sounding curve for short), therefore, the method is also called controllable source audio magnetotelluric sounding, and has the advantages of high working efficiency, large exploration depth range, high spatial resolution, small topographic influence, small high resistance shielding effect and the like, in the signal acquisition process of the CSAMT, local conductive inhomogeneities often exist on the earth surface or near the earth surface, the apparent resistivity curves of adjacent measuring points in an apparent resistivity-frequency log coordinate system move up and down along an apparent resistivity axis and are reflected on an apparent resistivity pseudo-section diagram, a steep dense zone with a small transverse range can appear, the underground abnormity caused by an ore body or a target object is covered, and misleading is caused to the judgment of the actual geological condition.
The correction method of the static effect can be divided into (1) static correction is carried out by a curve translation method; (2) performing static correction by a numerical analysis method; (3) static correction is carried out by a joint inversion method; (4) performing static correction by using the phase conversion data; (5) performing static correction by using magnetic field data; (6) direct two-dimensional and multi-dimensional inversion methods; (7) wavelet filtering methods, etc. The effect of the curve translation method depends on the cognitive level of a processor on data to a great extent, the judgment and the processing experience of the static effect play a very key role, and the abnormality which is not caused by the static effect is easily mistakenly taken as the static effect to be corrected; for some static effects with simple terrain background, especially static effects generated by terrain influence, a numerical simulation method can be adopted for identification and correction, and under actual field conditions, the static effects are often unknown and can only be used as reference research; the joint inversion method is equivalent to two geophysical methods which are repeated, so that the cost is high and the practicability is not strong; the wavelet transform is a variable window transform, which visually distinguishes static effect caused by surface unevenness from large structure abnormity in mathematics, and has higher resolution and correction precision.
The Hilbert-Huang Transform (HHT for short) method is a new signal processing method proposed by N.E. Huang and the like of the American national aerospace agency, has clear physical significance for non-stationary and non-linear signals, can obtain the characteristics of amplitude-time-frequency distribution of the signals, and has self-adaptability. Aiming at the problems of the traditional method in suppressing the static effect, the invention provides static correction of CSAMT data based on HHT.
Disclosure of Invention
The invention aims to provide a CSAMT data static correction method which has no problem of basis function selection, can suppress static effect and simultaneously reserve large structure abnormal information while realizing static correction.
In order to solve the above technical problem, the CSAMT data static correction method of the present invention comprises the following steps:
step one, collecting controllable source audio magnetotelluric discrete point-view resistivity data, fitting a controllable source audio magnetotelluric original signal s (t), and performing Empirical Mode Decomposition (EMP) on the signal to obtain n eigenmode functions c (n) of the signal1(t)、c2(t)、...ci(t)、...cn(t) and taking each eigenmode function as an IMF component;
step two, according to the cross correlation coefficient r between the visual resistivity values of all the frequency points of the observation point and the reference pointxyIdentifying static effect, if the static effect exists, reconstructing the original signal s (t) according to the step three;
thirdly, reconstructing the original signal s (t) according to each IMF component obtained in the first step to obtain a statically corrected controllable source audio magnetotelluric reconstruction signal s' (t);
the method for reconstructing the original signal s (t) is as follows:
(1) adopting a self-adaptive threshold value based on Gaussian white noise as a reconstruction threshold value, and setting an ith IMF component ci(T) a reconstruction threshold of Ti：
Wherein c represents a threshold coefficient, and the value range is 0.5-0.7; m represents the selected signal length, i.e. the number of discrete points for each IMF component, EiRepresenting the ith IMF component ci(t) an energy intensity calculated by the formula:
wherein,is the first IMF component c1(t) variance;
(2) by using h(ik)(t) denotes the i-th IMF component ci(t) at time tkM, the i-th IMF component ci(t) the corresponding reconstructed sub-signal is
The controllable source audio magnetotelluric reconstructed signal is s' (t):
the threshold coefficient is preferably 0.7.
In the first step, the method for obtaining a plurality of intrinsic mode function components by performing empirical mode decomposition on the apparent resistivity data of the CSAMT is as follows:
finding out all maximum value points and minimum value points of a controllable source audio magnetotelluric original signal s (t);
(II) respectively connecting all the maximum value points and minimum value points by adopting a cubic spline interpolation method to form an upper envelope line S1Lower envelope S2: subtracting the upper envelope S from the discrete point-view resistivity values ρ of the original signal S (t)1And the lower envelope S2Obtaining a plurality of discrete difference points and fitting the discrete difference points to obtain a first difference fitting curve h11(t); fitting the first difference to a curve h11(t) all the maximum and minimum points are connected to form the upper envelope H of the difference signal11Lower envelope H12(ii) a Examine the first difference fit curve h11(t) whether two conditions for IMF are satisfied, namely: the number of extreme points (including maximum and minimum) and the number of zero-crossing points are equal or differ by at most 1; difference signal upper envelope line H corresponding to any point11Upper discrete point value and lower envelope H of difference signal12The average value of the discrete point values above is 0; if yes, fitting the first difference value to a curve h11(t) as a function of the 1 st eigenmodec1(t), go to step three; otherwise, fitting the curve h with the first difference11(t) repeating the screening process instead of the original signal s (t) until a difference fitting curve satisfying two conditions of IMF is screened from the original signal s (t) and taken as the 1 st intrinsic mode function c1(t)；
(III) the rest signal r1(t)＝s(t)-c1(t) as a new signal, according to the steps (one) and (two), obtaining a second eigenmode function c2(t); repeating the steps to obtain n intrinsic mode functions c1(t)、c2(t)、...ci(t)、...cn(t) and an irresolvable residual component rn(t)：
n eigenmode functions c1(t)、c2(t)、...ci(t)、...cn(t) as IMF component.
In the second step, the static effect identification method comprises the following step of regarding the apparent resistivity values of all frequency points of the observation points as a group of data series xiThe apparent resistivity values corresponding to all frequency points of the same reference point are also a group of data series yiThe two groups of data are subjected to correlation matching, and the cross correlation coefficient r between the two groups of data is obtainedxy，
Wherein N is the number of frequency points,is the average value of all frequency point visual resistivity values of the observation points,reference toAverage value of all frequency point video resistivity values of points;
if rxyA static effect is considered to be present at > 0.85.
The invention may also include the steps of:
step four, each IMF component c obtained in the step one1(t)、c2(t)、...ci(t)、...cn(t) carrying out Hilbert transformation to obtain an instantaneous frequency, a Hilbert spectrum and a marginal spectrum of the Hilbert transformation; for any IMF component ci(t) instantaneous frequency of ωi(τ), the Hilbert spectrum H (ω, τ) of the original signal s (t) and the marginal spectrum H (ω) of the original signal s (t) are as follows:
in the formula, yi(τ) is the IMF component ci(t) Hilbert transformed signal,. tau.is the time variable of the Hilbert transformed signal,. theta.i(tau) is the phase of the phase,
H(ωi,τ)＝H(ωi(τ),τ)≡ai(τ)
wherein H (ω)iτ) is the IMF component ci(t) a corresponding one of the Hilbert spectra, ai(τ) is a time-varying amplitude;
the invention has the beneficial effects that: (1) the traditional static correction method (such as a wavelet method) has the problem of difficult basis function selection, but the method provided by the invention does not have the problem, and the limitation condition is avoided. (2) The CSAMT static correction provided by the invention can suppress the static effect and simultaneously keep the information of the abnormal body in the original data so as not to distort. And (3) comparing the Hilbert spectrum and the marginal spectrum h (omega) obtained in the step four with the HHT spectrum and the HHT marginal spectrum of the original data, and intuitively reflecting the correction condition of the static effect. Different frequencies of each time interval of the signal can be clearly shown in the time frequency spectrum, and the dominant frequency range existing in the signal can be seen from the marginal spectrum.
Description of the drawings:
FIG. 1: the invention discloses a flow chart of a CSAMT data static correction method.
FIG. 2: and (3) a schematic process diagram of screening IMF components by EMD decomposition of controllable source audio magnetotelluric data.
FIG. 3: HHT spectra of controlled source audio magnetotelluric data.
FIG. 4: HHT margin spectra of the controllable source audio magnetotelluric data.
FIG. 5: a uniform half-space forward model.
FIG. 6(a), FIG. 6(b), and FIG. 6(c) are respectively apparent resistivity contour maps of the uniform half-space forward model without static correction processing, static correction processing based on wavelet method, and static correction using the present invention.
FIG. 7: a uniform half-space forward model containing a cap layer.
FIG. 8(a), FIG. 8(b), and FIG. 8(c) are apparent resistivity contour maps of the uniform half-space forward model with cap layer without static correction, with wavelet-based static correction, and with the present invention.
The specific implementation mode is as follows:
the present invention will be described in further detail with reference to the following drawings and examples.
As shown in fig. 1, the CSAMT data static correction method of the present invention comprises the following steps:
step one, collecting controllable source audio magnetotelluric discrete point-view resistivity data, fitting a controllable source audio magnetotelluric original signal s (t), and performing Empirical Mode Decomposition (EMP) on the signal to obtain n eigenmode functions c (n) of the signal1(t)、c2(t)、...ci(t)、...cn(t) and taking each eigenmode function as an IMF component; as shown in fig. 2, the specific decomposition steps are as follows:
finding out all maximum value points and minimum value points of a controllable source audio magnetotelluric original signal s (t);
(II) respectively connecting all the maximum value points and minimum value points by adopting a cubic spline interpolation method to form an upper envelope line S1Lower envelope S2(ii) a Subtracting the envelope S from each discrete point view resistivity value ρ1Lower envelope S2Obtaining a plurality of discrete difference points by the mean value of the corresponding points (m (t) is a mean value fitting curve in the figure) and fitting the discrete difference points to obtain a first difference value fitting curve h11(t); fitting the first difference value to a curve h by a cubic spline interpolation method11(t) all the maximum and minimum points are connected to form the upper envelope H of the difference signal11Lower envelope H12(ii) a Examine the first difference fit curve h11(t) two bars of whether or not IMF is satisfiedThe piece, namely: 1. the number of extreme points (including maximum and minimum) and the number of zero-crossing points are equal or differ by at most 1; 2. difference signal upper envelope line H corresponding to any point11Upper point and lower envelope H of difference signal12The average value of the points above is 0; if yes, fitting the first difference value to a curve h11(t) as the 1 st eigenmode function c1(t), go to step two; otherwise, fitting the curve h with the first difference11(t) repeating the screening process instead of the original signal s (t) until a difference fitting curve satisfying two conditions of IMF is screened from the original signal s (t) and taken as the 1 st intrinsic mode function c1(t)；c1(t) represents the high frequency components in the original signal s (t);
(III) the rest signal r1(t)＝s(t)-c1(t) as a new signal, and continuously screening according to the steps (I) and (II) to obtain a second eigenmode function c2(t); repeating the steps to obtain n intrinsic mode functions c1(t)、c2(t)、...ci(t)、...cn(t) and an irresolvable residual component rn(t): residual component rn(t) represents a mean or trend term of the signal; the original signal s (t) is thus decomposed into n IMF components (i.e., n eigenmode functions c)1(t)、c2(t)、...ci(t)、...cn(t)) and a trend term (i.e., residual component rn(t)) sum:
n intrinsic mode functions c obtained by empirical mode decomposition1(t)、c2(t)、...ci(t)、...cn(t) are arranged in the order of frequency from high to low.
Step two, according to the cross correlation coefficient r between the visual resistivity values of all the frequency points of the observation point and the reference pointxyThe static effect is identified by the following specific method:
first, if there is a static effect in the controllable source audio magnetotelluric data, a sharp steep band may appear in the contour map thereof at the observation point where the static effect is present (as shown in fig. 6(a) and 8 (a)). In order to be able to verify that there is indeed a static effect at an observation point where a steep band appears, we need to calculate the cross-correlation coefficient of that observation point with its surrounding observation points. According to the characteristics of static effect, in a log-log coordinate, the curve form of an observation point affected by static state and the curve form of a reference point not affected by static state are unchanged, and the apparent resistivity values of all frequency points of the observation point are regarded as a group of data series x by combining the characteristic of continuous change of underground electrical propertyiThe apparent resistivity values corresponding to all frequency points of the same reference point are also a group of data series yiThe two groups of data are subjected to correlation matching, and the cross correlation coefficient r between the two groups of data is obtainedxyIf the curves are identical or similar in shape and the cross-correlation coefficient is large, it is considered that the data deviation caused by the static effect is the data deviation, and the data deviation should be corrected, otherwise, if the correlation coefficient is small, the data is judged to be the real data which is caused by the abnormality and reflects the subsurface electrical property. Generally, apparent resistivity data which is close to an observation point and has obvious numerical difference or regional background apparent resistivity value data obtained by other means is selected as reference point apparent resistivity data. Cross correlation coefficient rxyThe calculation formula is as follows:
wherein N is the number of frequency points,is the average value of all frequency point visual resistivity values of the observation points,the average value of all frequency point video resistivity values of the reference point;
if rxyThe static effect is considered to exist when the number is more than or equal to 0.85, and the three pairs of steps are adoptedReconstructing the original signal s (t);
thirdly, reconstructing the original signal s (t) according to each IMF component obtained in the first step to obtain a statically corrected controllable source audio magnetotelluric reconstruction signal s' (t);
the method for reconstructing the original signal s (t) is as follows:
(1) in the reconstruction threshold selection, if the reconstruction threshold is too small, the static effect may not be completely suppressed; if the reconstruction threshold is too large, the large structure abnormal information of the original signal s (t) can be suppressed, and the reconstruction threshold adopts a self-adaptive threshold based on Gaussian white noise in consideration of the characteristics of EMD decomposition signals. Let i' th IMF component ci(T) a reconstruction threshold of Ti：
Wherein c represents a threshold coefficient, and the value range is 0.5-0.7 (after a large number of experimental analyses, c is preferably 0.7); m represents the selected signal length, i.e. the number of discrete points selected for each IMF component, EiRepresenting the ith IMF component ci(t) an energy intensity calculated by the formula:
wherein,is the first IMF component c1(t) variance;
(2) by using h(ik)(t) denotes the i-th IMF component ci(t) at time tkM, the i-th IMF component ci(t) the corresponding reconstructed sub-signal is
The controllable source audio magnetotelluric reconstructed signal is s' (t):
step four: the HHT time frequency spectrum is also called as an HHT energy spectrum (HHT spectrum for short, as shown in figure 3), and can display the specific distribution rule of the controllable source audio magnetotelluric signal energy along with time frequency. It is equivalent to the fourier spectrum but has a higher frequency resolution than the fourier spectrum. The Hilbert marginal spectrum (shown in FIG. 4) is obtained by integrating the Hilbert spectrum, is a signal spectrum estimation method, provides total amplitude (energy) distribution of frequency, represents the amplitude accumulation effect of the whole signal on a time span, and overcomes the defect that Fourier transform is only effective on a stable signal.
For each IMF component c obtained in the step one1(t)、c2(t)、...ci(t)、...cn(t) carrying out Hilbert transformation to obtain an instantaneous frequency, a Hilbert spectrum and a marginal spectrum of the Hilbert transformation; for IMF component ci(t) the mathematical expression for the Hilbert transform is as follows:
in the formula, yi(τ) is the IMF component ci(t) Hilbert transformed signal, τ being the time variable of the Hilbert transformed signal, yiThe analytical signals of (τ) are as follows:
zi(τ)＝ci(t)+jyi(τ)
j is the imaginary unit of the analytic signal.
From the euler equation: e.g. of the typejθCos θ + jsin θ and yiThe analytic signal expression of (τ) is known as:
order toThen it can be derived:
zi(τ)＝ci(t)+jyi(τ)＝a(τ)·ejθ(τ)，
wherein a isi(τ) is a time-varying amplitude, and
similarly, the expression for the phase θ (τ) can be derived as follows:
also, instantaneous frequency ωi(τ) can be defined as:
time-varying amplitude aiThe time-frequency distribution of (τ) is defined as:
H(ωi，τ)＝H(ωi(τ)，τ)≡ai(τ)
h (omega, tau) is one of the spectral lines of the Hilbert spectrum.
And finally, accumulating Hilbert spectrums of all components together to obtain the Hilbert spectrums of the original signals:
the Hilbert marginal spectrum h (ω) can also be expressed as follows:
where h (ω) is the superposition of the amplitude distribution within the unit frequency.
The specific embodiment is as follows:
FIG. 5 is a built uniform half-space model, the thickness of the model is 900 meters, the total number of 14 frequency points is 1-8192 Hz, the line distance AB is 2000 meters, the total number of 41 measuring points is measured, and the distance between the measuring points is 50 meters. In this uniform half-space, an electrically conductive body, i.e. a static body, is arranged, which will cause a static effect. The static body is located near the measuring point 20, the top buried depth is 0.5 m, the size is 50m x 50m, and the resistivity is 500 omega m. An abnormal body is further arranged right below the static body, the top buried depth of the abnormal body is 600 meters, the size is 200m x 300m, the resistivity value is 1000 omega m, and the background resistivity of the model is 200 omega m.
FIG. 6(a) is an apparent resistivity contour map of the model shown in FIG. 5 without the static correction process, as shown by the steep contour lines between measurement points 20 and 22, and the presence of the static effect can be determined by calculating the cross-correlation coefficients at these several measurement points. Fig. 6(b) is an apparent resistivity contour map of a static correction based on a wavelet method, and it can be seen from the figure that the static effect is well suppressed, however, the position of the abnormal body is distorted and deviated from the position of the abnormal body set in the model. Fig. 6(c) is an apparent resistivity contour map subjected to static correction based on the present invention, and it is apparent from the figure that information on an anomalous body is accurately reflected while the static effect is effectively suppressed.
FIG. 7 is a uniform half-space model with a cover layer covering the surface, the cover layer thickness of the model is 50m, the uniform half-space thickness under the cover layer is 900 m, 14 frequency points are counted, the frequency points are respectively 1 to 8192Hz, the measuring line distance is 2000 m, 41 measuring points are counted, and the measuring point distance is 50 m. Two static bodies are provided in the model, one near each measuring point 10 and one near each measuring point 16. The top buried depth of the static body is 50 meters, the size is 50m x 50m, and the resistivity is 10 omega m and 500 omega m respectively. An abnormal body is further arranged between the measuring points 17 and 25, the top buried depth of the abnormal body is 650 meters, the size is 200m x 300m, the resistivity value is 100 omega m, and the background resistivity of the model is 50 omega m.
FIG. 8(a) is an apparent resistivity contour plot of the model shown in FIG. 7 without a static correction process, as shown by the presence of a steep contour near the measurement point 10 at which the cross-correlation coefficient was calculated to determine the presence of a static effect. In the vicinity of the measuring point 16, the cover layer is connected with the lower high-resistance body, so that effective judgment of information cannot be performed, and further data interpretation is seriously hindered. FIG. 8(b) is an apparent resistivity contour map statically corrected based on the wavelet method, from which it can be seen that the static effect is well suppressed, however, the information of subsurface anomalies is also suppressed. Fig. 8(c) is an apparent resistivity contour map subjected to static correction based on the present invention, and it is apparent from the figure that information on an anomalous body is accurately reflected while the static effect is effectively suppressed.
Claims (5)
1. A CSAMT data static correction method is characterized by comprising the following steps:
step one, collecting controllable source audio magnetotelluric discrete point-view resistivity data, fitting a controllable source audio magnetotelluric original signal s (t), and performing Empirical Mode Decomposition (EMP) on the signal to obtain n eigenmode functions c (n) of the signal1(t)、c2(t)、...ci(t)、...cn(t) and taking each eigenmode function as an IMF component;
step two, according to observation point and reference pointCross correlation coefficient r between apparent resistivity values of all frequency pointsxyIdentifying static effect, if the static effect exists, reconstructing the original signal s (t) according to the step three;
thirdly, reconstructing the original signal s (t) according to each IMF component obtained in the first step to obtain a statically corrected controllable source audio magnetotelluric reconstruction signal s' (t);
the method for reconstructing the original signal s (t) is as follows:
(1) adopting a self-adaptive threshold value based on Gaussian white noise as a reconstruction threshold value, and setting an ith IMF component ci(T) a reconstruction threshold of Ti：
Wherein c represents a threshold coefficient, and the value range is 0.5-0.7; m represents the selected signal length, i.e. the number of discrete points for each IMF component, EiRepresenting the ith IMF component ci(t) an energy intensity calculated by the formula:
wherein,is the first IMF component c1(t) variance;
(2) by using h(ik)(t) denotes the i-th IMF component ci(t) at time tkM, the i-th IMF component ci(t) the corresponding reconstructed sub-signal is
The controllable source audio magnetotelluric reconstructed signal is s' (t):
2. the CSAMT data static correction method according to claim 1, wherein the threshold coefficient c is 0.7.
3. The method for statically correcting CSAMT data according to claim 1, wherein in the first step, the empirical mode decomposition of the apparent resistivity data of CSAMT to obtain a plurality of intrinsic mode function components is as follows:
finding out all maximum value points and minimum value points of a controllable source audio magnetotelluric original signal s (t);
(II) respectively connecting all the maximum value points and minimum value points by adopting a cubic spline interpolation method to form an upper envelope line S1Lower envelope S2: subtracting the upper envelope S from the discrete point-view resistivity values ρ of the original signal S (t)1And the lower envelope S2Obtaining a plurality of discrete difference points and fitting the discrete difference points to obtain a first difference fitting curve h11(t); fitting the first difference to a curve h11(t) all the maximum and minimum points are connected to form the upper envelope H of the difference signal11Lower envelope H12(ii) a Examine the first difference fit curve h11(t) whether two conditions for IMF are satisfied, namely: the number of extreme points (including maximum and minimum) and the number of zero-crossing points are equal or differ by at most 1; difference signal upper envelope line H corresponding to any point11Upper discrete point value and lower envelope H of difference signal12The average value of the discrete point values above is 0; if yes, fitting the first difference value to a curve h11(t) as the 1 st eigenmode function c1(t), go to step three; otherwise, fitting the curve h with the first difference11(t) repeating the above screening process in place of the original signal s (t) until a difference fit curve satisfying two conditions of IMF is screened from the original signal s (t)Line and take it as the 1 st eigenmode function c1(t)；
(III) the rest signal r1(t)＝s(t)-c1(t) as a new signal, according to the steps (one) and (two), obtaining a second eigenmode function c2(t); repeating the steps to obtain n intrinsic mode functions c1(t)、c2(t)、...ci(t)、...cn(t) and an irresolvable residual component rn(t)：
n eigenmode functions c1(t)、c2(t)、...ci(t)、...cn(t) as IMF component.
4. The CSAMT data static correction method of claim 1, wherein in the second step, the static effect is identified by considering all the frequency point apparent resistivity values of the observation points as a group of data series xiThe apparent resistivity values corresponding to all frequency points of the same reference point are also a group of data series yiThe two groups of data are subjected to correlation matching, and the cross correlation coefficient r between the two groups of data is obtainedxy，
Wherein N is the number of frequency points,is the average value of all frequency point visual resistivity values of the observation points,the average value of all frequency point video resistivity values of the reference point;
if rxyA static effect is considered to be present at > 0.85.
5. The method for statically correcting CSAMT data according to claim 1, further comprising the steps of:
step four, each IMF component c obtained in the step one1(t)、c2(t)、...ci(t)、...cn(t) carrying out Hilbert transformation to obtain an instantaneous frequency, a Hilbert spectrum and a marginal spectrum of the Hilbert transformation; for any IMF component ci(t) instantaneous frequency of ωi(τ), the Hilbert spectrum H (ω, τ) of the original signal s (t) and the marginal spectrum H (ω) of the original signal s (t) are as follows:
in the formula, yi(τ) is the IMF component ci(t) Hilbert transformed signal,. tau.is the time variable of the Hilbert transformed signal,. theta.i(tau) is the phase of the phase,
H(ωi，τ)＝H(ωi(τ),τ)≡ai(τ)
wherein H (ω)iτ) is the IMF component ci(t) a corresponding one of the Hilbert spectra, ai(τ) is a time-varying amplitude;
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 | |
| Reninger et al. | Singular value decomposition as a denoising tool for airborne time domain electromagnetic data | |
| Pinnegar et al. | The S-transform with windows of arbitrary and varying shape | |
| Huang et al. | Time dependent intrinsic correlation analysis of temperature and dissolved oxygen time series using empirical mode decomposition | |
| US8090539B2 (en) | Method for spatial filtering of electromagnetic survey data | |
| Costabel et al. | Despiking of magnetic resonance signals in time and wavelet domains | |
| Varentsov et al. | System of electromagnetic field transfer operators for the BEAR array of simultaneous soundings: Methods and results | |
| Li et al. | Dictionary learning and shift-invariant sparse coding denoising for controlled-source electromagnetic data combined with complementary ensemble empirical mode decomposition | |
| Ji et al. | Noise reduction of grounded electrical source airborne transient electromagnetic data using an exponential fitting-adaptive Kalman filter | |
| CN109581516B (en) | Denoising method and system for data of curvelet domain statistic adaptive threshold value ground penetrating radar | |
| CN108345039B (en) | A method of eliminating adjacent frequency harmonic wave interference in ground nuclear magnetic resonance data | |
| CN109765624A (en) | A kind of frequency domain aviation electromagnetic data de-noising method based on variation mode decomposition | |
| Kremer et al. | Review of acquisition and signal processing methods for electromagnetic noise reduction and retrieval of surface nuclear magnetic resonance parameters | |
| Liu et al. | Investigating DEM error patterns by directional variograms and Fourier analysis | |
| Liu et al. | A modified empirical mode decomposition method for multiperiod time-series detrending and the application in full-waveform induced polarization data | |
| CN109188542A (en) | A kind of the remote of wave area correlation detection refers to magnetotelluric impedance computation method | |
| Chen et al. | Two-dimensional nonlinear geophysical data filtering using the multidimensional EEMD method | |
| CN111289800B (en) | Small-resistance vibration monitoring method based on generalized regression neural network | |
| Chen et al. | Natural logarithm transformed EEMD instantaneous attributes of reflection data | |
| CN111142134A (en) | Coordinate time series processing method and device | |
| Fan et al. | An elaborately designed virtual frame to level aeromagnetic data | |
| CN104570131A (en) | Method and device for estimating magnetotelluric parameters | |
| CN110109184B (en) | Passive field source type three-dimensional electric field exploration method based on multiple daily change points | |
| LIU et al. | Robust estimation method of sea magnetotelluric impedance based on correlative coefficient | |
| CN114966874B (en) | Forward-modeling-based frequency domain electromagnetic method time sequence data synthesis method |
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 |