The noise-reduction method of transient electromagnetic detecting echo signal
(1) technical field
The present invention relates to Transient Electric Magnetic Field Testing Technology, be specially a kind of noise-reduction method of the transient electromagnetic detecting echo signal based on empirical mode decomposition.
(2) background technology
As everyone knows, the principle of transient electromagnetic detecting is: the transient electrical magnetic wave is subject to the decay of different medium in each stratum in the process of underground propagation, and different medium is different to the transient electromagnetic wave attenuation value of different frequency.Therefore, as long as extract in the echoed signal characteristic frequency spectrum corresponding to the transient electromagnetic signal of different material, just can be finally inversed by different material in the stratum.Yet, geologic remote sensing signal right and wrong stably, and ground unrest is complicated, adopts conventional time-frequency combination analytic approach, can produce spurious signal and false frequency, finally causes system's false alarm rate too high.
The signal that the transient electrical magnetic wave collects when surveying the stratum unavoidably is being mingled with multiple noise.The echoed signal that the needs that comprise in the reception signal obtain is quite faint, particularly for the signal of the deep layer that receives, wherein the neighbourhood noise around the frequency of echoed signal and the energy Ratios is all much smaller, and the noise of this moment is compared with echoed signal, is very noisy.In order to obtain the accurately echoed signal of formation information, must receive to this signal and carry out noise reduction process, therefrom extract the echoed signal that needs.At present mostly being to convert thereof into first linear stationary signal for the method for the filtering of nonlinear and nonstationary signal and noise reduction process processes again.So just there is a priori problem, namely must knows the feature of nonlinear properties, just can convert thereof into linear stationary signal, then carry out noise reduction process.And the judgement of the feature prior imformation of nonlinear properties is great difficult points, and noise reduction process effect thereafter is not good yet.
It is typical non-stationary nonlinear properties that transient electromagnetic detects the echoed signal that receives.The instrument of present various utilization transient electromagnetic methods, just by simple integration stack, the number of times that broadly increases stack reaches the purpose of noise reduction.Although can eliminate the noise of a part, can realize extracting for the larger echoed signal of amplitude, that is to say that the echoed signal noise reduction for shallow-layer is to produce effect.
In the echoed signal for deep layer, the large manyfold of noise ratio echo signal, how no matter the method for integration stack increase stacking fold, final signal is still take noise as main, can not eliminate the noise of the echoed signal of deep layer because of the method for integration stack, also just can not extract with this method the echo signal of deep layer geological condition, can't accurately obtain the deep layer geologic feature.
In order to make transient electromagnetic detecting be used for the noise-reduction method that deep layer geology must be studied new reception signal.
(3) summary of the invention
The objective of the invention is to disclose a kind of noise-reduction method of the transient electromagnetic detecting echo signal based on empirical mode decomposition, transient electromagnetic detecting is received signal carry out Hilbert-Huang transform (the English Hilbert-Huang Transform of being, it is abbreviated as HHT) empirical mode decomposition (EMD), extreme value with Cubic Spline Functions Fitting transient electromagnetic signal, the fluctuation at Inhibitory signal two ends, find the eigenfunction orthogonal set intrinsic mode function component of conversion, decomposite complete geology spectrum information.
Transient electromagnetic detecting of the present invention receives the signal de-noising method, and step is as follows:
I, transient electromagnetic detecting is received signal carry out empirical mode decomposition (EMD)
All extreme points of picked up signal data x (t) at first, all local maximums are formed the coenvelope of data with cubic spline functions, in like manner, all local minimums are formed the lower envelope of data with cubic spline functions, upper lower envelope should cover all data points, and its average is denoted as m
1, from former data sequence, deduct this average and obtain first component h
1:
h
1=x(t)-m
1
Said process is the empirical mode decomposition of Hilbert-Huang transform, is also referred to as " sieve ".
In the ideal case, h
1It is an intrinsic mode function (IMF) component.But the match overshoot (overshoot) of envelope and owe the punching (undershoot) be ubiquitous, " sieve " can produce new extreme point, displacement or amplify already present extreme point, perhaps so that in the waveform a slight variation zoom into a Local Extremum.Another problem of the process of " sieve " is for nonlinear data, there is error between envelope average and the real local mean value, no matter data are through how many times " sieve ", and some asymmetric waveform still exist, and this mainly is owing to adopted the average approximation method of envelope.
The process of " sieve " has two purposes: remove the stack ripple and make waveform more symmetrical.In order to reach this two purposes, the process of " sieve " must be carried out repeatedly.In the process of the second time " sieve ", h
1Regard pending data as, its envelope average is m
11, then have
h
11=h
1-m
11
This process repeats k time,
h
1k=h
1(k-1)-m
1k
Until h
1kSatisfy till the condition of intrinsic mode function, so just from the original signal data, decomposite first intrinsic mode function component, be denoted as:
c
1=h
1k
The process of " sieve " has been removed the stack ripple, this instantaneous frequency that the intrinsic mode function component that obtains is obtained behind Hilbert-Huang transform is meaningful, simultaneously, the process of " sieve " also so that the very large adjacent wave deformation of amplitude variations get smoothly, this might remove significant amplitude fluctuation, obtain amplitude constant, only have warbled intrinsic mode function component.In order to obtain all significant intrinsic mode function components of frequency modulation (PFM) and amplitude modulation(PAM), the present invention determines " sieve " stopping of process standard, and the value of the standard deviation SD by calculating double screening result determines whether " sieve " stops:
Satisfied 0.2<SD<0.3 o'clock, the process of first leg " sieve " finishes, and obtains first intrinsic mode function component c
1, this component is the highest frequency component in the original signal data, this condition has been controlled the number of times of " sieve ", makes the intrinsic mode function component that obtains keep the information of amplitude modulation(PAM) in the original signal data.
From the original signal data, isolate first intrinsic mode function component c
1, obtain residue signal r
1:
x(t)-c
1=r
1
r
1Therefore still comprise the frequency information in the original signal data, it is repeated the process of above-mentioned " sieve " as new signal, obtain second intrinsic mode function component c
2, by that analogy;
r
1-c
2=r
2
r
n-1-c
n=r
n
As the intrinsic mode function component c that obtains at last
nOr residue signal r
nValue less than in advance setting value 0.2~0.3, perhaps last residue signal r
nBe monotonic quantity, in the time of can not " sieve " goes out the intrinsic mode function component again, the empirical mode decomposition process finishes.Be that the original signal data are n intrinsic mode function component and a residual volume r by empirical mode decomposition
nThe original signal data can be expressed as:
The pattern discrimination of II, signal
This step is according to the energy distribution of each intrinsic mode function component, distinguish echo signal leading with intrinsic mode function component noise dominant.
For noise, no matter be single noise, or the different mixed noise of rank, the energy distribution of its intrinsic mode function component is unified successively decreasing all.Therefore according to this characteristic as can be known: in the signal data of Noise, as first intrinsic mode function component c
1Energy less than the energy of other intrinsic mode function component, the intrinsic mode function component of respectively organizing of monotone decreasing is noise.No matter that is to say, be single noise signal or mixed noise signal, and the energy distribution characteristics of its intrinsic mode function component are same mode.Just can judge that from the energy monotone decreasing pattern of the intrinsic mode function component of signal this group component is noise dominant.Its concrete steps are as follows:
I, intrinsic mode function divide the calculating of energy
To each intrinsic mode function component c
iCalculate the mould of its square, as its ENERGY E
i, computing formula is:
The pattern discrimination of ii, signal
As first intrinsic mode function component c
1ENERGY E
1More than or equal to the energy of the intrinsic mode function component of other intrinsic mode function component, i.e. E
1When being in a ratio of maximum, first intrinsic mode function component c
1Be the signal dominant component; Otherwise, first intrinsic mode function component c
1Be the noise dominant component.
When q intrinsic mode function component ENERGY E
qAmplitude is less than 2 * 10
18(μ V)
2The time, signal is very weak, is more or less the same with the size of noise, if as the deletion of noise dominant component, will also have been rejected the signal section that contains the deep information.So q reaches later intrinsic mode function component all as the useful signal dominant component, carries out noise reduction process by filtering.
Second to q-1 intrinsic mode function component c
iENERGY E
iWith Next intrinsic mode function component c
I+1ENERGY E
I+1Compare, work as E
i〉=E
I+1, be monotone decreasing.Judging has several monotone decreasing groups, and the first intrinsic mode function component in each monotone decreasing group is the noise dominant component, and other is the useful signal dominant component.
III, cancelling noise
Be judged as the deletion of noise dominant component in front q-1 intrinsic mode function component in Step II.
IV, filtering
The intrinsic mode function component leading to each signal of Step II I gained carries out the filtering of least mean-square error (LMS) Filtering and smoothing.
V, signal restructuring
Each intrinsic mode function component after step IV processed re-starts cumulative, obtains removing the signal behind the noise.
The advantage of the noise-reduction method of transient electromagnetic detecting echo signal of the present invention is: adopt the empirical mode decomposition of adaptive Time Frequency Analysis, the solid-state pattern based on signal need not any prior imformation.This law is processed the different natural frequencies composition of transient electromagnetic detecting echo signal, rejects the high frequency noise part, and low-frequency noise is partly carried out filtering, then recombinate, thereby obtain cleaner echoed signal, especially for the deep layer echoed signal, its anti-acoustic capability is very optimistic.Be applicable to the noise reduction process of the transient electromagnetic detecting echo signal of nonlinear and nonstationary, solved the extraction problem of weak signal under the strong noise background.
(4) description of drawings
Fig. 1 is the noise-reduction method embodiment general flow chart of transient electromagnetic detecting echo signal;
Fig. 2 is the process flow diagram of pattern discrimination of the noise-reduction method embodiment signal of transient electromagnetic detecting echo signal;
Fig. 3 is the noise-reduction method embodiment LMS sef-adapting filter structural representation of transient electromagnetic detecting echo signal;
Fig. 4 is the noise-reduction method embodiment transient electromagnetic detecting echo signal oscillogram of transient electromagnetic detecting echo signal; Its ordinate is that voltage, unit are microfarad, and horizontal ordinate is that time, unit are second;
Fig. 5 is the oscillogram after the noise-reduction method embodiment transient electromagnetic detecting echo signal of transient electromagnetic detecting echo signal adds very noisy; Coordinate is identical with Fig. 4 in length and breadth for it;
Fig. 6 is that the noise-reduction method embodiment of transient electromagnetic detecting echo signal adds the oscillogram of transient electromagnetic detecting echo signal after this law is processed of making an uproar, and coordinate is identical with Fig. 4 in length and breadth for it.
(5) embodiment
The main-process stream that this transient electromagnetic detecting receives the signal de-noising embodiment of the method as shown in Figure 1, step is as follows:
I, transient electromagnetic detecting is received signal carry out empirical mode decomposition (EMD)
All extreme points of picked up signal data x (t) at first, all local maximums are formed the coenvelope of data with cubic spline functions, in like manner, all local minimums are formed the lower envelope of data with cubic spline functions, upper lower envelope should cover all data points, and its average is denoted as m
1, from former data sequence, deduct this average and obtain first component h
1:
h
1=x(t)-m
1
Said process is called " sieve ".
In the process of the second time " sieve ", h
1Regard pending data as, its envelope average is m
11, then have
h
11=h
1-m
11
This process repeats k time,
h
1k=h
1(k-1)-m
1k
Until h
1kSatisfy till the condition of intrinsic mode function, from the original signal data, decomposite first intrinsic mode function component, be denoted as:
c
1=h
1k
Calculate the value of the standard deviation SD of double screening result:
Satisfied 0.2<SD<0.3 o'clock, the process of first leg " sieve " finishes, and obtains first intrinsic mode function component c
1, from the original signal data, isolate first intrinsic mode function component c
1, obtain residue signal r
1:
x(t)-c
1=r
1
r
1Still comprise the frequency information in the original signal data, it is repeated the process of above-mentioned " sieve " as new signal, obtain second intrinsic mode function component c
2, by that analogy;
r
1-c
2=r
2
r
n-1-c
n=r
n
This routine intrinsic mode function component c
8Less than setting value 0.2 in advance, the empirical mode decomposition process finishes; The original signal data are 8 intrinsic mode function components and a residual volume r by empirical mode decomposition
8The original signal data are expressed as:
If obtain first residue signal r
nValue less than in advance setting value 0.2, perhaps last residue signal r
nDuring for monotonic quantity, same empirical mode decomposition process finishes.
The pattern discrimination of II, signal
The flow process of this step as shown in Figure 2, its concrete steps are as follows:
I, intrinsic mode function divide the calculating of energy
To each intrinsic mode function component c
iCalculate the mould of its square, as its ENERGY E
i, computing formula is:
The pattern discrimination of ii, signal
As shown in Figure 2, as first intrinsic mode function component c
1Energy when being in a ratio of maximum with other the energy of intrinsic mode function component respectively, c
1Be the signal dominant component; Otherwise, c
1Be the noise dominant component; This routine c
1Be the signal dominant component;
The 5th intrinsic mode function component ENERGY E of this example
5Amplitude is less than 2 * 10
18(μ V)
2So the 5th to the 8th intrinsic mode function component is as the useful signal dominant component.
Second to the 4th intrinsic mode function component c
iENERGY E
iWith Next intrinsic mode function component c
I+1ENERGY E
I+1Compare, the first intrinsic mode function component in each monotone decreasing group is the noise dominant component, and other is the useful signal dominant component.
This example is E
2〉=E
3, E
3<E
4, E
4〉=E
5, be 2 groups of monotone decreasings, so c
2, c
4Be noise dominant component, c
3Be the signal dominant component.
III, cancelling noise
Be judged as the deletion of noise dominant component in front four intrinsic mode function components in Step II.This example obtains 6 of useful signal dominant component, is c
1, c
3And c
5To c
8
IV, filtering
The intrinsic mode function component leading to each signal of Step II I gained carries out the filtering of least mean-square error (LMS) Filtering and smoothing.
I, least mean-square error filtering LMS filtering
The LMS sef-adapting filter as shown in Figure 3, z among the figure
-1Be chronotron, W
I1(kk) be weight coefficient, p=1-M, xx (kk) is input, and y (kk) is output, and ε (kk) is square error, and the basic element of character that consists of adaptive digital filter is adaptive linear combiner.The intrinsic mode function component that the M that Step II I an obtains signal is dominated is xx (kk-1) as the input of LMS wave filter ... .xx (kk-M), the output y (kk) of LMS wave filter are the linear combination behind the above-mentioned weighted input, namely
The LMS algorithm is obtained the weight coefficient W hour as ε (kk)
i, the solution formula of weight coefficient is as follows:
W
p(kk)=W
p(kk-1)+2u[xx(kk)-W
T(kk-1)xx(kk)]xx(kk)
Wherein the speed of convergence of the weight coefficient of the size of noise and setting is determined in the selection of u, generally gets u=0.001.
The LMS algorithm is followed the tracks of the leading intrinsic mode function component of each signal, removes adaptively the noise of environment.
Ii, smothing filtering
Affected by noise may the shake a little of the filtered part IMF component of LMS, can carry out smothing filtering this moment to it.
Smothing filtering mainly comes this point is carried out the wave amplitude correction according to the wave amplitude of the contiguous sampled point of certain signaling point, thereby waveform is carried out denoising.Can do simply on average also can be weighted on average j neighbor point as required to j neighbor point, j be 3~9 integer.This example is got the data point of 5 vicinities, and the basic calculating formula of smothing filtering is as follows:
In the formula, v is sampled data; Yy is the data after the smoothing processing; J is number of data points; N is on average counting; Hh
NnBe the weighted mean factor.
Wherein weighting factor satisfies following formula:
This example adopts five-spot triple smoothing, and its formula is as follows:
In the formula: b=3,4 ... .j-2.
In order to obtain better result, can increase according to actual needs for the level and smooth number of times of each component.
V, signal restructuring
Each intrinsic mode function component after step IV processed re-starts cumulative, obtains removing the signal behind the noise.
Figure 4 shows that the echoed signal oscillogram of transient electromagnetic detecting standard; Fig. 5 is the oscillogram after transient electromagnetic detecting echo signal adds very noisy, and echoed signal is submerged in the noise waveform and is beyond recognition; Fig. 6 is that this law is removed the oscillogram of making an uproar after processing, and three figure see that more clearly the added most noise of transient electromagnetic detecting echo signal all by filtering, safeguarded the feature of echoed signal.This shows that method of the present invention can be for the noise reduction process of the transient electromagnetic detecting echo signal that contains very noisy.
Above-described embodiment only is the specific case that purpose of the present invention, technical scheme and beneficial effect are further described, and the present invention is defined in this.All any modifications of within scope of disclosure of the present invention, making, be equal to replacement, improvement etc., all be included within protection scope of the present invention.