CN108693555A - Intelligent time-varying blind deconvolution wideband processing method and processing device - Google Patents
Intelligent time-varying blind deconvolution wideband processing method and processing device Download PDFInfo
- Publication number
- CN108693555A CN108693555A CN201810465628.7A CN201810465628A CN108693555A CN 108693555 A CN108693555 A CN 108693555A CN 201810465628 A CN201810465628 A CN 201810465628A CN 108693555 A CN108693555 A CN 108693555A
- Authority
- CN
- China
- Prior art keywords
- value
- time
- varying
- wavelet
- values
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
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
Abstract
The present invention provides a kind of intelligent time-varying blind deconvolution wideband processing method and processing device, wherein this approach includes the following steps:Go out primary earthquake wavelet from decaying earthquake extracting data;Randomly generate the Q value groups for including multiple stratum Q values;Time-varying wavelet matrix group is established according to the primary earthquake wavelet and the Q value groups;According to the decaying seismic data and the time-varying wavelet matrix group, multiple reflectance factors are obtained;Establish fitness function;According to the fitness function and the multiple reflectance factor, multiple adaptive values are obtained, wherein an adaptive value corresponds to a stratum Q value;Optimally interval quality factors Q values are determined according to the multiple adaptive value.Since the program does not determine optimal Q values by manually, but optimal Q values are obtained by intelligent search, reduces workload and error, to estimate optimal time-varying wavelet from unstable state data, while obtaining optimal time-varying blind deconvolution result.
Description
Technical field
The present invention relates to geophysical exploration High-resolution Processing, inverting and explain technical field, more particularly to a kind of intelligence
Time-varying blind deconvolution wideband processing method and processing device can be changed.
Background technology
In field of geophysical exploration, the waveform data of acquisition can regard the wavelet of ground excitation as by ground
Addition of waveforms when being received after lower propagation and interface reflection, is the result of seismic wavelet and reflectance factor convolution.Theoretically its
It is required that the seismic data of input is stable state (decaying and the frequency dispersion effect that do not consider stratum quality factor q).If seismic wave
Its Fourier frequency spectrum (amplitude spectrum and phase spectrum) does not change in communication process, and it is stable state to claim this seismic wave.On the contrary, ground
Seismic wave its frequency spectrum in communication process changes (amplitude is decayed and phase frequency dispersion), and it is unstable state to claim this seismic wave.By
Unstable state seismic wave groups at data be known as unstable state seismic data.Due to the viscoplasticity on stratum, actual seismic data are non-steady
State.Decaying and the frequency dispersion of seismic wave, physical significance are usually described with quality factor q in actual seismic data handling procedure
It is seismic wave after propagating a wavelength distance, the ratio between the energy stored originally and the energy consumed.
Traditional reflectance factor inversion method is based on stable state wavelet it is assumed that it is stable state theoretically to require the seismic data of input
's.For unstable state seismic data, if stable state deconvolution or reflectance factor inversion method is still used to carry out High-resolution Processing,
It can lead to even more serious image artifacts, further influence subsequent seismic interpretation work.Then use such methods from practical field
Inverting obtains reflection coefficient sequence in the unstable state seismic data of acquisition, it is necessary first to compensate the Q filter effects on stratum.Anti- Q filters
Wave is more a kind of method of common compensation stratum Q filter effects, generally includes amplitude compensation and phase compensation.Inverse Q filtering
Amplitude compensation is unstable intrinsically, unless considering in earnest very much during inverse Q filtering and devising stable
Amplitude compensation operator, amplitude compensation item will inevitably amplify high-frequency noise.At other unstable state data high-resolution
For reason method, accurate time-varying wavelet matrix is the precondition for obtaining reasonable inversion result.In general, Q model is inclined
Difference often leads to seismic wavelet and is had differences between the attenuation degree and actual value of different time position, varitron when to influence
The accuracy of wave matrix is further solved by influencing Waveform Matching so that the reflectance factor of estimation generates deviation.Therefore, accurately
Estimate that Q model and structure inversion equation are to obtain truly to descend the precondition of reflected impulse image in ground.There are also sweep at present
The method for retouching Q strategies to obtain best Q value, but Q value effects rely primarily on artificial cognition, and this brings a lot for inverting work
Workload and error.The trend sequencing that current seism processing is explained is intelligent, a set of independent of people to be badly in need of
The intelligentized deconvolution system of work intervention.
Invention content
An embodiment of the present invention provides a kind of intelligent time-varying blind deconvolution wideband processing method and processing devices, not by artificial
It determines optimal Q values, but optimal Q values is obtained by intelligent search, reduce workload and error, to from unstable state data
Optimal time-varying wavelet is estimated, while obtaining optimal time-varying blind deconvolution result.
The intelligence time-varying blind deconvolution wideband processing method includes:
Go out primary earthquake wavelet from decaying earthquake extracting data;
Randomly generate the Q value groups for including multiple stratum Q values;
Time-varying wavelet matrix group is established according to the primary earthquake wavelet and the Q value groups;
According to the decaying seismic data and the time-varying wavelet matrix group, multiple reflectance factors are obtained;
Establish fitness function;
According to the fitness function and the multiple reflectance factor, multiple adaptive values are obtained, wherein an adaptive value pair
Answer a stratum Q value;
Optimally interval quality factors Q values are determined according to the multiple adaptive value.
The intelligence time-varying blind deconvolution wideband processing unit includes:
Primary earthquake wavelet extraction module, for going out primary earthquake wavelet from decaying earthquake extracting data;
Q value group generation modules, for randomly generating the Q value groups for including multiple stratum Q values;
Time-varying wavelet matrix group establishes module, varitron when for being established according to the primary earthquake wavelet and the Q value groups
Wave matrix group;
Reflectance factor obtains module, for according to the decaying seismic data and the time-varying wavelet matrix group, obtaining more
A reflectance factor;
Fitness function establishes module, for establishing fitness function;
Adaptive value obtains module, for according to the fitness function and the multiple reflectance factor, obtaining multiple adaptations
Value a, wherein adaptive value corresponds to a stratum Q value;
Optimally interval quality factors Q values determining module, for according to the multiple adaptive value determine optimal stratum quality because
Sub- Q values.
In embodiments of the present invention, optimally interval quality factors Q values are determined by the fitness function of foundation, and it is existing
Technology is compared, and does not determine optimal Q values by manually, but obtains optimal Q values by intelligent search, reduces workload and mistake
Difference to estimate optimal time-varying wavelet from unstable state data, while obtaining optimal time-varying blind deconvolution result.
Description of the drawings
In order to more clearly explain the embodiment of the invention or the technical proposal in the existing technology, to embodiment or will show below
There is attached drawing needed in technology description to be briefly described, it should be apparent that, the accompanying drawings in the following description is only this
Some embodiments of invention for those of ordinary skill in the art without creative efforts, can be with
Obtain other attached drawings according to these attached drawings.
Fig. 1 is a kind of intelligent time-varying blind deconvolution wideband process flow figure provided in an embodiment of the present invention;
Fig. 2 is a kind of intelligent time-varying blind deconvolution wideband processing method particular flow sheet provided in an embodiment of the present invention;
Fig. 3 is that true time-varying wavelet under the conditions of a kind of nothing provided in an embodiment of the present invention is made an uproar and 20 different randoms are initial
The time-varying wavelet that model inversion obtains;
Fig. 4 is that true and 20 different random initial model invertings obtain under the conditions of a kind of nothing provided in an embodiment of the present invention is made an uproar
The reflection coefficient sequence obtained;
Fig. 5 be 20 with decay characteristics true and synthesis under the conditions of a kind of nothing provided in an embodiment of the present invention is made an uproar not
With the earthquake record data of random initial model inverting reconstruct;
Fig. 6 is the true time-varying wavelet and 20 different randoms under the conditions of provided in an embodiment of the present invention a kind of noisy 5%
The time-varying wavelet that initial model inverting obtains;
Fig. 7 be under the conditions of one kind noisy 5% provided in an embodiment of the present invention true reflection coefficient sequence and 20 differences with
The reflection coefficient sequence that machine initial model inverting obtains;
Fig. 8 is true record and synthesis with decay characteristics under the conditions of provided in an embodiment of the present invention a kind of noisy 5%
20 different random initial model invertings reconstruct earthquake record data;
Fig. 9 is a kind of actual seismic data design sketch provided in an embodiment of the present invention;
Figure 10 is a kind of Q model result figure using obtained by method proposed by the present invention provided in an embodiment of the present invention;
Figure 11 is that a kind of application method proposed by the present invention provided in an embodiment of the present invention carries out actual seismic data
The design sketch of inverting;
Figure 12 is a kind of intelligent time-varying blind deconvolution wideband processing device structure diagram provided in an embodiment of the present invention.
Specific implementation mode
Following will be combined with the drawings in the embodiments of the present invention, and technical solution in the embodiment of the present invention carries out clear, complete
Site preparation describes, it is clear that described embodiment is only a part of the embodiment of the present invention, instead of all the embodiments.Based on this
Embodiment in invention, every other reality obtained by those of ordinary skill in the art without making creative efforts
Example is applied, shall fall within the protection scope of the present invention.
In embodiments of the present invention, a kind of intelligent time-varying blind deconvolution wideband processing method is provided, as shown in Figure 1,
This method includes:
Step 101:Go out primary earthquake wavelet from decaying earthquake extracting data;
Step 102:Randomly generate the Q value groups for including multiple stratum Q values;
Step 103:Time-varying wavelet matrix group is established according to the primary earthquake wavelet and the Q value groups;
Step 104:According to the decaying seismic data and the time-varying wavelet matrix group, multiple reflectance factors are obtained;
Step 105:Establish fitness function;
Step 106:According to the fitness function and the multiple reflectance factor, multiple adaptive values are obtained, wherein one
Adaptive value corresponds to a stratum Q value;
Step 107:Optimally interval quality factors Q values are determined according to the multiple adaptive value.
When it is implemented, assuming there is n Q value particles i=1,2 ..., n, each Q values particle in the Q value group bodies randomly generated
It is the individual in D dimension spaces, the position of each Q values particle is expressed as xi=(xi1,xi2,...,xiD), wherein dimension is
1,2, D.Each Q values particle moves in D dimension spaces, and speed is expressed as vi=(vi1,vi2,...,viD).No
Same Q values particle is in different positions in space, there is different adaptive values relative to fitness function f (r).It is good in group
Position, that is, adaptive value minimum particle, use PgIt indicates;The desired positions P that i-th of Q value particle lives throughiIt indicates.Such as Fig. 2 institutes
Show, then step 102 specifically includes:Randomly generate the initial position x of each stratum Q valueiWith initial velocity vi;It determines
The initial global desired positions P of the Q value groupsg, the initial history desired positions P of each stratum Q valuei。
When it is implemented, step 103 specifically establishes time-varying wavelet matrix group as follows.
From classical Acoustic Wave-equation, it is as follows to establish its plane wave analytic solutions:
Wherein, s (t) indicates that earthquake record and plane wave, ω indicate that frequency content, r (t ') represent reflectance factor,Monochromatic wave is represented, t '=h/v represents the propagation time, and h represents depth, v representation speeds.The analytic solutions can regard different as
One linear weighted function of frequency monochromatic wave is superimposed, and weight coefficient is reflectance factor.
Then the inner real speed of formula (1) is substituted for complex velocity, plane wave analytic solutions under attenuation medium can be established such as
Under:
Wherein,It indicates that amplitude attenuation term, γ indicate decay time spectrum, depends on depth time
Variable;Indicate phase correction terms or wavelet shaping item, ω0It is to refer to angular frequency, Q (t ') is medium
Quality factor.
Discretization is carried out to formula (2), the mathematical connection that can be established between data and Q values and reflectance factor is as follows:
Or,
S=FHWA(Q)r (4)
Wherein,Indicate Fourier's operator,Indicate primary earthquake Wavelet Martrix,
Indicate that the time-varying damping matrix dependent on Q values, s indicate that decaying seismic data, r indicate reflectance factor.
Time-varying wavelet matrix group G is defined to indicate the bandpass filtering effect of seismic wavelet and the Q filter effects on stratum.Time-varying
The expression formula of Wavelet Martrix group G is:
G=FHWA(Q) (5)
Wherein, the n-th row of time-varying wavelet matrix group G are the primary earthquake wavelet w at the n-th momentn;HIndicate transposition conjugation.
When it is implemented, step 104 specifically includes:According to the decaying seismic data and the time-varying wavelet matrix group,
Time-varying blind deconvolution processing is carried out with sparse Bayesian inversion algorithm, obtains multiple reflectance factors.
Specifically, according to formula (4) and (5), noisy unstable state data d can be write a Chinese character in simplified form into:
Gr+n=d (6)
Wherein, r=[r1,...,rK]TThe model parameter for waiting for inverting is represented, n represents noise, and T indicates transposition.Due to ground
The band limit essence of data is shaken, it is ill posed directly to solve equation (6).For this purpose, carrying out time-varying using sparse Bayesian inversion method
Blind anti-pleat processing.
Assuming that n obedience mean values are 0, variance σ2Gaussian Profile,That is, in setting models parameter r and σ2
Under conditions of, the conditional probability of d, also referred to as likelihood function can be expressed as:
p(d|r,σ22 π σ of)=(2)-Kexp[-(d-Gr)T(d-Gr)/(2σ2)] (7)
Assumed that driving is solved by geology to obtain, management loading method needs to apply prerequisite to r, according to
Bayesian theory, the probability distribution for limiting r are the standardized normal distributions centered on zero, i.e.,
Wherein, h=[h1,h2,L,hK]TK mutually independent hyper parameters are represented, it is right that each hyper parameter controls its respectively
Answer the prior information of reflectance factor size.Indicate reflectance factor rkIt is 0 to obey mean value, and variance isGauss point
Cloth, wherein k=1,2 ..., K.
By bayesian criterion, in d, h and σ2Under the conditions of known, then there are the item of reflectance factor r in simultaneous formula (7) and (8)
Part probability or posterior probability Density Distribution:
Wherein,
Wherein, C indicates a constant;Σ indicates covariance;μ indicates mean value;H=diag (h1,h2,L,hK) it is one right
Angular moment battle array.The estimation of reflectance factor is provided by the mean μ of reflectance factor Posterior distrbutionp.In order to obtain time-varying blind deconvolution as a result,
First have to estimation hyper parameter h and σ2Optimum value.According to Bayesian frame, the edge distribution of hyper parameter is calculate by the following formula
p(d|h,σ2∫ p (the d| of)=- 2;r,σ2)p(r|H) dr=(2 π)K|Q|exp(dTQ-1d) (11)
Wherein,I indicates unit matrix.To formula
(11) ask minimum that can find out hyper parameter h.Here, a kind of quick iterative algorithm-Method Using Relevance Vector Machine has been used to obtain time-varying
Blind deconvolution result.The each iteration of the algorithm only updates an a reflectance factor pulse r or base vector Gk, but every time after update
Object function or formula (11) can all reduce.It iterates, calculates h using formula (11), formula (10) calculates μ, until front and back 2
The difference of the object function of secondary iteration reaches tolerable error.Finally, the μ of output is the sparse reflectance factor of estimation.
Substantially, by formula (11), iteration finds out hyper parameter h every time, is equivalent to the position for obtaining reflection Sparse Pulse r
(position refers to that each non-zero reflection coefficient pulse corresponds to time location in reflectance factor, because there is many zeros in reflectance factor,
So being sparse), and μ is found out by formula (10), is equivalent to and obtains the size of reflection Sparse Pulse r (size refers to non-zero
The amplitude size of reflectance factor pulse).Since Method Using Relevance Vector Machine is in interative computation, simply by add or delete base vector Gk
Come position and the size of real-time update pulse, the inversion operation of big matrix need not be carried out, there are receipts quickly in this way
Hold back speed.
When it is implemented, step 105 specifically establishes fitness function as follows.
The inversion result of same decaying seismic data, different Q value is different, and accurately the corresponding inversion result of Q values is sparse
Spend higher.For this phenomenon, present invention lpNorm carries out the sparsity assessment of reflection coefficient sequence.
For the l of some sequencepNorm, the size of value is by the non-zero number and amplitude for measuring whole element in sequence
Collective effect.And with the reduction of p, the weight shared by element non-zero number increases, and amplitude weight reduces.Due to seismic signal
Waveform Matching condition, when the seismic wavelet waveform that inaccurate Q values generate removes the waveform in matching physical record, the difference of phase
It is different to lead to not portray the waveform using single reflectance factor pulse, but need multiple superimposed pulses that can fit
The signal known, therefore there is continuous same polarity pulse in obtained reflectance factor inversion result, and it is not sparse.Even amplitude
It is faint, still result in the l of entire inversion resultpNorm increases.Therefore, lpNorm can accurately reflect its degree of rarefication
Size, by searching so that lpThe minimum value of Norm function can determine real reflectance factor and the combination of Q values.When for building
When the Q of mapping matrix G is exact value, corresponding reflectance factor inversion result is most sparse, lpNorm minimum.
L is selected in the present invention0.1Norm is as interpretational criteria and establishes following fitness function
F (r)=s ||r||0.1 (12)
Wherein, r is the reflection coefficient sequence obtained by management loading inverting.
When it is implemented, step 107 specifically includes:
The multiple adaptive value is compared between any two, obtains minimum adaptive value;
The minimum adaptive value is compared with default adaptive value threshold value, when minimum adaptive value is less than default adaptive value threshold
When value, the corresponding stratum Q value of the minimum adaptive value is optimally interval quality factors Q values.
Step 107 is exactly when not finding, i.e., to be preset when minimum adaptive value is more than to find current globally optimal solution
When adaptive value threshold value, the Q value groups are updated;Optimally interval quality factors Q values are redefined according to the updated Q value groups.
As shown in Figure 2.
When minimum adaptive value is more than default adaptive value threshold value, then need to update the Q value groups, including:According to described every
Initial position and initial velocity, the initial global desired positions of the Q value groups, each stratum product of a stratum Q value
The initial history desired positions of prime factor Q values update the Q value groups.
Specifically, update is to want:Update the desired positions P of groupgAnd the desired positions P that each Q value particles live throughi;
Update the position x of groupiWith speed vi;
The position x of group is updated according to following formulaiWith speed vi:
vi(k+1)=ω vi(k)+c1φ1[pi(k)-xi(k)]+c2φ2[pg(k)-xi(k)] (13)
xi(k+1)=xi(k)+avi(k+1) (14)
Wherein, k indicates update times;vi(k+1) speed of i-th of Q value when+1 update of kth is indicated;ω indicates inertia
Coefficient;vi(k) speed of i-th of Q value when kth time update is indicated;pi(k) the best position of i-th of Q value when kth time update is indicated
It sets;xi(k) position of i-th of Q value when kth time update is indicated;pg(k) desired positions of Q value groups when kth time update are indicated;c1With
c2For aceleration pulse, between 0~2;φ1And φ2It is independent normal distyribution function, mean value 0, variance 1;xi(k+1) table
The position of i-th of Q value when showing+1 update of kth;A is constraint factor, constant.
As shown in Fig. 2, when the condition is satisfied, terminate update, export best reflectance factor as a result, one of condition just
It is that the minimum adaptive value mentioned above is less than default adaptive value threshold value.Another condition is exactly that update times reach in preset times
Limit, the effect of the condition are that update is prevented to be absorbed in endless loop.It is referred to specifically, meeting condition:(1) it is obtained when by repeatedly update
Minimum adaptive value be less than default adaptive value threshold value and when update times are not up to the preset times upper limit, can terminate to update;
(2) when update times have had reached the preset times upper limit, but pass through the minimum adaptive value that repeatedly update obtains and be also greater than
When default adaptive value threshold value, also to terminate to update;(3) when the minimum adaptive value by repeatedly update acquisition is less than default adaptive value
When threshold value and update times have reached the preset times upper limit, it can terminate to update.
Embodiment
Verification is handled seismic data using intelligent time-varying blind deconvolution wideband processing method proposed by the present invention
Effect.
It tests to obtain a seismic wavelet by a reflectance factor, then uses the method for the present invention to this seismic wavelet
Inversion procedure is carried out, obtains optimal time-varying wavelet, as shown in Figure 3.Fig. 3 indicate without make an uproar under the conditions of true time-varying wavelet and adopt
The time-varying wavelet obtained with 20 different random initial model invertings.Wavelet at where the dotted line signifies that superficial part 100ms, band are real
Heart solid line indicates the wavelet at the 500ms of middle part, and the wavelet at the 900ms of deep is indicated with asterisk solid line.It can be seen that true and anti-
Drill the wavelet of acquisition it is shallow, in, deep layer all it is matched very well, and can be obtained under the conditions of 20 different random initial models
It is such good as a result, demonstrating the stability and accuracy of time-varying blind deconvolution method proposed by the present invention.
Fig. 4 indicates true reflection coefficient sequence under the conditions of without making an uproar and uses 20 different random initial model invertings
The reflection coefficient sequence of acquisition.Wherein right side the 1st is true reflection coefficient sequence, and left side 20 is 20 different random
The reflection coefficient sequence that inverting obtains under the conditions of initial model.It can be seen that the reflectance factor and true reflectance factor of inverting are basic
It coincide, it is effective to show time-varying blind deconvolution method proposed by the present invention.
Fig. 5 indicates that true and synthesis 20 different random initial models with decay characteristics are anti-under the conditions of without making an uproar
Drill the earthquake record data of reconstruct.As can be seen that the Waveform Matching carried out using time-varying blind deconvolution method proposed by the present invention
Very well.
True time-varying wavelet and 20 different random initial model invertings of Fig. 6 expressions under the conditions of noisy 5% obtain
Time-varying wavelet.As can be seen that as depth increases, estimates that the precision of wavelet is opposite and be deteriorated, but still it is acceptable.
Fig. 7 indicates that true reflection coefficient sequence and 20 different random initial model invertings obtain under the conditions of noisy 5%
Reflection coefficient sequence.As can be seen that the inversion result that is obtained using time-varying blind deconvolution method proposed by the present invention and true
Model is almost the same.
Fig. 8 shows 20 different randoms of true record and synthesis that decay characteristics are carried under the conditions of noisy 5% are initial
The earthquake record data of model inversion reconstruct, it can be seen that 20 obtained using time-varying blind deconvolution method proposed by the present invention
The earthquake record data of a different random initial model inverting reconstruct, with the presence of some and true record fine difference, still
The two can be described as almost the same, and it is effective to show time-varying blind deconvolution method proposed by the present invention.
Fig. 9 is original seismic data.In view of the variation of section horizontal direction than shallower, carried out just for the therein 60th
Then handling result is applied to other roads by processing.In addition, section, is divided into 8 portions by the characteristics of according to section from top to bottom
Point, the Q value value ranges of the time location and corresponding portion at this 8 part interfaces are given, are calculated.Wherein initial wavelet
It is to estimate from the earthquake record of preceding 200ms.Figure 10 is the Q value curves finally estimated.Figure 11 is the result of inverting.It can see
It arrives, it can be seen that the corresponding data information in this part 400-600ms is weaker in original section.400- in inversion result
The corresponding amplitude in this part 600ms has obtained effective compensation, and 600ms some detailed information below are all shown.It is whole same
Phase axis attenuates, and resolution ratio is got higher, and reflectance factor maintains preferable continuity in the horizontal.
Based on same inventive concept, a kind of intelligent time-varying blind deconvolution wideband processing is additionally provided in the embodiment of the present invention
Device, as described in the following examples.The principle and intelligence solved the problems, such as due to intelligent time-varying blind deconvolution wideband processing unit
It can change that time-varying blind deconvolution wideband processing method is similar, therefore the implementation of intelligent time-varying blind deconvolution wideband processing unit can be with
Referring to the implementation of intelligent time-varying blind deconvolution wideband processing method, overlaps will not be repeated.It is used below, term
The combination of the software and/or hardware of predetermined function may be implemented in " unit " or " module ".Although described in following embodiment
Device is preferably realized with software, but the realization of the combination of hardware or software and hardware is also that may and be contemplated.
Figure 12 is a kind of structure diagram of the intelligent time-varying blind deconvolution wideband processing unit of the embodiment of the present invention, is such as schemed
Shown in 12, including:
Primary earthquake wavelet extraction module 1201, for going out primary earthquake wavelet from decaying earthquake extracting data;
Q value groups generation module 1202, for randomly generating the Q value groups for including multiple stratum Q values;
Time-varying wavelet matrix group establishes module 1203, when for being established according to the primary earthquake wavelet and the Q value groups
Become wavelet matrix group;
Reflectance factor obtains module 1204, for according to the decaying seismic data and the time-varying wavelet matrix group, obtaining
Obtain multiple reflectance factors;
Fitness function establishes module 1205, for establishing fitness function;
Adaptive value obtains module 1206, for according to the fitness function and the multiple reflectance factor, obtaining multiple
Adaptive value a, wherein adaptive value corresponds to a stratum Q value;
Optimally interval quality factors Q values determining module 1207, for determining optimal stratum product according to the multiple adaptive value
Prime factor Q values.
The structure is illustrated below.
It is specifically used for when it is implemented, the reflectance factor obtains module 1204:
According to the decaying seismic data and the time-varying wavelet matrix group, when being carried out with sparse Bayesian inversion algorithm
Become blind deconvolution processing, obtains multiple reflectance factors.
It is specifically used for when it is implemented, the fitness function establishes module 1205:
Select l0.1Norm establishes fitness function as interpretational criteria.
When it is implemented, the optimally interval quality factors Q values determining module 1207 is specifically used for:
The multiple adaptive value is compared between any two, obtains minimum adaptive value;
The minimum adaptive value is compared with default adaptive value threshold value, when minimum adaptive value is less than default adaptive value threshold
When value, the corresponding stratum Q value of the minimum adaptive value is optimally interval quality factors Q values.
When it is implemented, the intelligence time-varying blind deconvolution wideband processing unit further includes:
Update module, for when minimum adaptive value is more than default adaptive value threshold value, updating the Q value groups;
The optimally interval quality factors Q values determining module 1207 is specifically used for:
Optimally interval quality factors Q values are redefined according to the updated Q value groups.
When it is implemented, the Q value groups generation module 1202 is specifically used for:
Randomly generate the initial position and initial velocity of each stratum Q value;
Determine the initial global desired positions of the Q value groups, the best position of initial history of each stratum Q value
It sets;
The update module is specifically used for:
According to the initial position of each stratum Q value and initial velocity, the Q value groups it is initial it is global most
The initial history desired positions of good position, each stratum Q value update the Q value groups.
When it is implemented, the update module is additionally operable to:
When update times reach the preset times upper limit or minimum adaptive value is less than default adaptive value threshold value, terminate update.
It is specifically used for when it is implemented, the time-varying wavelet matrix group establishes module 1203:
Time-varying wavelet matrix group is established according to following formula:
G=FHWA(Q);
Wherein, G indicates time-varying wavelet matrix group;F indicates Fourier's operator;W indicates primary earthquake Wavelet Martrix;A (Q) table
Show the time-varying damping matrix dependent on stratum Q value;H indicates transposition conjugation.
When it is implemented, the update module is specifically used for:
The Q value groups are updated according to following formula:
vi(k+1)=ω vi(k)+c1φ1[pi(k)-xi(k)]+c2φ2[pg(k)-xi(k)];
xi(k+1)=xi(k)+avi(k+1);
Wherein, k indicates update times;vi(k+1) speed of i-th of Q value when+1 update of kth is indicated;ω indicates inertia
Coefficient;vi(k) speed of i-th of Q value when kth time update is indicated;pi(k) the best position of i-th of Q value when kth time update is indicated
It sets;xi(k) position of i-th of Q value when kth time update is indicated;pg(k) desired positions of Q value groups when kth time update are indicated;c1With
c2For aceleration pulse, between 0~2;φ1And φ2It is independent normal distyribution function, mean value 0, variance 1;xi(k+1) table
The position of i-th of Q value when showing+1 update of kth;A is constraint factor, constant.
In conclusion find under study for action, when Q values are not punctual, had differences between inversion result and true model, mainly
It shows in pulse amplitude and non-zero pulses number.And this phenomenon of inversion result how to be utilized to capture true Q model, it is this
A key of swarm intelligence optimization fitness function is established in invention.For the inverting knot of same decaying seismic data different Q value
Fruit is different, while accurately there is the corresponding inversion result of Q values this phenomenon of higher degree of rarefication, the present invention to select l0.1Norm
Fitness function is established as interpretational criteria and carries out evaluation of result.When the Q for building time-varying wavelet matrix is exact value,
Corresponding management loading reflectance factor inversion result is most sparse, l0.1The fitness function of norm structure is minimum.Cause
This, the present invention uses the strategy of particle group optimizing method intelligentized search Q, it is proposed that intelligent time-varying blind deconvolution wideband processing
Method and system can not only estimate optimal time-varying wavelet under the premise of meeting wave propagation law from unstable state data,
And optimal time-varying blind deconvolution result can be obtained simultaneously.Therefore, improve existing (time-varying) higher-order spectra and (when
Become) efficiency and precision of blind anti-pleat technology.
It should be understood by those skilled in the art that, the embodiment of the present invention can be provided as method, system or computer program
Product.Therefore, complete hardware embodiment, complete software embodiment or reality combining software and hardware aspects can be used in the present invention
Apply the form of example.Moreover, the present invention can be used in one or more wherein include computer usable program code computer
The computer program production implemented in usable storage medium (including but not limited to magnetic disk storage, CD-ROM, optical memory etc.)
The form of product.
The present invention be with reference to according to the method for the embodiment of the present invention, the flow of equipment (system) and computer program product
Figure and/or block diagram describe.It should be understood that can be realized by computer program instructions every first-class in flowchart and/or the block diagram
The combination of flow and/or box in journey and/or box and flowchart and/or the block diagram.These computer programs can be provided
Instruct the processor of all-purpose computer, special purpose computer, Embedded Processor or other programmable data processing devices to produce
A raw machine so that the instruction executed by computer or the processor of other programmable data processing devices is generated for real
The device for the function of being specified in present one flow of flow chart or one box of multiple flows and/or block diagram or multiple boxes.
These computer program instructions, which may also be stored in, can guide computer or other programmable data processing devices with spy
Determine in the computer-readable memory that mode works so that instruction generation stored in the computer readable memory includes referring to
Enable the manufacture of device, the command device realize in one flow of flow chart or multiple flows and/or one box of block diagram or
The function of being specified in multiple boxes.
These computer program instructions also can be loaded onto a computer or other programmable data processing device so that count
Series of operation steps are executed on calculation machine or other programmable devices to generate computer implemented processing, in computer or
The instruction executed on other programmable devices is provided for realizing in one flow of flow chart or multiple flows and/or block diagram one
The step of function of being specified in a box or multiple boxes.
The foregoing is only a preferred embodiment of the present invention, is not intended to restrict the invention, for the skill of this field
For art personnel, the embodiment of the present invention can have various modifications and variations.All within the spirits and principles of the present invention, made by
Any modification, equivalent substitution, improvement and etc. should all be included in the protection scope of the present invention.
Claims (10)
1. a kind of intelligence time-varying blind deconvolution wideband processing method, which is characterized in that including:
Go out primary earthquake wavelet from decaying earthquake extracting data;
Randomly generate the Q value groups for including multiple stratum Q values;
Time-varying wavelet matrix group is established according to the primary earthquake wavelet and the Q value groups;
According to the decaying seismic data and the time-varying wavelet matrix group, multiple reflectance factors are obtained;
Establish fitness function;
According to the fitness function and the multiple reflectance factor, multiple adaptive values are obtained, wherein an adaptive value corresponds to one
A stratum Q value;
Optimally interval quality factors Q values are determined according to the multiple adaptive value.
2. intelligence time-varying blind deconvolution wideband processing method as described in claim 1, which is characterized in that according to the decaying
Seismic data and the time-varying wavelet matrix group, obtain multiple reflectance factors, including:
According to the decaying seismic data and the time-varying wavelet matrix group, cecutiency when being carried out with sparse Bayesian inversion algorithm
Deconvolution is handled, and obtains multiple reflectance factors.
3. intelligence time-varying blind deconvolution wideband processing method as described in claim 1, which is characterized in that establish fitness letter
Number, including:
Select l0.1Norm establishes fitness function as interpretational criteria.
4. intelligence time-varying blind deconvolution wideband processing method as described in claim 1, which is characterized in that according to the multiple
Adaptive value determines optimally interval quality factors Q values, including:
The multiple adaptive value is compared between any two, obtains minimum adaptive value;
The minimum adaptive value is compared with default adaptive value threshold value, when minimum adaptive value is less than default adaptive value threshold value
When, the corresponding stratum Q value of the minimum adaptive value is optimally interval quality factors Q values.
5. intelligence time-varying blind deconvolution wideband processing method as claimed in claim 4, which is characterized in that further include:
When minimum adaptive value is more than default adaptive value threshold value, the Q value groups are updated;
Optimally interval quality factors Q values are redefined according to the updated Q value groups.
6. intelligence time-varying blind deconvolution wideband processing method as claimed in claim 5, which is characterized in that randomly generate and include
The Q value groups of multiple stratum Q values, including:
Randomly generate the initial position and initial velocity of each stratum Q value;
Determine the initial global desired positions of the Q value groups, the initial history desired positions of each stratum Q value;
When minimum adaptive value is more than default adaptive value threshold value, the Q value groups are updated, including:
According to the initial best position of the overall situation of the initial position of each stratum Q value and initial velocity, the Q value groups
It sets, the initial history desired positions of each stratum Q value, updates the Q value groups.
7. such as intelligent time-varying blind deconvolution wideband processing method described in claim 5 or 6, which is characterized in that further include:
When update times reach the preset times upper limit and/or minimum adaptive value is less than default adaptive value threshold value, terminate update.
8. intelligence time-varying blind deconvolution wideband processing method as described in claim 1, which is characterized in that according to primary earthquake
Wavelet and stratum quality factor q value group establish time-varying wavelet matrix group, including:
G=FHWA(Q);
Wherein, G indicates time-varying wavelet matrix group;F indicates Fourier's operator;W indicates primary earthquake Wavelet Martrix;A (Q) indicate according to
Rely the time-varying damping matrix in Q values;HIndicate transposition conjugation.
9. intelligence time-varying blind deconvolution wideband processing method as claimed in claim 6, which is characterized in that according to described each
Initial global desired positions, each stratum quality of the initial position and initial velocity of stratum Q value, the Q value groups
The initial history desired positions of factor Q value update the Q value groups according to following formula:
vi(k+1)=ω vi(k)+c1φ1[pi(k)-xi(k)]+c2φ2[pg(k)-xi(k)];
xi(k+1)=xi(k)+avi(k+1);
Wherein, k indicates update times;vi(k+1) speed of i-th of Q value when+1 update of kth is indicated;ω indicates inertia coeffeicent;
vi(k) speed of i-th of Q value when kth time update is indicated;pi(k) desired positions of i-th of Q value when kth time update are indicated;xi
(k) position of i-th of Q value when kth time update is indicated;pg(k) desired positions of Q value groups when kth time update are indicated;c1And c2For
Aceleration pulse, between 0~2;φ1And φ2It is independent normal distyribution function, mean value 0, variance 1;xi(k+1) the is indicated
The position of i-th of Q value when k+1 update;A is constraint factor, constant.
10. a kind of intelligence time-varying blind deconvolution wideband processing unit, which is characterized in that including:
Primary earthquake wavelet extraction module, for going out primary earthquake wavelet from decaying earthquake extracting data;
Q value group generation modules, for randomly generating the Q value groups for including multiple stratum Q values;
Time-varying wavelet matrix group establishes module, for establishing time-varying wavelet square according to the primary earthquake wavelet and the Q value groups
Battle array group;
Reflectance factor obtains module, for according to the decaying seismic data and the time-varying wavelet matrix group, obtaining multiple anti-
Penetrate coefficient;
Fitness function establishes module, for establishing fitness function;
Adaptive value obtains module, for according to the fitness function and the multiple reflectance factor, obtaining multiple adaptive values,
In, an adaptive value corresponds to a stratum Q value;
Optimally interval quality factors Q values determining module, for determining optimally interval quality factors Q according to the multiple adaptive value
Value.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810465628.7A CN108693555B (en) | 2018-05-16 | 2018-05-16 | Intelligent time-varying blind deconvolution wideband processing method and processing device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810465628.7A CN108693555B (en) | 2018-05-16 | 2018-05-16 | Intelligent time-varying blind deconvolution wideband processing method and processing device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108693555A true CN108693555A (en) | 2018-10-23 |
| CN108693555B CN108693555B (en) | 2019-08-02 |
Family
ID=63847547
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810465628.7A Active CN108693555B (en) | 2018-05-16 | 2018-05-16 | Intelligent time-varying blind deconvolution wideband processing method and processing device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108693555B (en) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110988985A (en) * | 2019-12-18 | 2020-04-10 | 北京邮电大学 | Seismic signal detection method based on waveform characteristics |
| CN111060961A (en) * | 2019-12-27 | 2020-04-24 | 中国石油大学(北京) | Quality factor determination method, device and system based on multi-information constraint inversion |
| CN111880218A (en) * | 2020-07-13 | 2020-11-03 | 西南石油大学 | Inversion wavelet dictionary construction method based on quality factor |
| CN112415587A (en) * | 2019-08-21 | 2021-02-26 | 中国石油化工股份有限公司 | Reservoir seismic wave attenuation characteristic analysis method and reservoir reflection coefficient inversion method |
| CN112666603A (en) * | 2019-10-16 | 2021-04-16 | 中国石油化工股份有限公司 | Lp norm constraint-based phase identification method and system |
| CN113219530A (en) * | 2020-02-05 | 2021-08-06 | 中国石油天然气股份有限公司 | Unsteady-state blind deconvolution method and device |
| CN113341463A (en) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | Pre-stack seismic data non-stationary blind deconvolution method and related components |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090204330A1 (en) * | 2005-02-18 | 2009-08-13 | Leon Thomsen | System and method for using time-distance characteristics in acquisition, processing, and imaging of t-csem data |
| CN106324675A (en) * | 2016-10-09 | 2017-01-11 | 中国石油大学(华东) | Broad earthquake wave impedance low-frequency information prediction method and system |
| CN107356964A (en) * | 2017-07-05 | 2017-11-17 | 西安交通大学 | Q value estimation and compensation method of the S-transformation domain based on variation principle |
-
2018
- 2018-05-16 CN CN201810465628.7A patent/CN108693555B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090204330A1 (en) * | 2005-02-18 | 2009-08-13 | Leon Thomsen | System and method for using time-distance characteristics in acquisition, processing, and imaging of t-csem data |
| CN106324675A (en) * | 2016-10-09 | 2017-01-11 | 中国石油大学(华东) | Broad earthquake wave impedance low-frequency information prediction method and system |
| CN107356964A (en) * | 2017-07-05 | 2017-11-17 | 西安交通大学 | Q value estimation and compensation method of the S-transformation domain based on variation principle |
Non-Patent Citations (1)
| Title |
|---|
| 程亮 等: "步进迭代法井地联合地震资料拓频处理", 《石油地球物理勘探》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112415587A (en) * | 2019-08-21 | 2021-02-26 | 中国石油化工股份有限公司 | Reservoir seismic wave attenuation characteristic analysis method and reservoir reflection coefficient inversion method |
| CN112666603A (en) * | 2019-10-16 | 2021-04-16 | 中国石油化工股份有限公司 | Lp norm constraint-based phase identification method and system |
| CN110988985A (en) * | 2019-12-18 | 2020-04-10 | 北京邮电大学 | Seismic signal detection method based on waveform characteristics |
| CN110988985B (en) * | 2019-12-18 | 2020-11-17 | 北京邮电大学 | Seismic signal detection method based on waveform characteristics |
| CN111060961A (en) * | 2019-12-27 | 2020-04-24 | 中国石油大学(北京) | Quality factor determination method, device and system based on multi-information constraint inversion |
| CN111060961B (en) * | 2019-12-27 | 2020-11-20 | 中国石油大学(北京) | Quality factor determination method, device and system based on multi-information constraint inversion |
| CN113219530A (en) * | 2020-02-05 | 2021-08-06 | 中国石油天然气股份有限公司 | Unsteady-state blind deconvolution method and device |
| CN113219530B (en) * | 2020-02-05 | 2023-09-26 | 中国石油天然气股份有限公司 | Unsteady state blind deconvolution method and device |
| CN111880218A (en) * | 2020-07-13 | 2020-11-03 | 西南石油大学 | Inversion wavelet dictionary construction method based on quality factor |
| CN113341463A (en) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | Pre-stack seismic data non-stationary blind deconvolution method and related components |
| CN113341463B (en) * | 2021-06-10 | 2023-05-26 | 中国石油大学(北京) | Non-stationary blind deconvolution method for pre-stack seismic data and related components |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108693555B (en) | 2019-08-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108693555B (en) | Intelligent time-varying blind deconvolution wideband processing method and processing device | |
| US8892413B2 (en) | Convergence rate of full wavefield inversion using spectral shaping | |
| US11360223B2 (en) | System and method for improved full waveform inversion | |
| EP3129809B1 (en) | Seismic adaptive focusing | |
| CN105425289B (en) | The method and apparatus for determining low frequency wave impedance | |
| US20180164453A1 (en) | Method for Improved Geophysical Investigation | |
| CN103293552A (en) | Pre-stack seismic data retrieval method and system | |
| CN113945982B (en) | Method and system for removing low frequency and low wave number noise to generate enhanced images | |
| CN110895348B (en) | Method, system and storage medium for extracting low-frequency information of seismic elastic impedance | |
| Robinson et al. | Digital imaging and deconvolution: The ABCs of seismic exploration and processing | |
| EP3548932A2 (en) | Seismic acquisition geometry evaluation using full-waveform inversion | |
| Dhara et al. | Physics-guided deep autoencoder to overcome the need for a starting model in full-waveform inversion | |
| Ivanov et al. | Traveltime parameters in tilted orthorhombic medium | |
| CN104730572A (en) | Diffracted wave imaging method and device based on L0 semi-norm | |
| CN112711065A (en) | Pre-stack seismic inversion method and device | |
| CN108680968A (en) | Complex structural area seismic prospecting data collecting observation system evaluation method and device | |
| Kazei et al. | Mapping full seismic waveforms to vertical velocity profiles by deep learning | |
| Greer et al. | Improving migration resolution by approximating the least-squares Hessian using nonstationary amplitude and frequency matching | |
| CN114114421B (en) | Deep learning-based guided self-learning seismic data denoising method and device | |
| CN109857979A (en) | A kind of log approximating method based on wavelet analysis and BP neural network | |
| CN111781635A (en) | Seabed four-component elastic wave Gaussian beam depth migration method and device | |
| US9383464B2 (en) | Seismic imaging apparatus without edge reflections and method for the same | |
| CN112415601A (en) | Method and device for determining surface quality factor Q value | |
| Zhu | Seismic modeling, inversion, and imaging in attenuating media | |
| Santos et al. | Fast Marchenko multiples elimination on CMP processing |
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 |