CN107367760B - Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm - Google Patents
Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm Download PDFInfo
- Publication number
- CN107367760B CN107367760B CN201710525466.7A CN201710525466A CN107367760B CN 107367760 B CN107367760 B CN 107367760B CN 201710525466 A CN201710525466 A CN 201710525466A CN 107367760 B CN107367760 B CN 107367760B
- Authority
- CN
- China
- Prior art keywords
- estimated
- function
- green
- source wavelet
- domain
- 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.)
- Expired - Fee Related
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 107
- 238000000034 method Methods 0.000 title claims abstract description 98
- 238000001228 spectrum Methods 0.000 title claims abstract description 44
- 230000008859 change Effects 0.000 claims abstract description 33
- 230000001133 acceleration Effects 0.000 claims abstract description 13
- 239000011159 matrix material Substances 0.000 claims description 63
- 238000013507 mapping Methods 0.000 claims description 8
- 238000004088 simulation Methods 0.000 claims description 8
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 8
- 230000009466 transformation Effects 0.000 claims description 5
- 230000004044 response Effects 0.000 abstract description 18
- 230000003595 spectral effect Effects 0.000 abstract 1
- 230000006870 function Effects 0.000 description 58
- 238000010586 diagram Methods 0.000 description 28
- 238000012545 processing Methods 0.000 description 10
- 241001323321 Pluto Species 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 6
- 238000012360 testing method Methods 0.000 description 5
- 230000005284 excitation Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 210000001367 artery Anatomy 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 239000010410 layer Substances 0.000 description 3
- 239000003208 petroleum Substances 0.000 description 3
- 238000004080 punching Methods 0.000 description 3
- 210000003462 vein Anatomy 0.000 description 3
- 241001269238 Data Species 0.000 description 2
- 230000001186 cumulative effect Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 239000013535 sea water Substances 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 230000017105 transposition Effects 0.000 description 2
- 102000002274 Matrix Metalloproteinases Human genes 0.000 description 1
- 108010000684 Matrix Metalloproteinases Proteins 0.000 description 1
- 230000001154 acute effect Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000006854 communication Effects 0.000 description 1
- 238000011109 contamination Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000013079 data visualisation Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 239000003989 dielectric material Substances 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 239000011229 interlayer Substances 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 229940050561 matrix product Drugs 0.000 description 1
- 238000003825 pressing Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000002945 steepest descent method Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000013518 transcription Methods 0.000 description 1
- 230000035897 transcription Effects 0.000 description 1
- JLYXXMFPNIAWKQ-UHFFFAOYSA-N γ Benzene hexachloride Chemical compound ClC1C(Cl)C(Cl)C(Cl)C(Cl)C1Cl JLYXXMFPNIAWKQ-UHFFFAOYSA-N 0.000 description 1
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. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/56—De-ghosting; Reverberation compensation
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Complex Calculations (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a kind of based on the surface-related multiple and higher-order spectra method and system that accelerate linear Bregman algorithm, feedback model based on description primary wave and surface-related multiple relationship, this method carries out sparse inversion first with the linear Bregman algorithm of acceleration, to obtain the impulse response unrelated with surface, then the source wavelet with big gun collection spatial position change is updated using a step linear spectral estimation formulas, and the impulse response and source wavelet unrelated with surface are updated by cycle alternation, final estimation obtains primary wave and surface-related multiple.Compared to traditional multi-parameter, the surface-related multiple estimation method based on Chang Zibo hypothesis and complicated algorithm, implement estimation method and system of the invention, the source wavelet with big gun collection spatial position change can be estimated, and high-efficient, the time-consuming precision for less, estimating result of estimation method and system is high.
Description
Technical field
The present invention relates to seismic data processing aspects in field of seismic exploration more particularly to seismic exploration technique, more specifically
Ground is said, is related to a kind of based on the surface-related multiple and higher-order spectra method and system that accelerate linear Bregman algorithm.
Background technique
Multiple reflection (multiple wave) relevant issues are most one of distinct issues generally existing in reflection shooting.
Since seawater surface is a strong reflection interface, surface-related multiple problem is particularly acute in marine seismic prospectiong.Now
Offshore oil natural gas-combustible ice resource exploration and exploitation have become the task of top priority of domestic and international major petroleum-based energy company.I
State's south sea petroleum natural gas-combustible ice is resourceful, is the hot spot and focus of current seafari.It is produced in current petroleum industry
In, a basic assumption of conventional seismic data processing technology is: the reflectance data body of input is only (primary by primary reflection
Wave) composition.Multiple wave can be mistakenly considered a part for being primary wave or being blended in primary wave.Multiple wave is to one
The processing and imaging of subwave cause to interfere strongly, influence resource exploration, so multiple wave estimation technique has fabulous practical work
Industry application prospect.
The multiple wave estimation or drawing method that industry uses at present are broadly divided into following a few classes: one kind using multiple wave with
The feature difference multiple suppression of primary wave, referred to as filter method；One kind predicts multiple wave from seismic data first, then will
It is adaptively subtracted from initial data, referred to as prediction subtractive method；There are also a kind of inverse-scattering series method multiple suppressions.
Filter method is premised on multiple wave and primary wave have apparent feature difference in particular transform domain it is assumed that passing through
Various mathematical algorithms are designed to obtain this species diversity and Multiple attenuation.The main feature difference that filter method utilizes includes multiple wave
Periodicity and separability.Class method is filtered when being applied to the multiple attenuation of some complicated structures, multiple wave can be encountered
Feature difference very little between primary wave does not either have the case where difference, at this moment filters the applicable precondition of class method and obtains
Less than satisfaction, multiple wave pressing result is unobvious.
Back scattering series predicts multiple wave, and there is also some problems: computationally intensive, at high cost；It is difficult to predict remote offset distances
Multiple wave；The predictive ability of multiple wave also has certain limitation in complex dielectrics.Among these, calculate it is at high cost so that this method very
Hardly possible is applied in the processing of real data.
Berkhout proposes the feedback model theory for describing complicated multiple wave in nineteen eighty-two, and it is multiple to have established feedback iteration
The mathematical physics basis of wave drawing method, but known source wavelet is needed when Multiple attenuation.In order to illustrate this problem, Yi Jibian
In subsequent explanation, we introduce the feedback model formula closely related with the present invention first, namely
Wherein, D is the expression matrix form of the time-domain of collected input data,It is D after Fourier's direct transform
One slice of the matrix expression of frequency domain, S are the matrix expression of the time-domain of source wavelet to be estimated,For S warp
One of the matrix expression of frequency domain slice after Fourier's direct transform, G are Green's function (namely the underground arteries and veins without multiple wave
Punching response) time-domain matrix expression,For a slice of G matrix expression of frequency domain after Fourier's direct transform
We indicate the data or variable of frequency domain, r with ^ symbolsurfFor the reflection coefficient of air and the contact surface of water, it is assumed that rsurf
=-1, whereinWithFrequency domain product (namely time domain convolution) represent estimation primary wave frequency slice, namelyWhereinFor a slice of the frequency domain of primary wave.Primary wave involved in the present invention refers to except surface is multiple
All waves except wave, namely primary wave here contain interbed multiple (it is unrelated with surface).Table
Show the slice namely data of the frequency domain of the surface-related multiple of estimationIt is propagated downwards as secondary focus, with underground arteries and veins
Punching response Green's functionFrequency domain product (time domain convolution) predict the slice of surface-related multiple frequency domain.By equation
1 formula can be seen that the precision of source wavelet of the accurate surface-related multiple estimation dependent on estimation.
Verschuur etc. proposed surface-related multiple compact technique SRME in 1992, was based on one by introducing
What secondary wave energy minimum was assumed adaptively subtracts each other, and source wavelet is estimated from real data.In order to avoid due to solving focus
Wave and the nonlinear problem generated, Berkhout and Verschuur were in proposition iteration SRME method in 1997.Wave energy is most
This small hypotheses condition is not that can meet in all cases.In addition, in practical applications, SRME there is also
Adaptive prediction subtracts each other the problem of damage useful signal.
There are problems that damage useful signal, van Groenestijn and Verschuur are subtracted each other in prediction to improve SRME
In proposition sparse inversion primary wave estimation techniques in 2009 namely EPSI.The problem of in order to improve EPSI, Lin and
Herrmann proposed steady sparse inversion primary wave estimation technique namely REPSI in 2013.EPSI and REPSI be with
The more relevant prior art of the present invention, their technical solution and there are the problem of, disadvantage will be explained in detail below
It states.
Further, since the source wavelet that every big gun is excited is not quite similar, and current institute when practical field earthquake data acquisition
Some multiple wave estimation methods assume that all source wavelets are identical, influence multiple wave estimated accuracy, and the present invention can estimate
With the source wavelet of big gun collection spatial position change.
The computing resource (memory, processor) that each iteration expends when due to multiple wave estimation is more, it is longer to calculate the time,
Other than improving multiple wave, the estimated accuracy of source wavelet, the efficiency of improvement method, simplifies and calculates the rate of convergence of boosting algorithm
The development trend of the complexity of method, the adjustment parameter of reduction method and multiple wave estimation technique.
The prior art related to the present invention is described below.
The technical solution of the prior art one
Van Groenestijn and Verschuur proposed sparse inversion one in 2009 in the theoretical basis of SRME
Subwave estimation method EPSI, this method are driving with fitting data residual error, are handed over by steepest descent method during iterative inversion
For update underground impulse response G and two unknown parameters of source wavelet S, and then reconstruct primary waveAnd its corresponding table
Face related multipleCompared with SRME, EPSI avoids the step of adaptively subtracting each other, and preferably protects primary
Wave useful signal.
The EPSI method of van Groenestijn and Verschuur are when carrying out source wavelet estimation, it is assumed that all shakes
Source wavelet is identical, does not provide the estimation scheme (namely normal higher-order spectra) with big gun collection spatial position change source wavelet.This
Invention can be estimated with big gun collection spatial position change source wavelet.
In order to guarantee underground impulse response G in the sparsity of time domain, EPSI during updating impulse response G, according to
Time is ascending, amplitude is by by force to weak sequential update pulse.Also in order to promoting the time domain of underground impulse response G sparse
Property, EPSI method devise numerous artificial adjustment parameters, for example, EPSI it needs to be determined that: a time window, for determining Green's letter
Several more new ranges (namely being updated in time window, the zero setting outside time window)；It is reflected in Green's function when each iteration
The number of lineups, etc..As Lin and Herrmann in paper in 2013 pointed by: even if being accurately known that underground is anti-
The number (we can not can know that in practice) in firing area face, EPSI method can not also make full use of this prior information.EPSI
Underground impulse response update, be easy by pulse sliced window and pick up the parameters such as quantity and influenced, and then influence multiple wave
Estimated accuracy.Intrinsically, the degree of rarefication of Green's function is controlled by artificially in EPSI method, rather than by mathematics
Optimization algorithm controls.This disadvantage limits application of the EPSI method in actual industrial.
The prior art two related to the present invention
Based on EPSI identical theoretical (namely 1 formula of feedback model), Lin and Herrmann are in 2013 by the underground EPSI arteries and veins
The sparse constraint condition of punching response is revised as a norm constraint, and using spectrum mapping one norm algorithm of gradientTo solve
To sparse underground impulse response Green's function, the stability of original EPSI method is improved to a certain extent.Lin and
Herrmann is named as Robust EPSI (REPSI) in proposed method in 2013.What Lin and Herrmann was proposed
REPSI method, it is the same with EPSI, and alternately update underground impulse response G and two unknown parameters of source wavelet S, Jin Erchong
Structure primary waveAnd its corresponding surface-related multipleThe maximum difference of EPSI and REPSI exists
In: the degree of rarefication majority of Green's function is controlled by artificial adjustment parameter in EPSI, and REPSI is then usedMathematics
Optimization algorithm solves to obtain sparse underground impulse response Green's function.
The REPSI method of Lin and Herrmann also assumes that all source wavelet phases when carrying out source wavelet estimation
Together, the estimation scheme (namely normal higher-order spectra) with big gun collection spatial position change source wavelet is not provided yet.The present invention can
Estimation is with big gun collection spatial position change source wavelet.In addition, REPSI is changed when carrying out normal higher-order spectra using LSQR algorithm
In generation, solves and obtains source wavelet, considerably increases the calculation amount of method, reduces efficiency.And a step of the invention can be estimated to obtain
Source wavelet.
Another significant drawbacks of the REPSI method of Lin and Herrmann are: it relies on an extremely complex spectrum mapping
One norm algorithm of gradientTo solve to obtain sparse underground impulse response Green's function.It is an open source
MATLAB program can be downloaded by access https: //github.com/mpf/spgl1.The complexity master of algorithm
It is embodied in: numerous adjustment parameters, numerous subfunctions, upper Kilo Lines of Code.Schlumberger (Schlumberger) oilfield services are public
Department once attempts reorganizationIt is applied in actual industrial, but for fear ofThe complexity of algorithm, it is final and failed.The complexity of algorithm also reduces the efficiency of REPSI method simultaneously, hinders the reality of REPSI method in the industry and answer
With.
Summary of the invention
The technical problem to be solved in the present invention is that for technology existing for above-mentioned EPSI in the prior art and REPSI
Defect provides a kind of based on the surface-related multiple and higher-order spectra method and system that accelerate linear Bregman algorithm.
Wherein one side according to the present invention, the present invention are to solve its technical problem, are provided a kind of linear based on acceleration
The surface-related multiple and higher-order spectra method of Bregman algorithm include following step:
S1, Green's function to be estimated and source wavelet to be estimated are alternately updated using following step S11 and S12
Until reach preset times, wherein carrying out the source wavelet to be estimated of step S11 for the first time has preset initial value:
S11, according to the input data of source wavelet to be estimated and acquisition, using accelerating linear Bregman algorithm more
New Green's function to be estimated；
S12, according to updated Green's function to be estimated and the input data, updated using following formula wait estimate
The source wavelet of meter；
When source wavelet to be estimated is not with big gun collection spatial position change:
When source wavelet to be estimated is with big gun collection spatial position change:
S2, according to after final updated as a result, obtaining the estimation knot of at least one of source wavelet and surface-related multiple
Fruit；
Wherein, subscript est is indicated to be estimated,The element on the diagonal line of source wavelet matrix slice is represented, j is indicated
It is numbered with big gun collection spatial position or big gun collection, vec () indicates a matrix becoming column vector, subscriptRepresentative takes a frequency
It is sliced the jth column of matrix, subscriptAll column for taking a frequency slice matrix are represented,For the square of the frequency domain of input data
The slice of battle array expression formula,For the slice of the matrix expression of the frequency domain of source wavelet, rsurfFor the contact surface of air and water
Reflection coefficient,For the slice of the matrix expression of the frequency domain of Green's function.
Further, of the invention based on the surface-related multiple and higher-order spectra method that accelerate linear Bregman algorithm
In, it is to be estimated using accelerating linear Bregman algorithm to update according to source wavelet to be estimated and the input data of acquisition
Green's function specifically includes:
According to formula AndIt is iterated, successively finds out y1、x1、y2、x2、y3…、Wherein,It is held for this
The vector expression-form of the time-domain of the update result of Green's function to be estimated after row step S11, lmaxStep is executed for this
The maximum number of iterations of the linear Bregman algorithm of acceleration when rapid S11；
In formula, x0、With y0It is equal to the vector expression shape of the time-domain of the preset initial value of Green's function to be estimated
Formula, observation data vector b are the vector expression-form of the time-domain of the input data,What t was represented is step
It is long, expression formula are as follows:
Mapping function ∏σ(Axl- b) is defined as:
Soft-threshold function shrink (x, μ) is defined as:
Shrink (x, μ)=max (| x |-μ, 0) sign (x),
Wherein, σ represents the noise level factor, and max () function representation is maximized, | | function representation takes absolute value, sign
(x) symbol value function is indicated；
Matrix operator A indicates to carry out the operational set of operations described below: by Green's function to be estimated by the vector of time-domain
Sheet form, which carries out Fourier transformation, becomes the vector sheet form of frequency domain, by Green's function to be estimated by the vector table of frequency domain
Form becomes the matrix sheet form of frequency domain, according to formulaIt is corresponding to calculate separately input data
Each slice of the expression matrix form of the frequency domain of analogue dataAccording to each sliceCarry out inversefouriertransform
Obtain the vector expression-form of the time-domain of the corresponding analogue data of input data.
Further, of the invention based on the surface-related multiple and higher-order spectra method that accelerate linear Bregman algorithm
In, in step S2,
The estimated result of source wavelet is obtained according to the result for updating source wavelet to be estimated for the last time；
The estimated result of surface-related multiple is according to m=rsurfG*d is obtained, and m is the vector table of the time-domain of surface-related multiple
Up to form, g is the vector expression-form of the time-domain of Green's function, and d is the vector expression-form of the time-domain of input data.
Further, of the invention based on the surface-related multiple and higher-order spectra method that accelerate linear Bregman algorithm
In, if Green's function to be estimated has reached the maximum update times, source wavelet to be estimated is not up to maximum update time
Number, then when carrying out step S11 again, lmaxThe maximum update times for being updated to source wavelet to be estimated, which subtract, to be updated
Number after difference, and before every time execute step S11 when lmaxIt is equal；If source wavelet to be estimated has reached most
Big update times and Green's function to be estimated is not up to maximum update times, then step S12 of subsequent progress, and it is each
Execute l when step S11maxIt is equal.
The present invention is to solve its technical problem, and it is multiple to additionally provide a kind of surface based on the linear Bregman algorithm of acceleration
Wave and higher-order spectra system include:
Alternately updating unit, for using following Green's function updating unit and source wavelet updating unit alternating more
New Green's function to be estimated and source wavelet to be estimated are until reach preset times, wherein calling Green's function to update for the first time
The source wavelet to be estimated of unit has preset initial value:
Green's function updating unit, it is linear using accelerating according to source wavelet to be estimated and the input data of acquisition
Bregman algorithm updates Green's function to be estimated；
Source wavelet updating unit, according to updated Green's function to be estimated and the input data, under
It states formula and updates source wavelet to be estimated；
When source wavelet to be estimated is not with big gun collection spatial position change:
When source wavelet to be estimated is with big gun collection spatial position change:
Estimated result acquiring unit, for according to it is after final updated as a result, obtain source wavelet to be estimated and to
The estimated result of at least one of the surface-related multiple of estimation；
Wherein, subscript est is indicated to be estimated,The element on the diagonal line of source wavelet matrix slice is represented, j is indicated
It is numbered with big gun collection spatial position or big gun collection, vec () indicates a matrix becoming column vector, subscriptRepresentative takes a frequency
It is sliced the jth column of matrix, subscriptAll column for taking a frequency slice matrix are represented,For the square of the frequency domain of input data
The slice of battle array expression formula,For the slice of the matrix expression of the frequency domain of source wavelet, rsurfFor the contact surface of air and water
Reflection coefficient,For the slice of the matrix expression of the frequency domain of Green's function.
Further, of the invention based on the surface-related multiple and higher-order spectra system that accelerate linear Bregman algorithm
In, it is to be estimated using accelerating linear Bregman algorithm to update according to source wavelet to be estimated and the input data of acquisition
Green's function specifically includes:
According to formulaAndIt is iterated, successively finds out y1、x1、y2、x2、y3…、Wherein,It is held for this
The vector expression-form of the time-domain of the update result of Green's function to be estimated after row step S11, lmaxStep is executed for this
The maximum number of iterations of the linear Bregman algorithm of acceleration when rapid S11；
In formula, x0、With y0It is equal to the vector expression of the time-domain of the preset initial value of Green's function to be estimated
Form, observation data vector b are the vector expression-form of the time-domain of the input data,tlRepresent be
Step-length, expression formula are as follows:
Mapping function ∏σ(Axl- b) is defined as:
Soft-threshold function shrink (x, μ) is defined as:
Shrink (x, μ)=max (| x |-μ, 0) sign (x),
Wherein, σ represents the noise level factor, and max () function representation is maximized, | | function representation takes absolute value, sign
(x) symbol value function is indicated；
Matrix operator A indicates to carry out the operational set of operations described below: by Green's function to be estimated by the vector of time-domain
Sheet form, which carries out Fourier transformation, becomes the vector sheet form of frequency domain, by Green's function to be estimated by the vector table of frequency domain
Form becomes the matrix sheet form of frequency domain, according to formulaIt is corresponding to calculate separately input data
Each slice of the expression matrix form of the frequency domain of analogue dataAccording to each sliceCarry out anti-Fourier's change
Get the vector expression-form of the time-domain of the corresponding analogue data of input data in return.
Further, of the invention based on the surface-related multiple and higher-order spectra system that accelerate linear Bregman algorithm
In, in estimated result acquiring unit,
The estimated result of source wavelet is obtained according to the result for updating source wavelet to be estimated for the last time；
The estimated result of surface-related multiple is according to m=rsurfG*d is obtained, and m is the vector table of the time-domain of surface-related multiple
Up to form, g is the vector expression-form of the time-domain of Green's function, and d is the vector expression-form of the time-domain of input data.
Further, of the invention based on the surface-related multiple and higher-order spectra system that accelerate linear Bregman algorithm
In, if Green's function to be estimated has reached the maximum update times, source wavelet to be estimated is not up to maximum update time
Number, then when carrying out step S11 again, lmaxThe maximum update times for being updated to source wavelet to be estimated, which subtract, to be updated
Number after difference, and before every time execute step S11 when lmaxIt is equal；If source wavelet to be estimated has reached most
Big update times and Green's function to be estimated is not up to maximum update times, then step S12 of subsequent progress, and it is each
Execute l when step S11maxIt is equal.
Implement to have below in the estimation method and system of the invention based on the linear Bregman algorithm of acceleration beneficial to effect
Fruit: estimation method and system of the invention can estimate the source wavelet with big gun collection spatial position change, and estimation method and be
High-efficient, the time-consuming precision for less, estimating result of system is high.
Detailed description of the invention
Present invention will be further explained below with reference to the attached drawings and examples, in attached drawing:
Fig. 1 is that the linear Bregman algorithm of acceleration of the invention solves Ax=b to obtain mono- preferred embodiment of sparse solution x
Flow chart；
Fig. 2 is the flow chart of the estimation source wavelet and one preferred embodiment of surface-related multiple of estimation method of the invention；
Fig. 3 is the schematic diagram of error between true solution, the solution that some algorithm is found out and the solution and true solution that find out；
Fig. 4 is comparison diagram of the data residual error that finds out of some algorithm with the change curve of the number of iterations；
The schematic diagram data comprising multiple wave of Fig. 5 (a) input；
Fig. 5 (b) is the primary wave schematic diagram data that the present invention estimates；
Fig. 5 (c) is the surface-related multiple schematic diagram data that the present invention estimates；
Fig. 6 (a) is that zero big gun of data in Fig. 5 (a) examines diagrammatic cross-section；
Fig. 6 (b) is the surface-related multiple schematic diagram data that the present invention estimates；
Fig. 6 (c) is the primary wave schematic diagram data of one embodiment of the invention estimation；
The primary wave schematic diagram data that the REPSI method that Fig. 6 (d) is Lin and Herrmann is estimated；
Fig. 7 (a) is the schematic diagram of the source wavelet of estimation；
Fig. 7 (b) is the frequency domain amplitude spectrum schematic diagram of the source wavelet with the variation of big gun collection of one embodiment of the invention estimation；
Fig. 8 is the REPSI method of Lin and Herrmann, original linear Bregman algorithm, accelerates linear Bregman algorithm
Inversion error with the number of iterations change curve
Fig. 9 (a) is the open release of world SMAART JV tissue for testing the Pluto 1.5 of multiple wave drawing method
Schematic diagram data；
Fig. 9 (b) is the primary wave schematic diagram data that the present invention estimates；
Fig. 9 (c) is the surface-related multiple schematic diagram data that the present invention estimates；
Figure 10 (a) is that zero big gun of data in Fig. 9 (a) examines diagrammatic cross-section；
Figure 10 (b) is the primary wave schematic diagram data that the present invention estimates；
The surface-related multiple schematic diagram data that Figure 10 (c) present invention estimates；
Figure 11 (a) is the source wavelet schematic diagram changed with big gun collection estimated from 1.5 data of Pluto；
Figure 11 (b) is that the frequency domain amplitude spectrum for the source wavelet with the variation of big gun collection estimated from 1.5 data of Pluto shows
It is intended to；
Figure 12 (a) is the data diagrammatic cross-section of the intermediate big gun of data in Fig. 9 (a)；
Figure 12 (b) is the primary wave schematic diagram data that the present invention estimates；
Figure 12 (c) is the single order surface-related multiple schematic diagram data that the present invention estimates；
Figure 12 (d) is the second-order surface multiple wave schematic diagram data that the present invention estimates；
Figure 12 (e) is the three rank surface-related multiple schematic diagram datas that the present invention estimates；
Figure 12 (f) is the quadravalence surface-related multiple schematic diagram data that the present invention estimates；
Figure 12 (g) is the five rank surface-related multiple schematic diagram datas that the present invention estimates；
Figure 12 (h) is the schematic diagram of Figure 12 (b) cumulative summation of data into Figure 12 (g)；
Figure 12 (i) is the schematic diagram of the difference of Figure 12 (a) and Figure 12 (h)；
Figure 13 (a) is big gun number in 1.5 data of Pluto between the schematic diagram of the zero big gun inspection office section of 701-880；
Figure 13 (b) is the primary wave schematic diagram data that the present invention estimates；
Figure 13 (c) is the surface-related multiple schematic diagram data that the present invention estimates.
Noun explanation
Primary wave (Primaries): refer to that the seismic wave of excitation only occurs one in subsurface interface during earth-layer propagation
The wave detector that secondary reflection is embedded in earth's surface is received, referred to as primary reflection, abbreviation primary wave.Involved in the present invention
Primary wave refers to all waves in addition to surface-related multiple, namely primary wave here contains interbed multiple (itself and surface
It is unrelated).
Multiple wave (Multiples): refer to that the seismic wave of excitation is (empty in subsurface interface or Free Surface in communication process
The contact surface of gas and water) occur just to be embedded in the wave detector of earth's surface after multiple reflections and is received, referred to as multiple reflection, abbreviation
Multiple wave.
Surface-related multiple (Surface-related multiples) or surface-related multiple (Surface
Multiples): referring to multiple reflection relevant to Free Surface (contact surface of air and water), summary in the present invention is
Surface-related multiple.
Interbed multiple (Internal multiples, IM): refer between subsurface reflective boundary and produced (between stratum)
Raw multiple reflection.
Source wavelet: the wave that epicenter excitation generates is referred to.
Feedback model: feedback model.
Surface-Related Multiple Elimination (SRME): surface-related multiple compacting.
Estimation of Primaries by Sparse Inversion (EPSI): sparse inversion primary wave estimation.
Robust Estimation of Primaries by Sparse Inversion (EPSI): steady sparse anti-
Drill primary wave estimation.
: Spectral-ProjectedSpectrum mapping one norm algorithm of gradient.
LSQR:Least-squares QR, least square QR- factorization algorithm.
Linearized Bregman (LB): linear Bregman algorithm.
Accelerated Linearized Bregman (ALB): accelerate linear Bregman algorithm.
MATLAB: be MathWorks company, the U.S. produce business mathematics software, for algorithm development, data visualization,
The advanced techniques computational language and interactive environment that data analysis and numerical value calculate.MATLAB is matrix&laboratory
Two contaminations mean matrix factory (matrix labotstory).
Signal-to-Noise Ratios (SNR): signal-to-noise ratio.
Specific embodiment
For a clearer understanding of the technical characteristics, objects and effects of the present invention, now control attached drawing is described in detail
A specific embodiment of the invention.In order to facilitate the understanding of the present invention, now above-mentioned and following each expression formula letters are made such as
Lower explanation, alphabetical s, g, d, p, m of small letter are the vector expression-form for indicating time-domain, and S, G, D, P, M are the matrix of time-domain
Expression-form.For the slice of the expression matrix form of the frequency domain of S, G, D, P, M, subscript est is indicated
It is to be estimated,For the vector expression-form of the frequency domain of s, g, d, p, m, S, G, D, P, M are three-dimensional matrice,
And its slice is then two-dimensional matrix, each three-dimensional matrice can be split as multiple two-dimensional matrixes.Any one two-dimensional matrix Z=
[Z1Z2Z3…Zn] vector expression-form can be expressed as:
Wherein Z1To ZnThe column for indicating two-dimensional matrix Z, conversely, expression matrix can be obtained by vector expression-form
Form.
The main formulas that primary wave, surface-related multiple and source wavelet estimation utilize is the mathematical physics relationship of feedback model
Formula namely formula (1).The source wavelet s of excitation is passed to underground, carries out time domain convolution (star with underground impulse response Green's function g
Number * indicates convolution), generate primary wave p, mathematical relationship expression formula are as follows:
P=g*s (2)
The primary wave p of generation encounters seawater and reflects with air contact surfaces, propagates downwards, and then pass again as focus
Enter underground, carries out time domain convolution with underground impulse response g, generate single order multiple wave m1
m1=g*rsurfP=rsurfg*g*s (3)
The single order multiple wave m of generation1It encounters Free Surface to reflect, and then is passed to underground again as focus, with ground
Lower impulse response g carries out time domain convolution, generates second-order multiples m2
m2=g*rsurfm1=(rsurf)2g*g*g*s (4)
The surface-related multiple of higher order time generates formula and so on.
Observe the summation that obtained total data d is primary wave p and all order surface-related multiples, mathematical expression
Formula are as follows:
And then have
Above-mentioned feedback model mathematical physics relational expression can summarize are as follows:
Expression formula (1) is frequency domain expression matrix form of the expression formula 7 under two dimension or three-dimensional situation, for formula 7
In parameter be revised as the expression matrix form of its corresponding frequency domain and also set up, for convenient for comparison, by its transcription in this:
The main object of the present invention is: according to expression formula 7-8, known variables: source wavelet s and Green's function g are solved, into
And the primary wave p=g*s and surface-related multiple m=r estimatedsurfg*d。
In the case of the present invention illustrates the source wavelet of given estimation first, inverting Sparse Pulse responds the side of Green's function
Method.
In the source wavelet s of given estimationestIn the case where, according to expression formula (7)-(8), we write out inverting gestWhen institute
Formula, namely:
HereIt is the representative forward operator of simplification, core is in calculation formula 8
Specifically, matrix operatorIndicate the operational set of progress operations described below: by Green's function to be estimated by the vector of time-domain
Sheet form g, which carries out Fourier transformation, becomes the vector sheet form of frequency domainBy Green's function to be estimated by the vector of frequency domain
Sheet formBecome the matrix sheet form of frequency domainAccording to formulaCalculate separately input data pair
Each slice of the expression matrix form of the frequency domain for the analogue data answeredAccording to the corresponding mould of calculated input data
Each slice of the expression matrix form of the frequency domain of quasi- dataIt is corresponding that progress inversefouriertransform obtains input data
The vector expression-form of the time-domain of analogue data
How to solve to obtain sparse Green's function g for the ease of the present invention is further explained, 9 formulas are rewritten as to following mark
Quasi- form:
Ax=b (10)
Wherein,X=g, b=d.
In order to seek the sparse solution of Green's function g, the present invention solves 10 formulas using the linear Bregman algorithm of acceleration.Accelerate
Linear Bregman algorithm is intended to solve:
Expression solves x and expression formula is made to obtain minimum value, and subject to is indicated in limited conditions.
Wherein, μ >=0 is degree of rarefication controlling elements, | | x | |1With | | x | |2Respectively represent vectorNorm andNorm, σ
Represent the noise level factor in data.Following algorithm 1 gives the detailed iterative step for accelerating linear Bregman algorithm.
It is that the linear Bregman algorithm of acceleration of the invention solves Ax=b to obtain the process of sparse solution x with reference to Fig. 1, Fig. 1
The flow chart of figure namely algorithm 1, algorithm 1 specifically include following step:
A) parameters described below: forward operator A is obtained, data vector b, initial value vector x are observed0, threshold value μ, greatest iteration time
Number lmax。
B) it initializes: iteration variable l=0, by x0It is assigned to y0
C) repeat following each d to g step lmaxIt is secondary
D) it calculates
E) it calculates
F) it calculates
G) l is updated to l+1
H) sparse solution is obtainedSparse solution x is that this carries out accelerating linear Bregman algorithm
Update gestResult afterwards.
In above-mentioned algorithm 1, the t of step dlWhat is represented is step-length, expression formula are as follows:
Wherein, subscript H represents complex conjugate transposition.Mapping function ∏σ(Axl- b) be for consider data present in noise,
Is defined as:
And soft-threshold function shrink (x, μ) is defined as:
Shrink (x, μ)=max (| x |-μ, 0) sign (x) (14)
Wherein, max () function representation is maximized, | | function representation takes absolute value, and sign (x) indicates symbol value
Function.
In above-mentioned algorithm 1, αlDefinition expression formula are as follows:
This method passes through weighting system αlCarry out the rate of convergence of accelerating algorithm.
The present invention is following will to be illustrated to estimate source wavelet in the case where inverting obtains Sparse Pulse response Green's function.
When inverting obtains gestAfterwards, according to formula 8, we estimate to solve the problems, such as when source wavelet are as follows:
At this timeFor unknown number to be solved.We carry out source wavelet with a simple derivation to illustrate the present invention
The scheme of estimation.Estimate source wavelet, i.e., to estimate the Fourier spectrum of source wavelet frequency domain.It is assumed that two vector of complex values a
Meet lower relation of plane with b:
Ac=b (17)
Wherein c is source wavelet Fourier transform spectrum complex coefficient to be asked.And then we have:
aHAc=aHb (18)
And then obtain,
The present invention is to estimate the source wavelet with big gun collection spatial position change according to derivation above.According to 16 He of formula
Formula 19 has when the source wavelet of estimation is not with big gun collection spatial position change (namely normal higher-order spectra):
When the source wavelet of estimation is with big gun collection spatial position change (namely becoming higher-order spectra), have:
Wherein, S(j,j)The element on source wavelet diagonal of a matrix is represented, j indicates to number with big gun collection spatial position or big gun collection
One matrix is become column vector by variation, vec () expression, and H represents complex conjugate transposition, subscriptRepresentative takes a frequency slice
The jth of matrix arranges, subscriptRepresent all column for taking a frequency slice matrix.When carrying out source wavelet estimation, the present invention
Any assumed condition is not needed, namely does not need to do any hypothesis to the phase of wavelet.In addition, not considering repeatedly compared to tradition
The higher-order spectra problem of wave, the multiple wave item when present invention carries out higher-order spectra help to overcome existing for traditional higher-order spectra
(amplitude and phase) nonuniqueness problem.
Fig. 2 is the flow chart of a preferred embodiment of estimation method of the invention, which includes the following steps:
A) following data: data d, threshold value μ, g are obtainedestGreatest iteration update times kmax, linear Bregman is accelerated to calculate
Method inner iterative number lmax, sestGreatest iteration update times mmax。
B) it initializes: iteration variable k=0, m=0, b=d, gestAnd sestIt is initialized as 0.
C) judge k < kmaxIt is whether true, if so, step c-i is then repeated, it is no to then follow the steps j.
D) according to x0=gest、sestAnd b=d, call algorithm 1 to update to obtain gest。
E) k is updated to k+lmax。
F) judge m < mmaxIt is whether true, if so, then follow the steps g-i.
G) according to b and updated gest, formula 20 or 21 formulas is called to update to obtain sest.It is calculated by formula 20 or 21
Each diagonal matrixIn diagonal line on each element, according to calculated allUpdate sest。
H) judge m=mmaxWhether -1 is true, by k if setting upmax- k is assigned to lmax。
I) m+1 is assigned to m.
J) according to after final updated as a result, obtaining the estimated result of source wavelet and surface-related multiple.
There are two unknown numbers in the above process: the slice s of source waveletestWith the slice g of Green's functionest, anti-solving
Traditional blind deconvolution method is similar in problem process, i.e., fixed sest, carry out inverting gest；Then g is fixedest, to estimate sest。
In inverting gestWhen, previous inversion result has been used in each iteration more new capital.In the sparse g of invertingestXie Shi generally requires to change
Generation several times, because according to sparse Renew theory, generally requiring iteration several times could restrain to obtain an optimal sparse solution.It is obtaining
One sparse solution gestAfterwards, then source wavelet estimation is carried out.Estimate obtained sestIt is immediately used for inverting and updates gest.Entire iteration
Major cycle terminates meeting some stopping criterion backed off after random.Alternately update gestWith sestPreset times be decided by kmaxAnd
mmax, work as gestReach k after maximum update timesmaxAnd sestWhen not reaching maximum update times, g is being executedestLast time
After update, a s is executedestUpdate complete gestWith sestAlternating update, wherein carrying out accelerating every time linear
The l of Bregman algorithmmaxIt is equal；Work as sestReach k after maximum update timesmaxAnd gestWhen not reaching maximum update times,
Execute sestLast time update after, only carry out last time and accelerate linear Bregman algorithm, last time executes acceleration line
The l of property Bregman algorithmmaxIt is assigned gestMaximum update times and the difference of number that has been updated, and before every
It is secondary to execute the l for accelerating linear Bregman algorithmmaxIt is equal.Wherein carry out accelerating the l of linear Bregman algorithm every timemaxIt is equal most
G is once updated afterwardsestThat is the estimated result of g, last time update sestThe as estimated result of s, it is following to use gestIndicate g's
Estimated result, sestIndicate the estimated result of s.
Obtaining gestAnd sestAfterwards, it calculates and obtains primary wave p=gest*sest。
By gestAnd input dataObtain surface-related multiple m=rsurfgest*d。
According to secondary source principles, Green's function g is being obtainedestAfter the primary wave p of acquisition, by feedback iteration model meter
It calculates and obtains the slice of the surface-related multiple time domain of different orders, for example, single order surface-related multiple m1Calculation formula are as follows:
m1=rsurfgest* p=rsurfgest*gest*uest (22)
The slice m of the time domain of second-order surface multiple wave2Calculation formula are as follows:
m2=rsurfgest*m1=(rsurf)2gest*gest*gest*sest (23)
The calculation formula of surface-related multiple and so on of other orders.
Computationally intensive part when due to multiple wave estimation, which is consumed, is doing matrix-matrix product in each frequency slice, namely
8 formulas are calculated, so in order to be further reduced the operation time of multiple wave estimation, the present invention limits characteristic according to the band of seismic data
(namely the frequency content of seismic data is concentrated mainly in the frequency range of a certain finite length), the present invention is in multiple wave, shake
Merely with the dominant frequency band that seismic data confidence level and signal-to-noise ratio are high when the higher-order spectra of source, the calculating consumption of method is greatly reduced
When, this is different from conventional method and all (low including confidence level and signal-to-noise ratio) frequency contents is utilized.
After finding out primary wave or surface-related multiple, so that it may obtain underground structure figure respectively using them.
Here give following specific example to illustrate the feasibility of technical solution of the present invention, validity and excellent
More property.
(1) sparse recovery numerical example, to illustrate that the present invention accelerates the superiority of linear Bregman algorithm.
(2) In A Salt-dome Model data example.The data are small by the DELPHI multiple wave research of TU Delft Polytechnics
Group simulates generation, is one of the signal data for being used for open test multiple wave algorithm in the world.It is once in van
Occur in papers in 2013 of the papers in 2009 of Groenestijn and Verschuur, Lin and Herrmann and for multiple
The test of wave estimation method.
(3) 1.5 data example of Pluto.The data are by Subsalt Multiples Attenuation and in the world
Reduction Team Joint Venture (SMAART JV) microstructure modeling generates and discloses release, and is used in the world
One of the signal data of open test multiple wave algorithm.
Sparse recovery numerical example
The data of true model and simulation in this numerical example are to be randomly generated.We have been randomly generated one first
Sparse vector has 30 nonzero values the length is 512 in 512 points.For 512 unknown numbers, we are randomly generated
450 observation data, namely: to 512 unknown numbers, 450 equation groups are established, this equation group system owes fixed, because
Equation number is less than unknown number number.For all methods for participating in comparison, we set maximum number of iterations as 20.Due to mould
Quasi- data not Noise, can isolate influence of the noise to all algorithms.It is rightAlgorithm we use its developer pushed away
Recommend the parameter of default.Here we evaluate the algorithm for participating in comparison with Signal to Noise Ratio (SNR) and runing time, wherein the meter of signal-to-noise ratio
Calculate formula are as follows:
Here x represents true solution,Represent the solution acquired.
From Fig. 3, (the first row representative is true sparse model namely true solution；What the second row represented is by least square
The result that QR decomposition method solves；The third line represent be byThe result that algorithm solves；Fourth line represent beThe result that algorithm solves is the same as the error really solved；The linear Bregman algorithm that fifth line represents solves to obtain
As a result；What the 6th row represented is the result that the present invention accelerates linear Bregman algorithm to solve；What the 7th row represented is this hair
The bright result for accelerating linear Bregman algorithm to solve is with the error really solved) in we can see that least square QR points
There are many random noises for the solution that resolving Algorithm acquires.The precision that algorithm acquires solution can not show a candle to linear Bregman algorithm and this
Invention accelerates linear Bregman algorithm.The error that algorithm acquires solution accelerates linear Bregman algorithm than the present invention
Error is much larger.The solution that the present invention accelerates linear Bregman algorithm to acquire has highest Signal to Noise Ratio (SNR).Fig. 3 is shown simultaneously
The present invention accelerate the runing time of linear Bregman algorithm thanThe runing time of algorithm is few very much (about to be reduced
90%).
(representative of black fine line is the data residual error of least square QR- factorization algorithm with the change of the number of iterations from Fig. 4
Change curve；Black fine dotted line represent beThe data residual error of algorithm with the number of iterations change curve；Black thick dashed line generation
Table is the data residual error of linear Bregman algorithm with the change curve of the number of iterations；That black heavy line represents is the present invention
Accelerate the data residual error of linear Bregman algorithm with the change curve of the number of iterations) show solution procedure in Data Matching miss
From the point of view of difference is with the variation tendency of the number of iterations, linear Bregman convergence speed of the algorithm is better than least square QR- factorization algorithm
WithAlgorithm.And the present invention accelerates linear Bregman convergence speed of the algorithm to be further better than linear Bregman algorithm.
Numerical examples in Fig. 3 and Fig. 4 confirm the validity for accelerating linear Bregman algorithm in the present invention and superior
Property, i.e., algorithm is simple, rate of convergence is high.
In A Salt-dome Model data example
The data that the DELPHI multiple wave research group simulation of TU Delft Polytechnics generates are shown in Fig. 5 a, the data
It is input by rate pattern and density model, using finite-difference forward modeling, time sampling interval is 4ms.When forward simulation
Source wavelet used is Ricker wavelet, and every big gun wavelet is the same namely Chang Zibo.The shallow-layer of model is the water of about 200 meters of depths
Layer.In order to which the REPSI method of same Lin and Herrmann compares, total the number of iterations is set during processing as 71.
The result of REPSI method is provided by Lin and Herrmann.The present invention frequency range used in processing is 0-70Hz (former
Entire frequency range begin as 0-125Hz).
In from Fig. 5 to Fig. 6 from the point of view of processing result of the present invention, the present invention has obtained preferable effect.In original input data
In Fig. 6 a, multiple wave is mixed in together with primary wave, interfere with each other, differentiate it is unclear.After the method for the present invention is handled, multiple wave and
Primary wave is accurately separated.The weak lineups for belonging to tomography in Fig. 6 a initial data are difficult to identify, processed by the invention
Afterwards, it is able to clearly identify in fig. 6 c, show.Due to the present invention is directed to handle multiple wave relevant to surface, and then interlayer is more
Subwave can be comprised in the primary wave obtained after processing, see Fig. 6 c.
By comparison result (Fig. 6 c) of the present invention, the REPSI methods and results (Fig. 6 d) of Lin and Herrmann, it is not difficult to find that
The primary wave that the present invention is handled is more accurate, and multiple wave is able to preferably estimate, suppress；And Lin and Herrmann
Some multiple wave lineups still can be seen in REPSI methods and results (Fig. 6 d), as shown in Fig. 6 d arrow.
Result (Fig. 7 of the REPSI method estimation for source wavelet, Lin and the Herrmann that the present invention estimates in comparison diagram 7a
(a) what solid black lines represented is the source wavelet that the present invention estimates in, and what dotted line represented is the side REPSI of Lin and Herrmann
The source wavelet of method estimation), it can be seen that the result of the two is consistent, this demonstrates the correct of the method for the present invention from side
Property.That when being generated due to this digital simulation is the Chang Zibo not changed with big gun collection, so when the present invention carries out change higher-order spectra
The result of generation is more consistent, sees Fig. 7 b.
From Fig. 8 (solid black lines represent be Lin and Herrmann REPSI method inversion error with the number of iterations
Variation；When what grey filled lines represented is using original linear Bregman algorithm (in the case of normal higher-order spectra) inversion error with repeatedly
The variation of generation number；What black dotted lines represented is using inverting (under varitron wave estimation condition) when original linear Bregman algorithm
Error with the number of iterations variation；When what dash-dotted gray line represented is acceleration linear Bregman algorithm (under varitron wave estimation condition)
Inversion error with the number of iterations variation) in iterative process from the point of view of Data Matching error change, the convergence speed of the method for the present invention
Degree will be faster than the REPSI method of Lin and Herrmann.Convergence rate when normal higher-order spectra and convergence when becoming higher-order spectra are fast
It spends almost the same.The present invention accelerates linear Bregman convergence speed of the algorithm to be better than the convergence of original linear Bregman algorithm
Speed.
Runing time in following table statistics indicate that: in the case where similarly calculating environment and computing resource, the fortune of the method for the present invention
The row time will be far smaller than the REPSI method runing time of Lin and Herrmann.Calculating time-consuming when change higher-order spectra is than normal
Calculating time-consuming when higher-order spectra is slightly more.The runing time of the method for the present invention will be far smaller than Lin's and Herrmann
The main reason for REPSI method runing time, is as follows:
(1) present invention accelerates the simplicity of linear Bregman algorithm steps, this is used with Lin and Herrmann
The complexity of algorithm forms sharp contrast；
(2) present invention is when carrying out source wavelet estimation, it is only necessary to calculation formula 20 or 21 this step for being simple and efficient；And
The REPSI method of Lin and Herrmann calls least square QR- factorization algorithm when carrying out source wavelet estimation, needs repeatedly
Iteration；
(3) present invention limits characteristic according to the band of seismic data, in multiple wave, source wavelet estimation merely with earthquake number
According to the high dominant frequency band 0-70Hz of confidence level and signal-to-noise ratio (see Fig. 6 b)；And the method for Lin and Herrmann has used entire frequency
It is time-consuming to produce the calculating being largely not necessarily to by band 0-125Hz.
To sum up, we utilize same data, same computing resource environment, with Lin and Herrmann in 2013 Nian Suoti
REPSI method compares out, illustrates advantages of the present invention.
| Method | Runing time |
| The REPSI method of Lin and Herrmann | 21 points 6 seconds |
| The method of the present invention (normal higher-order spectra situation) | 5 points 57 seconds |
| The method of the present invention (becomes higher-order spectra situation) | 6 points 4 seconds |
1.5 data example of Pluto
In order to further show effectiveness of the invention, we are applied to another and have convictive Pluto
1.5 data.This data is generated by SMAART JV microstructure modeling in the world and is disclosed release, is more for open test in the world
One of the signal data of subwave algorithm.Source wavelet used is Ricker wavelet when the data forward simulation, every big gun wavelet
Equally namely Chang Zibo.Time sampling interval is 8 milliseconds.The original of the open release of SMAART JV tissue is directly utilized in we
Data do not do any pretreatment to initial data in addition to cutting off direct wave.Original input data from Fig. 6 a and Fig. 7 a
From the point of view of, there is a large amount of multiple wave energy in the data.The frequency range used in treatment process of the present invention is 0-50Hz.
From the point of view of processing result in Fig. 9 and Figure 10, in original input data Fig. 9 a and Figure 10 a, the same primary wave of multiple wave
It is mixed in together, interfere with each other, differentiate it is unclear, and the present invention realize accurate primary wave, multiple wave estimation and separate (especially
It is the signified place of arrow in Fig. 7).Due to when this digital simulation generates source wavelet used be do not change with big gun collection, so
It is more consistent that the present invention carries out the result generated when change higher-order spectra, sees Figure 11.
Figure 12 illustrates the scheme of another evaluation multiple wave estimation technique validity.The essence that Figure 12 is generated is feedback
Iterative model, it may be assumed that the primary wave (Figure 12 b) that the present invention estimates gets off to predict that single order surface is related more as secondary focus incomingly
Subwave (Figure 12 c), single order surface-related multiple (Figure 12 c) get off to predict that second-order surface is related more as secondary focus incomingly
Subwave (Figure 12 d), second-order surface related multiple (Figure 12 d) get off to predict that as secondary focus, three rank surfaces are related more incomingly
Subwave (Figure 12 e), three rank surface-related multiples (Figure 12 e) get off to predict that quadravalence surface is related more as secondary focus incomingly
Subwave (Figure 12 f), quadravalence surface-related multiple (Figure 12 f) get off to predict that as secondary focus, five rank surfaces are related more incomingly
Subwave (Figure 12 g), etc. and so on.For this data, (Figure 12 c is extremely for the surface-related multiple for the different orders predicted
Figure 12 g) it can be identified by identification easily in original input data (Figure 12 a).Primary wave estimated by the present invention, difference
The same original input data (Figure 12 a) of total data (Figure 12 h) that the cumulative summation of the surface-related multiple of order generates is intimate
It is completely the same, difference (Figure 12 i) very little between the two.
Figure 13 further illustrates the processing result that the present invention is applied to other big gun collection in 1.5 data of Pluto, can see
Primary wave and multiple wave are able to more accurately estimate and separate out.
The embodiment of the present invention is described with above attached drawing, but the invention is not limited to above-mentioned specific
Embodiment, the above mentioned embodiment is only schematical, rather than restrictive, those skilled in the art
Under the inspiration of the present invention, without breaking away from the scope protected by the purposes and claims of the present invention, it can also make very much
Form, all of these belong to the protection of the present invention.
Claims (6)
1. a kind of based on the surface-related multiple and higher-order spectra method that accelerate linear Bregman algorithm, which is characterized in that under including
State step:
S1, alternately update Green's function to be estimated and source wavelet to be estimated using following step S11 and S12 until
Reach preset times, wherein carrying out the source wavelet to be estimated of step S11 for the first time has preset initial value:
S11, according to the input data of source wavelet to be estimated and acquisition, using accelerate linear Bregman algorithm update to
The Green's function of estimation；
S12, according to updated Green's function to be estimated and the input data, updated using following formula to be estimated
Source wavelet；
When source wavelet to be estimated is not with big gun collection spatial position change:
When source wavelet to be estimated is with big gun collection spatial position change:
S2, according to after final updated as a result, obtaining the estimated result of at least one of source wavelet and surface-related multiple；
Wherein, subscript est is indicated to be estimated,The element on the diagonal line of source wavelet matrix slice is represented, j is indicated with big gun
Collect spatial position or big gun collection number, vec () indicates a matrix becoming column vector, subscriptRepresentative takes a frequency slice
The jth of matrix arranges, subscriptAll column for taking a frequency slice matrix are represented,For the matrix table of the frequency domain of input data
Up to the slice of formula,For the slice of the matrix expression of the frequency domain of source wavelet, rsurfFor the reflection of air and the contact surface of water
Coefficient,For the slice of the matrix expression of the frequency domain of Green's function；
It is to be estimated using accelerating linear Bregman algorithm to update according to source wavelet to be estimated and the input data of acquisition
Green's function specifically include:
According to formulaAnd
It is iterated, successively finds out y1、x1、y2、x2、y3…、Wherein,For this execute step S11 after wait estimate
The vector expression-form of the time-domain of the update result of the Green's function of meter, lmaxFor this execute step S11 when acceleration it is linear
The maximum number of iterations of Bregman algorithm；
In formula, x0、With y0It is equal to the vector expression-form of the time-domain of the preset initial value of Green's function to be estimated,
The vector expression-form for the time-domain that data vector b is the input data is observed,tlWhat is represented is step-length,
Its expression formula are as follows:
Mapping function ∏σ(Axl- b) is defined as:
Soft-threshold function shrink (x, μ) is defined as:
Shrink (x, μ)=max (| x |-μ, 0) sign (x),
Wherein, σ represents the noise level factor, and max () function representation is maximized, | | function representation takes absolute value, sign (x) table
Show that symbol value function, μ indicate threshold value；
Matrix operator A indicates to carry out the operational set of operations described below: by Green's function to be estimated by the vector table shape of time-domain
Formula, which carries out Fourier transformation, becomes the vector sheet form of frequency domain, by Green's function to be estimated by the vector sheet form of frequency domain
The matrix sheet form for becoming frequency domain, according to formulaCalculate separately the corresponding simulation of input data
Each slice of the expression matrix form of the frequency domain of dataAccording to each sliceInversefouriertransform is carried out to obtain
To the vector expression-form of the time-domain of the corresponding analogue data of input data.
2. surface-related multiple according to claim 1 and higher-order spectra method, which is characterized in that in the step S2,
The estimated result of source wavelet is obtained according to the result for updating source wavelet to be estimated for the last time；
The estimated result of surface-related multiple is according to m=rsurfG*d is obtained, and m is that the vector of the time-domain of surface-related multiple expresses shape
Formula, g are the vector expression-form of the time-domain of Green's function, and d is the vector expression-form of the time-domain of input data.
3. surface-related multiple according to claim 1 and higher-order spectra method, which is characterized in that if Green's letter to be estimated
Number has reached the maximum update times and source wavelet to be estimated is not up to maximum update times, then carries out step S11 again
When, lmaxThe maximum update times for being updated to source wavelet to be estimated subtract the difference after the number being updated, and before
L when step S11 is executed every timemaxIt is equal；If source wavelet to be estimated has reached the maximum update times and lattice to be estimated
Woods function is not up to maximum update times, then step S12 of subsequent progress, and l when each execution step S11maxPhase
Deng.
4. it is a kind of based on the surface-related multiple and higher-order spectra system that accelerate linear Bregman algorithm, characterized by comprising:
Alternately updating unit, for use following Green's function updating unit and source wavelet updating unit alternately update to
The Green's function of estimation and source wavelet to be estimated are until reach preset times, wherein calling Green's function updating unit for the first time
Source wavelet to be estimated have preset initial value:
Green's function updating unit, it is linear using accelerating according to source wavelet to be estimated and the input data of acquisition
Bregman algorithm updates Green's function to be estimated；
Source wavelet updating unit, according to updated Green's function to be estimated and the input data, using following public affairs
Formula updates source wavelet to be estimated；
When source wavelet to be estimated is not with big gun collection spatial position change:
When source wavelet to be estimated is with big gun collection spatial position change:
Estimated result acquiring unit, for according to after final updated as a result, obtaining source wavelet to be estimated and to be estimated
At least one of surface-related multiple estimated result；
Wherein, subscript est is indicated to be estimated,The element on the diagonal line of source wavelet matrix slice is represented, j is indicated with big gun
Collect spatial position or big gun collection number, vec () indicates a matrix becoming column vector, subscriptRepresentative takes a frequency slice
The jth of matrix arranges, subscriptAll column for taking a frequency slice matrix are represented,For the matrix table of the frequency domain of input data
Up to the slice of formula,For the slice of the matrix expression of the frequency domain of source wavelet, rsurfFor the reflection of air and the contact surface of water
Coefficient,For the slice of the matrix expression of the frequency domain of Green's function；
It is to be estimated using accelerating linear Bregman algorithm to update according to source wavelet to be estimated and the input data of acquisition
Green's function specifically include:
According to formulaAnd
It is iterated, successively finds out y1、x1、y2、x2、y3…、Wherein,Green's function, which is executed, for this updates list
The vector expression-form of the time-domain of the update result of Green's function to be estimated after member, lmaxGreen's function is executed for this
The maximum number of iterations of the linear Bregman algorithm of acceleration when updating unit；
In formula, x0、With y0It is equal to the vector expression-form of the time-domain of the preset initial value of Green's function to be estimated,
The vector expression-form for the time-domain that data vector b is the input data is observed,tlWhat is represented is step-length,
Its expression formula are as follows:
Mapping function ∏σ(Axl- b) is defined as:
Soft-threshold function shrink (x, μ) is defined as:
Shrink (x, μ)=max (| x |-μ, 0) sign (x),
Wherein, σ represents the noise level factor, and max () function representation is maximized, | | function representation takes absolute value, sign (x) table
Show that symbol value function, μ indicate threshold value；
Matrix operator A indicates to carry out the operational set of operations described below: by Green's function to be estimated by the vector table shape of time-domain
Formula, which carries out Fourier transformation, becomes the vector sheet form of frequency domain, by Green's function to be estimated by the vector sheet form of frequency domain
The matrix sheet form for becoming frequency domain, according to formulaCalculate separately the corresponding simulation of input data
Each slice of the expression matrix form of the frequency domain of dataAccording to each sliceInversefouriertransform is carried out to obtain
To the vector expression-form of the time-domain of the corresponding analogue data of input data.
5. surface-related multiple according to claim 4 and higher-order spectra system, which is characterized in that the estimated result obtains
In unit,
The estimated result of source wavelet is obtained according to the result for updating source wavelet to be estimated for the last time；
The estimated result of surface-related multiple is according to m=rsurfG*d is obtained, and m is that the vector of the time-domain of surface-related multiple expresses shape
Formula, g are the vector expression-form of the time-domain of Green's function, and d is the vector expression-form of the time-domain of input data.
6. surface-related multiple according to claim 4 and higher-order spectra system, which is characterized in that if Green's letter to be estimated
Number has reached the maximum update times and source wavelet to be estimated is not up to maximum update times, then carries out Green's function again
When updating unit, lmaxThe maximum update times for being updated to source wavelet to be estimated subtract the difference after the number being updated
Value, and before every time execute Green's function updating unit when lmaxIt is equal；If source wavelet to be estimated has reached the maximum more
New number and Green's function to be estimated is not up to maximum update times, then source wavelet updating unit of subsequent progress,
And every time execute Green's function updating unit when lmaxIt is equal.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710525466.7A CN107367760B (en) | 2017-06-27 | 2017-06-27 | Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710525466.7A CN107367760B (en) | 2017-06-27 | 2017-06-27 | Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107367760A CN107367760A (en) | 2017-11-21 |
| CN107367760B true CN107367760B (en) | 2019-04-02 |
Family
ID=60305797
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710525466.7A Expired - Fee Related CN107367760B (en) | 2017-06-27 | 2017-06-27 | Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107367760B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110858000B (en) * | 2018-08-24 | 2021-07-02 | 中国石油天然气股份有限公司 | Seismic data reconstruction method and device, computer equipment and storage medium |
| CN111694057B (en) * | 2020-06-03 | 2021-03-23 | 西安交通大学 | Method, storage medium and equipment for suppressing surge noise of seismic data |
| CN111694056B (en) * | 2020-06-03 | 2021-03-02 | 西安交通大学 | Method, storage medium and equipment for suppressing abnormal noise of seismic data |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103278849A (en) * | 2013-05-24 | 2013-09-04 | 中国石油天然气集团公司 | Method and system for performing wavelet estimation on the basis of seismic data and logging information |
| CN105334537A (en) * | 2015-10-26 | 2016-02-17 | 中国石油大学(华东) | Primary wave and multiple wave separation method based on alternative splitting Bregman iterative algorithm |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101182839B1 (en) * | 2010-08-26 | 2012-09-14 | 서울대학교산학협력단 | Method and Apparatus for Time domain Reverse Time Migration with Source Estimation |
-
2017
- 2017-06-27 CN CN201710525466.7A patent/CN107367760B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103278849A (en) * | 2013-05-24 | 2013-09-04 | 中国石油天然气集团公司 | Method and system for performing wavelet estimation on the basis of seismic data and logging information |
| CN105334537A (en) * | 2015-10-26 | 2016-02-17 | 中国石油大学(华东) | Primary wave and multiple wave separation method based on alternative splitting Bregman iterative algorithm |
Non-Patent Citations (4)
| Title |
|---|
| "A linearized Bregman method for compressive waveform inversion";Xintao Chai 等;《SEG International Exposition and 86th Annual Meeting》;20161231;第1449-1454页 |
| "SVD加速的线性Bregman算法";孙涛 等;《计算机应用研究》;20140731;第31卷(第7期);第2001-2003页 |
| "一种快速复数线性Bregman迭代算法及其在ISAR成像中的应用";李少东 等;《中国科学:信息科学》;20151231;第45卷(第9期);第1179-1196页 |
| "动态地震子波估计";彭才 等;《大庆石油地质与开发》;20071031;第26卷(第5期);第125-128页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107367760A (en) | 2017-11-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chen et al. | Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization | |
| Xue et al. | Amplitude-preserving iterative deblending of simultaneous source seismic data using high-order radon transform | |
| CN107367760B (en) | Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm | |
| Zhou | A POCS method for iterative deblending constrained by a blending mask | |
| NZ575497A (en) | Iterative inversion of data from simultaneous geophysical sources | |
| Zhou et al. | Robust noise attenuation based on nuclear norm minimization and a trace prediction strategy | |
| Zhou et al. | Sparse dictionary learning for seismic noise attenuation using a fast orthogonal matching pursuit algorithm | |
| US20220342103A1 (en) | Noise Attenuation Methods Applied During Simultaneous Source Deblending and Separation | |
| CN111766628A (en) | Preconditioned time domain elastic medium multi-parameter full waveform inversion method | |
| Yang et al. | Mini-batch optimized full waveform inversion with geological constrained gradient filtering | |
| CN107390261B (en) | Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithm | |
| CN112882115A (en) | Magnetotelluric signal denoising method and system based on GWO optimized wavelet threshold | |
| Zhou et al. | Deblending of simultaneous-source data using iterative seislet frame thresholding based on a robust slope estimation | |
| Aghamiry et al. | Complex-valued imaging with total variation regularization: an application to full-waveform inversion in visco-acoustic media | |
| Doulgeris et al. | Separation of blended impulsive sources using an iterative approach | |
| Bai et al. | Seismic signal enhancement based on the low‐rank methods | |
| Kontakis et al. | Using a hybrid focal: curvelet transform for deblending | |
| Bai et al. | Iterative deblending of simultaneous-source data using smoothed singular spectrum analysis | |
| Bai et al. | Nonstationary least-squares decomposition with structural constraint for denoising multi-channel seismic data | |
| Dondurur et al. | Swell noise suppression by Wiener prediction filter | |
| Zhang et al. | A robust method for random noise suppression based on the Radon transform | |
| Krebs et al. | Fast full wave seismic inversion using source encoding | |
| Yuan et al. | Ground roll attenuation based on an empirical curvelet transform | |
| Wu et al. | Iterative deblending based on the modified singular spectrum analysis | |
| Lopez et al. | SRME and estimation of primaries by sparse inversion: a hybrid approach |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190402 Termination date: 20200627 |