CN109613609B - A kind of combination signal decomposition method based on elongated degree processing - Google Patents
A kind of combination signal decomposition method based on elongated degree processing Download PDFInfo
- Publication number
- CN109613609B CN109613609B CN201910038470.XA CN201910038470A CN109613609B CN 109613609 B CN109613609 B CN 109613609B CN 201910038470 A CN201910038470 A CN 201910038470A CN 109613609 B CN109613609 B CN 109613609B
- Authority
- CN
- China
- Prior art keywords
- data
- continuation
- value
- indicate
- point
- 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.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
Abstract
The present invention relates to a kind of combination signal decomposition methods based on elongated degree processing, belong to seismic data process field, specifically includes the following steps: the data that one of length to be processed is N are extracted from original seismic data first, each resection length step in the length L and calculating process of the continuation of data-oriented both ends, the mode of both ends continuation is judged according to the numerical values recited at data endpoint and establishes new extension data, and empirical mode decomposition is carried out to new data；After often obtaining single order intrinsic mode function, the data of step length are cut away at both ends, data continue empirical mode decomposition after being then based on excision；Simultaneously when each rank intrinsic mode function calculates, the extreme value point data needed for envelope calculates is made of the extreme point of extension data after cutting off and its extreme point two parts after continuation.The present invention can effectively suppress the end effect problem in empirical mode decomposition, be of great significance to the precision that seismic signals decompose is improved.
Description
Technical field
The invention belongs to seismic data process fields, and in particular to a kind of combination signal point based on elongated degree processing
Solution method.
Background technique
Due to the heterogeneity of subsurface formations medium, vertical, horizontal variation, seismic wavelet after Different Strata decays,
The characteristics such as frequency, amplitude and phase are not quite similar.Nonlinear and non local boundary value problem, signal are belonged to for general seismic data
Feature shows as time-varying and block overlap of frequency bands more, this causes useful signal to be frequently found in some or certain several frequency bands, therefore
Separating and extract useful signal from original signal is one of important research content of seismic data process.Empirical mode decomposition
It (EMD) is a kind of signal spectrum decomposition method proposed by Chinese American HuangE. et al. in 1998.The core concept of this method
It is that the signal of any complexity can be broken down into the sum of limited intrinsic mode function (IMF) and a residual components.Because
This method is the decomposition realized on the basis based on signal feature itself, it has compared with the method based on Fourier transformation
Stronger adaptivity and Time Frequence Analysis ability, therefore with very high practical in analyzing non-linear and non-stationary signal
Value.
Although empirical mode decomposition is deemed appropriate to decompose non-linear and non-stationary signal, due to its propose time compared with
It is short, lack theory support, there are many problem to be solved, such as signal boundary is easy to generate in application experience mode decomposition
Signal decomposition problem of dtmf distortion DTMF, i.e. end effect.Its producing cause is when carrying out signal decomposition based on EMD method, due to endpoint
Place is non-extreme point, and will appear at endpoint when calculating signal envelope cannot complete envelope and as caused by cubic spline interpolation
The problems such as overshoot and owe punching, this meeting is so that the IMF function error after decomposing is excessive, so that the accuracy to signal decomposition result produces
Raw serious influence.
Domestic and foreign scholars have been presented for the conventional method of a variety of compacting end effects at present, such as characteristic wave method, Waveform Matching
Method, mirror image circumference symmetric extension method, extreme point end effect method, continuation method neural network based etc., but these methods exist
There is a degree of deficiency in.For example, the problem of characteristic wave method, is that the both ends of non-stationary nonlinear properties have
Wave character it is not necessarily identical, so selecting suitable characteristic wave relatively difficult；Waveform Matching method calculates in practical applications
Efficiency is relatively low；Mirror image circumference symmetric extension rule directly handles endpoint as extreme point, iterates and calculated
Cheng Zhonghui causes the distortion of endpoint value；Extreme point end effect method is usually using the extreme point nearest apart from endpoint as mirror point pair
Signal carries out continuation, also needs to carry out truncation, if data are shorter, end effect to signal when endpoint is not extreme point
Compacting can fail；That there are computational efficiencies is low and to big data quantity and signal period for continuation method neural network based
The problems such as dependence.
The problem of for current empirical mode decomposition method, needs to propose more effectively solve empirical modal
The method of end effect problem in decomposition, enables the needs for meeting seismic data process.
Summary of the invention
The technical problem to be solved in the present invention is that providing a kind of combination signal decomposition method based on elongated degree processing.?
Seismic data is carried out in EMD decomposable process, doing end extending to signal is essentially all to make by the feature of legacy data
Prediction, for non-linear non-stable seismic data, if the variation at endpoint has no rule and changes acutely, by various again
The result that prediction technique obtains may all be futile.Due to predict can not entirely accurate, at the aft terminal of possible successive ignition
Value dissipate, and will inwardly pollute entire data.In consideration of it, the present invention proposes the outside class symmetric extension one at endpoint
Point initial data makes it become to extend data, and extracts intrinsic mode function to extending data and carrying out empirical mode decomposition.Often
To after single order intrinsic mode function, the data of certain length are cut away at extension data both ends, this segment data is forecasting inaccuracy
True data make this inaccuracy not to be substituted into the solution of lower single order intrinsic mode function, thus utmostly after excision
Guarantee M signal be not contaminated, be then based on excision after extend data continue empirical mode decomposition；Simultaneously in each rank
When intrinsic mode function calculates, the extreme value point data needed for envelope calculates is by the extreme point of extension data and its continuation after cutting off
Extreme point two parts composition afterwards.Empirical Mode is realized by way of the processing of the elongated degree of data and continuation method combination in this way
Data extension at the endpoint that state is decomposed, is a kind of effective end effect compacting means.
The present invention takes following technical scheme:
A kind of combination signal decomposition method based on elongated degree processing, it is characterised in that including
1) setting man-made explosion generates seismic wave, receives reflection seismic waves by wave detector in earth's surface, obtains earthquake note
Record；If initial data is a certain section of signal to be processed in earthquake record, data length N, initial data is expressed as xi,
Middle i=0,1,2 ..., N-1；
2) length for setting the continuation of data both ends is L, and the maximum order of intrinsic mode function is M；Establish new data
To extend data, i=0,1,2 ..., N+2*L-1 are indicated initial data xiBoth ends respectively extend the data after L length, then
To new dataThe point of middle corresponding i=L, L+1, L+2 ..., N+L-1 carry out assignment, i.e.,
3) according to initial data xiThe size of data judges the continuation mode at data endpoint and to new data at endpoint
Assignment is carried out at endpoint, by taking left end point as an example:
Wherein, i=0,1,2 ..., L-1, ximinIndicate the minimum value of initial data, ximaxIndicate the maximum of initial data
Value, x0Indicate initial data xiValue at i=0, δ ∈ [0,1] are weight coefficient；Work as x0Value between ximinWith ximaxBetween when,
Point does class central symmetry continuation centered on the point；Work as x0Amplitude it is larger when, do symmetric extension using the point as extreme point；Right end
Continuation mode can similarly obtain at point；
4) to extension dataIt carries out intrinsic mode function in empirical mode decomposition to seek, is calculating every single order natural mode
During state function, it is both needed to prolong the extreme point position at data both ends and extreme value data before calculating extreme point envelope every time
It opens up, the extreme value point data needed for such envelope calculates is by the extreme point of extension data after cutting off and its extreme point two after continuation
Part forms；
5) when seeking out I rank intrinsic mode functionLater, wherein I indicates current intrinsic mode function order, initially
Value is 1 and is less than maximum order to be M；It indicates to extend dataI rank intrinsic mode function, i=0,1,2 ..., N+2*L-
1, fromI=L, L+1, L+2 ... are designated as under middle selection, the data of N+L-1 are assigned to IMFi I, i=0,1 ..., N-1, wherein
IMFi IIndicate initial data xiI rank intrinsic mode function；
6) the I rank intrinsic mode function that will be calculated in step 5)Number before calculating the intrinsic mode function
According toIn cut, obtain extend data redundancy componentExtending data redundancy component
Both ends the data of step length be individually subtracted obtainWhereinIt indicates to extend data redundancy componentBoth ends be individually subtracted
Data after step length, step indicate the data length cut away, value step=L/M, with season L=L-step；With
As extension dataContinue to calculate m+1 rank intrinsic mode function, then for calculating the extension of m+1 rank intrinsic mode function
DataIts length shortens 2*step compared with when calculating m rank intrinsic mode function；
4) 5) 6) step is repeated, until empirical mode decomposition terminates and from extending data redundancy componentI=L, L+ are designated as under middle selection
The data of 1, L+2 ..., N+L-1 are assigned to initial data residual components；
Further, the step of continuation being carried out to the extreme point position at data both ends and extreme value data described in step 4) packet
It includes:
(1) from the extension data described in step 4)Extreme value point data in extract the n information being worth at the proximal points,
Position data including n extreme point is lk, k=1,2 ..., n, extreme value data are vk, k=1,2 ..., n；The m pole to continuation
The position data of value point is lk, k=n+1, n+2 ..., n+m, extreme value data are vk, k=n+1, n+2 ..., n+m, wherein n is indicated
Given information point number, m indicate the number to continuation point, lkFor position data, vkFor extreme value data；
(2) here with lkContinuation calculate for, vkContinuation calculate and lkIt is identical；Non- continuation initial data is defined first
For xk (0), k=1,2,3 ..., n, wherein xk (0)Value and lk, the value correspondence of k=1,2,3 ..., n；If the number to continuation point is
M, the data after defining continuation areWork as k=n+1, when n+2 ..., n+m,Value and number
According to lk, k=n+1, n+2 ..., the value correspondence of n+m；
(3) to data initialization, xminIndicate xk (0)Minimum value, xmaxIndicate xk (0)Maximum value,WhereinIndicate initialization initial data, k=1,2,3 ..., n；
(4) it utilizesIt calculatesWherein k=1,2 ..., n,Indicate Accumulation Initialization original number
According to；
(5) it utilizesIt calculatesWherein k=2,3 ..., n,Indicate mean data；
(6) basisMatrix B and Y are established, wherein a indicates mean data system undetermined
It counts, the nonhomogeneous undetermined coefficient of b expression, the expression formula of matrix B and Y are
Wherein B indicates that coefficient matrix, Y indicate constant term；B and Y is substituted intoLeast square solution is calculated to ask
A and b is obtained, whereinIndicate initial data parameter vector,
(7) willSubstitute into equation
Calculating acquiresHereIt indicates to continuation value Accumulation Initialization initial data；
(8) basisIt acquiresHereIt indicates to initialize to continuation value former
Beginning data are rightIt is counter to be initialized, it calculatesHereIt indicates to continuation value initial data,Calculating formula meet following formula, thus obtain
(9) initial data x is calculatedk (0)WithError information φ, and the average calculation error residual errorError information φ
Calculating formula be
WhereinΔ(0)Indicate absolute error data andMeet
Mean error residual error
(10) α is given, whenAnd φnWhen < α is set up, error meets condition；WhenOr φn> α is calculated and is missed
Difference dataWherein
(11) initial data is initialized according to calculatingIt is original to continuation value Accumulation Initialization
DataProcess, by errorRegard asAccording to step 4) to asking in 7)
Take Accumulation Initialization initial dataTo continuation value Accumulation Initialization initial dataAnd initial data parameter vectorSide
Method respectively obtains and errorCorresponding Accumulation Initialization error informationTo continuation value Accumulation Initialization
Error informationAnd error information parameter vector WhereinIndicate cumulative initial
Change error information,To continuation value Accumulation Initialization error information, aeIndicate error mean data to be provided coefficient, beIndicate error
The nonhomogeneous undetermined coefficient of data,Calculating formula be
(12) utilization is found outIt is rightValue optimize, calculating formula is
δ in formulaeIt is correction weights coefficient, a indicates mean data undetermined coefficient, and b indicates nonhomogeneous undetermined coefficient, aeIt indicates
Error mean data to be provided coefficient, beIndicate the nonhomogeneous undetermined coefficient of error information,Indicate the cumulative just to continuation value of optimization
Beginningization initial data, k=0,1,2,3 ..., n；After acquiring optimization using above formula9) and 10) then step is repeated, directly
Meet condition to error；
(13) it is acquired using formula (6)And according toIt calculates Indicate optimization to continuation
Value initialization initial data, k=0,1,2 ..., n are right further according to formula (3)It is counter to be initialized, it obtainsValue, this
InIndicate optimization to continuation value initial data, k=0,1,2 ..., n；It, will as k=nValue be assigned to
According to the corresponding relationship in step (2), to l in the position data of continuationn+1Value calculating finish；Similarly, work as k=n+1, n+
When 2 ..., n+m-1, according to step (3) to (13), to l in the position data of continuationk+1Value pass through { lk-n+1,lk-n+2,…,
lkCalculated as initial data；Position data l to continuationk, k=n+1, n+2 ..., m value in n+m is extreme value
Value after point position Data extension；Extreme value data are v after similarly calculating continuationk, k=n+1, n+2 ..., the value of n+m, therefore,
Data l to continuationkAnd vk, k=n+1, n+2 ..., n+m by known n point data lkAnd vk, k=1,2 ..., n are by prolonging
It opens up and is calculated.
The beneficial effect of the present invention compared with prior art
The invention proposes a kind of combination signal decomposition methods based on elongated degree processing.It is non-due to subsurface formations medium
Homogenieity, vertical, horizontal variation, seismic wavelet is after Different Strata decays, the characteristics such as frequency, amplitude and phase phase not to the utmost
Together.Therefore seismic signal is usually expressed as nonlinear and non local boundary value problem.The present invention proposes the outside class symmetric extension one at endpoint
Part initial data makes it become to extend data, and carries out empirical mode decomposition extraction intrinsic mode function to data are extended.Often
After obtaining single order intrinsic mode function, the data of certain length are cut away at extension data both ends, this segment data is to predict not
Accurate data make this inaccuracy not to be substituted into the solution of lower single order intrinsic mode function, thus maximum journey after excision
The guarantee M signal of degree is not contaminated, and extension data continue empirical mode decomposition and extract intrinsic after being then based on excision
Mode function；Simultaneously when each rank intrinsic mode function calculates, the extreme value point data needed for envelope calculates is by extending after cutting off
The extreme point of data and its extreme point two parts composition after continuation.Pass through the elongated degree processing of data and continuation method in this way
Combined mode realizes Data extension at the endpoint of empirical mode decomposition, is a kind of effective end effect compacting means.
For nonlinear and non local boundary value problem, the method for the present invention is decomposed compared with the empirical mode decomposition method based on conventional continuation method
Error is obviously reduced, more acurrate to the processing of phase information at signal end, and discomposing effect is significantly improved.Therefore the present invention
It can be realized more accurate seismic data resolution process, this is established for the processing and inversion interpretation work of latter earthquake exploration
Good basis.
Detailed description of the invention
Fig. 1 is initial data xiWaveform, i=0,1,2 ..., 3000；
Fig. 2 is initial data xi(three figure lines are respectively the waveform of each component data from top to bottomWith )；
Fig. 3 be do not do end effect processing standard empirical mode decomposition method decomposition result (solid line be this method handle
As a result, IMF1, IMF2, IMF3 and IMF4 respectively decompose obtained the 1st, 2,3 rank intrinsic mode function of initial data and residue
Component)；
Fig. 4 is to carry out end effect processing using data end extending described in step 3), does not do the step 6) data and cuts
Cutting the decomposition method decomposition result that processing and step (1)-(13) the extreme point continuation are handled, (solid line is this method processing knot
Fruit, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components；Dotted line
Respectively correspond initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000)；
Fig. 5 is to be carried out using the data cutting of step 6) in data end extending described in step 3) and claim 1
End effect processing, do not do step (1)-(13) the extreme point continuation decomposition method decomposition result (solid line be this method at
Reason is as a result, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components；
Dotted line respectively corresponds initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000)；
Fig. 6 is to carry out end effect processing using step (1)-(13) the extreme point continuation, does not do the step 3) number
According to end extending handle and step 6) the data cutting process decomposition method decomposition result (solid line be this method processing knot
Fruit, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components；Dotted line
For initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000)；
Fig. 7 is to carry out endpoint using the extreme point continuation of data end extending described in step 3) and step (1)-(13)
Effect processing, do not do step 6) the data cutting process decomposition method decomposition result (solid line be this method processing result,
IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components；Dotted line is figure
Initial data x shown in 2iEach component dataWithWaveform, i=0,1,2 ..., 3000)；
For a kind of combination signal decomposition method decomposition result based on elongated degree processing, (solid line is this method processing knot to Fig. 8
Fruit, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components；Dotted line
For initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000)；
Specific embodiment
Technical solution of the present invention is further explained below by embodiment combination attached drawing, but protection of the invention
Range is not limited by any form.
Embodiment 1
The technical problem to be solved in the present invention is that providing a kind of combination signal decomposition method based on elongated degree processing.Tool
Body implementation process mainly includes two steps: 1) outside a part of initial data of class symmetric extension at endpoint makes it become new
Data cut away the data of certain length at both ends after often obtaining an intrinsic mode function, this segment data be prediction
The data of inaccuracy make this inaccuracy not to be substituted into the solution of lower single order intrinsic mode function, thus maximum after excision
The guarantee M signal of degree is not contaminated；2) when each rank intrinsic mode function calculates, the extreme value needed for envelope calculates is counted
According to being to carry out continuation by the extreme point to data after excision in step 1) to obtain.
Specific implementation process of the invention is described below:
1) setting man-made explosion generates seismic wave, receives reflection seismic waves by wave detector in earth's surface, obtains earthquake note
Record.It is that the seismic wave of earth's surface, t ∈ [0,3000], x (t) are reflected into through stratum after artificial epicenter excitation that x (t) is enabled in the present embodiment
Expression formula be
X (t)=15sin ((t/80)1.5π)+3sin((t/100)0.8π)+t/500+5, t ∈ [0,3000]
Since the working method of wave detector is discrete sampling, seismic reflection to earth's surface is detected what device received
Earthquake record is the point data of discretization.If the sample frequency of signal is 1Hz, data length N=3001, then x (t) is in earthquake
The value for being expressed as t in record is t=0, and the value of the series of discrete point of 1,2 ..., N-1 is denoted as xi, i=0,1,2 ..., N-
1, x hereiFor initial data, and as the data for needing resolution process, waveform is as shown in Figure 1；Under ideal conditions, x
(t) two sinusoidal component x are broken down into1(t)=15sin ((t/80)1.5π) and x2(t)=3sin ((t/100)0.8π) and one
A linear residual components x3(t)=t/500+5.To each component x1(t)、x2(t) and x3(t) each point is respectively obtained after discretization
Measuring data isAndThe waveform of each component data is as shown in Figure 2.
2) length for setting the continuation of data both ends is L, and the maximum order of intrinsic mode function is M, in the present embodiment
Length value is that 300, M value is 4；Establish new dataTo extend data, i=0,1,2 ..., N+2*L-1, expression will be former
Beginning data xiBoth ends respectively extend the data after L length, then to new dataMiddle corresponding i=L, L+1, L+2 ..., the point of N+L-1
Assignment is carried out, i.e.,
3) according to initial data xiThe size of data judges the continuation mode at data endpoint and to new data at endpoint
Assignment is carried out at endpoint, by taking left end point as an example:
Wherein, i=0,1,2 ..., L-1, ximinIndicate the minimum value of initial data, ximaxIndicate the maximum of initial data
Value, x0Indicate initial data xiValue at i=0, δ ∈ [0,1] are weight coefficient；Work as x0Value between ximinWith ximaxBetween when,
Point does class central symmetry continuation centered on the point；Work as x0Amplitude it is larger when, do symmetric extension using the point as extreme point；Right end
Continuation mode can similarly obtain at point；
4) to extension dataIt carries out intrinsic mode function in empirical mode decomposition to seek, is calculating every single order natural mode
During state function, it is both needed to prolong the extreme point position at data both ends and extreme value data before calculating extreme point envelope every time
It opens up, the extreme value point data needed for such envelope calculates is by the extreme point of extension data after cutting off and its extreme point two after continuation
Part forms；
5) when seeking out I rank intrinsic mode functionLater, wherein I indicates current intrinsic mode function order, initially
Value is 1 and is less than maximum order to be M；It indicates to extend dataI rank intrinsic mode function, i=0,1,2 ..., N+2*L-
1, fromI=L, L+1, L+2 ... are designated as under middle selection, the data of N+L-1 are assigned to IMFi I, i=0,1 ..., N-1, wherein
IMFi IIndicate initial data xiI rank intrinsic mode function；
6) the I rank intrinsic mode function that will be calculated in step 5)Number before calculating the intrinsic mode function
According toIn cut, obtain extend data redundancy componentExtending data redundancy component's
The data that step length is individually subtracted in both ends obtainWhereinIt indicates to extend data redundancy componentBoth ends be individually subtracted
Data after step length, step indicate the data length cut away, value step=L/M, with season L=L-step；With
As extension dataContinue to calculate m+1 rank intrinsic mode function, then for calculating the extension of m+1 rank intrinsic mode function
DataIts length shortens 2*step compared with when calculating m rank intrinsic mode function；
4) 5) 6) step is repeated, until empirical mode decomposition terminates and from extending data redundancy componentI=L, L+ are designated as under middle selection
The data of 1, L+2 ..., N+L-1 are assigned to initial data residual components；
Wherein include: to the step of extreme point position at data both ends and progress continuation of extreme value data described in step 4)
(1) from the extension data described in step 4)Extreme value point data in extract the n information being worth at the proximal points,
Position data including n extreme point is lk, k=1,2 ..., n, extreme value data are vk, k=1,2 ..., n；The m pole to continuation
The position data of value point is lk, k=n+1, n+2 ..., n+m, extreme value data are vk, k=n+1, n+2 ..., n+m, wherein n is indicated
Given information point number, m indicate the number to continuation point, lkFor position data, vkFor extreme value data；N value is 6 in the present embodiment,
M value is 4；
(2) here with lkContinuation calculate for, vkContinuation calculate and lkIt is identical；Non- continuation initial data is defined first
For xk (0), k=1,2,3 ..., n, wherein xk (0)Value and lk, the value correspondence of k=1,2,3 ..., n；If the number to continuation point is
M, the data after defining continuation areWork as k=n+1, when n+2 ..., n+m,Value and number
According to lkMiddle k=n+1, n+2 ..., the value of n+m point is corresponding；
(3) to data initialization, xminIndicate xk (0)Minimum value, xmaxIndicate xk (0)Maximum value,WhereinIndicate initialization initial data, k=1,2,3 ..., n；
(4) it utilizesIt calculatesWherein k=1,2 ..., n,Indicate Accumulation Initialization original number
According to；
(5) it utilizesIt calculatesWherein k=2,3 ..., n,Indicate mean data；
(6) basisMatrix B and Y are established, wherein a indicates mean data system undetermined
It counts, the nonhomogeneous undetermined coefficient of b expression, the expression formula of matrix B and Y are
Wherein B indicates that coefficient matrix, Y indicate constant term；B and Y is substituted intoLeast square solution is calculated to ask
A and b is obtained, whereinIndicate initial data parameter vector,
(7) willSubstitute into equation
Calculating acquiresHereIt indicates to continuation value Accumulation Initialization initial data；
(8) basisIt acquiresHereIt indicates to initialize to continuation value former
Beginning data are rightIt is counter to be initialized, it calculatesHereIt indicates to continuation value initial data,Calculating formula meet following formula, thus obtain
(9) initial data x is calculatedk (0)WithError information φ, and the average calculation error residual errorError information φ
Calculating formula be
WhereinΔ(0)Indicate absolute error data andMeet
Mean error residual error
(10) α is given, whenAnd φnWhen < α is set up, error meets condition；WhenOr φn> α is calculated and is missed
Difference dataWherein
(11) initial data is initialized according to calculatingTo continuation value Accumulation Initialization original number
According toProcess, by errorRegard asAccording to step 4) to seeking in 7)
Accumulation Initialization initial dataTo continuation value Accumulation Initialization initial dataAnd initial data parameter vectorMethod
It respectively obtains and errorCorresponding Accumulation Initialization error informationIt is missed to continuation value Accumulation Initialization
Difference dataAnd error information parameter vectorWhereinIndicate cumulative initial
Change error information,To continuation value Accumulation Initialization error information, aeIndicate error mean data to be provided coefficient, beIndicate error
The nonhomogeneous undetermined coefficient of data,Calculating formula be
(12) utilization is found outIt is rightValue optimize, calculating formula is
δ in formulaeIt is correction weights coefficient, a indicates mean data undetermined coefficient, and b indicates nonhomogeneous undetermined coefficient, aeIt indicates
Error mean data to be provided coefficient, beIndicate the nonhomogeneous undetermined coefficient of error information,Indicate the cumulative just to continuation value of optimization
Beginningization initial data, k=0,1,2,3 ..., n；After acquiring optimization using above formula9) and 10) then step is repeated, directly
Meet condition to error；
(13) it is acquired using formula (6)And according toIt calculates Indicate optimization to continuation
Value initialization initial data, k=0,1,2 ..., n are right further according to formula (3)It is counter to be initialized, it obtainsValue, hereIndicate optimization to continuation value initial data, k=0,1,2 ..., n；It, will as k=nValue be assigned toRoot
According to the corresponding relationship in step (2), to l in the position data of continuationn+1Value calculating finish；Similarly, work as k=n+1, n+2 ...,
When n+m-1, according to step (3) to (13), to l in the position data of continuationk+1Value pass through { lk-n+1,lk-n+2,…,lkMake
It is calculated for initial data；Position data l to continuationk, k=n+1, n+2 ..., m value in n+m is extreme value point
Value after setting Data extension；Extreme value data are v after similarly calculating continuationk, k=n+1, n+2 ..., the value of n+m, therefore, wait prolong
The data l opened upkAnd vk, k=n+1, n+2 ..., n+m by known n point data lkAnd vk, k=1,2 ..., n pass through continuation meter
It obtains.
By above sequence of steps, i.e., using a kind of combination signal point based on elongated degree processing of the present invention
Solution method decomposes the seismic data of extraction shown in FIG. 1.
In order to illustrate the necessity of end effect processing, signal shown in FIG. 1 is based on not doing end effect processing here
Standard empirical mode decomposition method handled, quadravalence intrinsic mode function is obtained in signal one after decomposition, as shown in figure 3,
Wherein IMF1, IMF2, IMF3 and IMF4 respectively decompose the 1st, 2,3 obtained rank intrinsic mode functions and residual components；Comparison
Suppressed if Fig. 2 and Fig. 3 can be seen that without end effect, obtained intrinsic mode function and residual components with it is original
There are significant differences for the signal component of data.
Since the present invention now compares and answers in order to illustrate the advantage of method after combination for a kind of combination signal decomposition method
Effect differently:
1. carrying out end effect processing using data end extending described in step 3), step 6) the data cut place is not done
The extreme point continuation processing of reason and step (1)-(13), decomposition result is as shown in figure 4, wherein solid line is this method processing knot
Fruit, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components, dotted line
Respectively correspond initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000；
2. carrying out endpoint using the data cutting of step 6) in data end extending described in step 3) and claim 1
Effect processing, does not do step (1)-(13) the extreme point continuation, and decomposition result is as shown in figure 5, wherein solid line is at this method
Reason as a result, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components,
Dotted line respectively corresponds initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000；
3. carrying out end effect processing using step (1)-(13) the extreme point continuation, the step 3) data terminal is not done
Point continuation processing and step 6) the data cutting process, decomposition result is as shown in fig. 6, wherein solid line is this method processing knot
Fruit, IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components, dotted line
Respectively correspond initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000；
4. carrying out endpoint effect using the extreme point continuation of data end extending described in step 3) and step (1)-(13)
Should handle, not do step 6) the data cutting process, decomposition result as shown in fig. 7, wherein solid line be this method processing result,
IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank intrinsic mode functions and residual components, dotted line difference
Corresponding initial data x shown in Fig. 2iEach component dataWithWaveform, i=0,1,2 ..., 3000；
5. the combination decomposition method carries out end effect processing according to the method for the present invention, decomposition result as shown in figure 8, its
Middle solid line is this method processing result, and IMF1, IMF2 and IMF3 respectively decompose obtained initial data the 1st, 2 rank natural mode of vibration
Function and residual components, dotted line respectively correspond initial data x shown in Fig. 2iEach component dataWithWaveform, i=
0,1,2,…,3000；
End effect processing method used in above method comparison is as shown in table 1.Distinct methods are provided in table 1 simultaneously
In the effect of two aspects of the waveform goodness of fit and error.Wherein, error calculation formula is as follows:
S indicates error, IMF in formulai IWithRespectively indicate initial data I rank intrinsic mode function and i-th number of components
According to the value at i point, when I is equal to 3, IMFi IIt indicates to decompose obtained initial data residual components.
1 distinct methods treatment effect of table compares
As it can be seen from table 1 the calculating error of the method for the invention decomposition result is 123.1, error is significantly less than it
His method.In addition, the method for the invention handles essence to data phase at endpoint in terms of data phase degree of agreement at endpoint
Degree is high, good with data phase degree of agreement at original signal endpoint.Thus explanation passes through the elongated degree processing of data and continuation side
The mode of method combination realizes Data extension at the endpoint of empirical mode decomposition, is a kind of effective end effect compacting hand
Section.
Claims (2)
1. a kind of combination signal decomposition method based on elongated degree processing, it is characterised in that the method includes
1) setting man-made explosion generates seismic wave, receives reflection seismic waves by wave detector in earth's surface, obtains earthquake record；If
Initial data is a certain section of signal to be processed in earthquake record, and data length N, initial data is expressed as xi, wherein i=
0,1,2,…,N-1；
2) length for setting the continuation of data both ends is L, and the maximum order of intrinsic mode function is M；Establish new dataTo extend
Data, i=0,1,2 ..., N+2*L-1 are indicated initial data xiBoth ends respectively extend the data after L length, then to extension
DataThe point of middle corresponding i=L, L+1, L+2 ..., N+L-1 carry out assignment, i.e.,
3) according to initial data xiAt endpoint the size of data judge the continuation mode at data endpoint and to new data endpoint at
Assignment is carried out, by taking left end point as an example:
Wherein, i=0,1,2 ..., L-1, ximinIndicate the minimum value of initial data, ximaxIndicate the maximum value of initial data, x0
Indicate initial data xiValue at i=0, δ ∈ [0,1] are weight coefficient；Work as x0Value between ximinWith ximaxBetween when, with this
Point does class central symmetry continuation centered on point；Work as x0Amplitude it is larger when, do symmetric extension using the point as extreme point；At right endpoint
Continuation mode can similarly obtain；
4) to extension dataIt carries out intrinsic mode function in empirical mode decomposition to seek, is calculating every single order natural mode of vibration letter
In several processes, it is both needed to carry out continuation to the extreme point position at data both ends and extreme value data before calculating extreme point envelope every time,
Extreme value point data needed for envelope calculates in this way is by the extreme point of extension data after cutting off and its extreme point two after continuation
It is grouped as；
5) when seeking out I rank intrinsic mode functionLater, wherein I indicates current intrinsic mode function order, initial value 1
And being less than maximum order is M；It indicates to extend dataI rank intrinsic mode function, i=0,1,2 ..., N+2*L-1, fromI=L, L+1, L+2 ... are designated as under middle selection, the data of N+L-1 are assigned to IMFi I, i=0,1 ..., N-1, wherein IMFi I
Indicate initial data xiI rank intrinsic mode function；
6) the I rank intrinsic mode function that will be calculated in step 5)Data before calculating the intrinsic mode function
In cut, obtain extend data redundancy componentExtending data redundancy componentBoth ends
The data that step length is individually subtracted obtainWhereinIt indicates to extend data redundancy componentBoth ends step long is individually subtracted
Data after degree, step indicate the data length cut away, value step=L/M, with season L=L-step；WithAs prolonging
Long dataContinue to calculate m+1 rank intrinsic mode function, then for calculating the extension data of m+1 rank intrinsic mode functionIts length shortens 2*step compared with when calculating m rank intrinsic mode function；It repeats
4) 5) 6) step, until empirical mode decomposition terminates and from extending data redundancy componentI=L, L+1, L+ are designated as under middle selection
The data of 2 ..., N+L-1 are assigned to initial data residual components.
2. a kind of combination signal decomposition method based on elongated degree processing, it is characterised in that data two described in the step 4)
The extreme point position at end and extreme value data carry out the step of continuation and include:
(1) from the extension data described in step 4)Extreme value point data in extract the n information being worth at the proximal points, including n
The position data of a extreme point is lk, k=1,2 ..., n, extreme value data are vk, k=1,2 ..., n；The m extreme point to continuation
Position data be lk, k=n+1, n+2 ..., n+m, extreme value data are vk, k=n+1, n+2 ..., n+m, wherein known to n expression
Information point number, m indicate the number to continuation point, lkFor position data, vkFor extreme value data；
(2) here with lkContinuation calculate for, vkContinuation calculate and lkIt is identical；Defining non-continuation initial data first is xk (0), k=1,2,3 ..., n, wherein xk (0)Value and lk, the value correspondence of k=1,2,3 ..., n；If the number to continuation point is m, fixed
Data after adopted continuation areWork as k=n+1, when n+2 ..., n+m,Value and data
lk, k=n+1, n+2 ..., the value correspondence of n+m；
(3) to data initialization, xminIndicate xk (0)Minimum value, xmaxIndicate xk (0)Maximum value,Its
InIndicate initialization initial data, k=1,2,3 ..., n；
(4) it utilizesIt calculatesWherein k=1,2 ..., n,Indicate Accumulation Initialization initial data；
(5) it utilizesIt calculatesWherein k=2,3 ..., n,Indicate mean data；
(6) basisMatrix B and Y are established, wherein a indicates mean data undetermined coefficient, b table
Show nonhomogeneous undetermined coefficient, the expression formula of matrix B and Y are
Wherein B indicates that coefficient matrix, Y indicate constant term；B and Y is substituted intoIt calculates least square solution and acquires a
And b, whereinIndicate initial data parameter vector,
(7) willSubstitute into equation
Calculating acquiresHereIt indicates to continuation value Accumulation Initialization initial data；
(8) basisIt acquiresHereIt indicates to initialize original number to continuation value
According to rightIt is counter to be initialized, it calculatesHereIt indicates to continuation value initial data,
Calculating formula meet following formula, thus obtain
(9) initial data x is calculatedk (0)WithError information φ, and the average calculation error residual errorThe meter of error information φ
Formula is
Wherein△(0)Indicate absolute error data andMeetIt is average
Error residual error
(10) α is given, whenAnd φnWhen < α is set up, error meets condition；WhenOr φn> α calculates error informationWherein
(11) initial data is initialized according to calculatingTo continuation value Accumulation Initialization initial dataProcess, by errorRegard asAccording to step 4) to seeking tiring out in 7)
Add initialization initial dataTo continuation value Accumulation Initialization initial dataAnd initial data parameter vectorMethod point
It does not obtain and errorCorresponding Accumulation Initialization error informationTo continuation value Accumulation Initialization error
DataAnd error information parameter vector WhereinIndicate that Accumulation Initialization misses
Difference data,To continuation value Accumulation Initialization error information, aeIndicate error mean data to be provided coefficient, beIndicate error information
Nonhomogeneous undetermined coefficient,Calculating formula be
(12) utilization is found outIt is rightValue optimize, calculating formula is
δ in formulaeIt is correction weights coefficient, a indicates mean data undetermined coefficient, and b indicates nonhomogeneous undetermined coefficient, aeIndicate error
Mean data undetermined coefficient, beIndicate the nonhomogeneous undetermined coefficient of error information,Indicate optimization to continuation value Accumulation Initialization
Initial data, k=0,1,2,3 ..., n；After acquiring optimization using above formula9) and 10) then step is repeated, until accidentally
Difference meets condition；
(13) it is acquired using formula (6)And according toIt calculates Indicate optimization at the beginning of continuation value
Beginningization initial data, k=0,1,2 ..., n are right further according to formula (3)It is counter to be initialized, it obtainsValue, here
Indicate optimization to continuation value initial data, k=0,1,2 ..., n；It, will as k=nValue be assigned toAccording to step
Suddenly the corresponding relationship in (2), to l in the position data of continuationn+1Value calculating finish；Similarly, work as k=n+1, n+2 ..., n+m-
When 1, according to step (3) to (13), to l in the position data of continuationk+1Value pass through { lk-n+1,lk-n+2,…,lkAs former
Beginning data are calculated；Position data l to continuationk, k=n+1, n+2 ..., m value in n+m is extreme point positional number
According to the value after continuation；Extreme value data are v after similarly calculating continuationk, k=n+1, n+2 ..., the value of n+m, therefore, to continuation
Data lkAnd vk, k=n+1, n+2 ..., n+m by known n point data lkAnd vk, k=1,2 ..., n are calculated by continuation.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910038470.XA CN109613609B (en) | 2019-01-16 | 2019-01-16 | A kind of combination signal decomposition method based on elongated degree processing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910038470.XA CN109613609B (en) | 2019-01-16 | 2019-01-16 | A kind of combination signal decomposition method based on elongated degree processing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109613609A CN109613609A (en) | 2019-04-12 |
| CN109613609B true CN109613609B (en) | 2019-09-24 |
Family
ID=66017658
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910038470.XA Active CN109613609B (en) | 2019-01-16 | 2019-01-16 | A kind of combination signal decomposition method based on elongated degree processing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109613609B (en) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101852871A (en) * | 2010-05-25 | 2010-10-06 | 南京信息工程大学 | Short-term climate forecasting method based on empirical mode decomposition and numerical value set forecasting |
| CN103496625A (en) * | 2013-10-17 | 2014-01-08 | 太原理工大学 | Multi-rope friction lifter load identification method based on vibration analysis |
| US8660848B1 (en) * | 2010-08-20 | 2014-02-25 | Worcester Polytechnic Institute | Methods and systems for detection from and analysis of physical signals |
| CN105116442A (en) * | 2015-07-24 | 2015-12-02 | 长江大学 | Lithologic oil-gas reservoir weak-reflection seismic signal reconstruction method |
| CN105760347A (en) * | 2016-02-04 | 2016-07-13 | 福建工程学院 | HHT end effect restraining method based on data/extreme value joint symmetric prolongation |
| CN108845230A (en) * | 2018-06-22 | 2018-11-20 | 国网陕西省电力公司电力科学研究院 | A kind of sub-synchronous oscillation random time-dependent modal identification method |
-
2019
- 2019-01-16 CN CN201910038470.XA patent/CN109613609B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101852871A (en) * | 2010-05-25 | 2010-10-06 | 南京信息工程大学 | Short-term climate forecasting method based on empirical mode decomposition and numerical value set forecasting |
| US8660848B1 (en) * | 2010-08-20 | 2014-02-25 | Worcester Polytechnic Institute | Methods and systems for detection from and analysis of physical signals |
| CN103496625A (en) * | 2013-10-17 | 2014-01-08 | 太原理工大学 | Multi-rope friction lifter load identification method based on vibration analysis |
| CN105116442A (en) * | 2015-07-24 | 2015-12-02 | 长江大学 | Lithologic oil-gas reservoir weak-reflection seismic signal reconstruction method |
| CN105760347A (en) * | 2016-02-04 | 2016-07-13 | 福建工程学院 | HHT end effect restraining method based on data/extreme value joint symmetric prolongation |
| CN108845230A (en) * | 2018-06-22 | 2018-11-20 | 国网陕西省电力公司电力科学研究院 | A kind of sub-synchronous oscillation random time-dependent modal identification method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109613609A (en) | 2019-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Boashash | Estimating and interpreting the instantaneous frequency of a signal. II. Algorithms and applications | |
| CN107272062B (en) | A kind of Q estimation methods of underground medium of data-driven | |
| Weedon | The detection and illustration of regular sedimentary cycles using Walsh power spectra and filtering, with examples from the Lias of Switzerland | |
| CN108267784A (en) | A kind of seismic signal random noise compression process method | |
| CN110674841B (en) | Logging curve identification method based on clustering algorithm | |
| CN104793253B (en) | Aviation electromagnetic data de-noising method based on mathematical morphology | |
| CA2456459C (en) | Method for enhancing depth and spatial resolution of one and two dimensional residual surfaces derived from scalar potential data | |
| Goldberg et al. | The response of simple nonlinear systems to a random disturbance of the earthquake type | |
| CN103956756A (en) | Electric system low-frequency oscillating mode identification method | |
| CN104933235B (en) | A kind of method for merging coastal waters multi-satellite sea level height abnormal data | |
| CN109884697A (en) | Glutenite sedimentary facies earthquake prediction method based on complete overall experience mode decomposition | |
| CN109613609B (en) | A kind of combination signal decomposition method based on elongated degree processing | |
| Wilson | Deep water wave generations by moving wind systems | |
| CN108427146A (en) | A kind of method and system of Milankovitch Cycles in identification sedimentary formation | |
| CN109633781B (en) | Geological property acquisition method and device, electronic equipment and storage medium | |
| CN109856672B (en) | Transient wave packet extracting method, storage medium and terminal based on depth wave-number spectrum | |
| CN108535774A (en) | A kind of method and device quickly identifying seismic phase using controlled source excitation seismic signal | |
| Yokoyama et al. | Sixty year variation in a time series of the geomagnetic Gauss coefficients between 1910 and 1983 | |
| CN111626377A (en) | Lithofacies identification method, device, equipment and storage medium | |
| CN109444957B (en) | Hydrate reservoir earthquake prediction method and processing terminal based on decision tree evaluation | |
| CN106443771B (en) | Improve the method and its velocity inversion method of converted wave seismic data resolution | |
| CN113050191B (en) | Shale oil TOC prediction method and device based on double parameters | |
| Kitagawa et al. | Signal extraction problems in seismology | |
| CN114298106A (en) | Characteristic wave identification method in roadbed rolling, rolling state discrimination method and application thereof | |
| Pan et al. | An approach to reserve estimation enhanced with 3-D seismic data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |