CN107390261A - Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithms - Google Patents
Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithms Download PDFInfo
- Publication number
- CN107390261A CN107390261A CN201710524014.7A CN201710524014A CN107390261A CN 107390261 A CN107390261 A CN 107390261A CN 201710524014 A CN201710524014 A CN 201710524014A CN 107390261 A CN107390261 A CN 107390261A
- Authority
- CN
- China
- Prior art keywords
- estimated
- function
- msub
- mrow
- source wavelet
- 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
- 238000000034 method Methods 0.000 title claims abstract description 78
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 64
- 238000001228 spectrum Methods 0.000 title claims abstract description 39
- 230000008859 change Effects 0.000 claims abstract description 26
- 239000011159 matrix material Substances 0.000 claims description 63
- 238000013507 mapping 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 17
- 230000003595 spectral effect Effects 0.000 abstract 1
- 230000006870 function Effects 0.000 description 58
- 238000010586 diagram Methods 0.000 description 26
- 238000012545 processing Methods 0.000 description 10
- 241001323321 Pluto Species 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 5
- 230000003044 adaptive effect Effects 0.000 description 3
- 235000013399 edible fruits Nutrition 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000010410 layer Substances 0.000 description 3
- 239000003208 petroleum Substances 0.000 description 3
- 241001269238 Data Species 0.000 description 2
- 230000001186 cumulative effect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 239000013535 sea water Substances 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 238000012360 testing method Methods 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
- YTAHJIFKAKIKAV-XNMGPUDCSA-N [(1R)-3-morpholin-4-yl-1-phenylpropyl] N-[(3S)-2-oxo-5-phenyl-1,3-dihydro-1,4-benzodiazepin-3-yl]carbamate Chemical compound O=C1[C@H](N=C(C2=C(N1)C=CC=C2)C1=CC=CC=C1)NC(O[C@H](CCN1CCOCC1)C1=CC=CC=C1)=O YTAHJIFKAKIKAV-XNMGPUDCSA-N 0.000 description 1
- 230000001154 acute effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 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
- 238000007405 data analysis Methods 0.000 description 1
- 238000013079 data visualisation 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
- 230000005284 excitation Effects 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 238000003384 imaging method Methods 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
- 238000004519 manufacturing process Methods 0.000 description 1
- 229940050561 matrix product Drugs 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000004080 punching Methods 0.000 description 1
- 238000012797 qualification Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000007787 solid Substances 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
- 210000003462 vein Anatomy 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/003—Seismic data acquisition in general, e.g. survey design
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Geophysics (AREA)
- Acoustics & Sound (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Geology (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Physics (AREA)
- Algebra (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Complex Calculations (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a kind of surface-related multiple and higher-order spectra method and system based on linear Bregman algorithms, feedback model based on description primary wave and surface-related multiple relation, this method carries out sparse inversion first with linear Bregman algorithms, to obtain the impulse response unrelated with surface, then the source wavelet using step linear spectral estimation formulas renewal with big gun collection spatial position change, and the impulse response unrelated with surface and source wavelet are updated by cycle alternation, final estimation obtains primary wave and surface-related multiple.Surface-related multiple method of estimation compared to traditional multi-parameter, based on Chang Zibo hypothesis and complicated algorithm, implement the method for estimation and system of the present invention, the source wavelet with big gun collection spatial position change, and the efficiency high of method of estimation and system, the precision height for taking less, estimating result can be estimated.
Description
Technical field
In terms of the present invention relates to seismic data processing in field of seismic exploration, more particularly to seismic exploration technique, more specifically
Ground is said, is related to a kind of surface-related multiple and higher-order spectra method and system based on linear Bregman algorithms.
Background technology
Multiple reflection (more subwaves) relevant issues are most one of distinct issues of generally existing in reflection shooting.
Because 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 aboundresources, it is the focus and focus of current seafari.Produced in current petroleum industry
In, a basic assumption of conventional seismic data processing technology is：The reflectance data body of input only by primary reflection (once
Ripple) composition.More subwaves can be mistakenly considered the part for being primary wave or being blended in primary wave.More subwaves are so as to one
The processing and imaging of subwave cause to disturb strongly, influence resource exploration, so more subwave estimation techniques have fabulous actual work
Industry application prospect.
The estimation of more subwaves or drawing method that industrial quarters uses at present are broadly divided into following a few classes：One kind using more subwaves with
The feature difference multiple suppression of primary wave, referred to as filter method；One kind predicts more subwaves from geological data first, then will
It is adaptively subtracted from initial data, referred to as predicts subtractive method；Also a kind of inverse-scattering series method multiple suppression.
Filter method by more subwaves and primary wave in particular transform domain have obvious feature difference premised on it is assumed that passing through
Various mathematical algorithms are designed to obtain this species diversity and Multiple attenuation.The principal character difference that filter method utilizes includes more subwaves
Periodicity and separability.Class method is filtered in the multiple attenuation applied to some complicated structures, more subwaves can be run into
Feature difference very little between primary wave does not either have discrepant situation, at this moment filters the applicable precondition of class method and obtains
Less than satisfaction, multiple wave pressure DeGrain.
Back scattering series predicts more subwaves, and there is also some problems：It is computationally intensive, cost is high；It is difficult to predict remote offset distance
More subwaves；The predictive ability of more subwaves also has certain limitation in complex dielectrics.Among these, calculate cost height and cause this method very
Hardly possible is applied among the processing of real data.
The feedback model that Berkhout proposes the complicated more subwaves of description in nineteen eighty-two is theoretical, and it is multiple to have established feedback iteration
The mathematical physics basis of ripple drawing method, but source wavelet known to needs during Multiple attenuation.In order to illustrate this problem, Yi Jibian
In follow-up 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 the input data collected,For D after Fourier's direct transform frequency
One section of the matrix expression in rate domain, S are the matrix expression of the time-domain of source wavelet to be estimated,It is S through in Fu
One of the matrix expression of frequency domain section after leaf direct transform, G are that (namely the underground pulse without more subwaves rings Green's function
Should) time-domain matrix expression,For G me is cut into slices after Fourier's direct transform for one of the matrix expression of frequency domain
Represent the data of frequency domain or variable, r with ^ symbolssurfFor air and the reflectance factor of the contact surface of water, it is assumed that rsurf
=-1, whereinWithFrequency domain product (namely time domain convolution) represent estimation primary wave frequency section, namelyWhereinFor a section of the frequency domain of primary wave.Involved primary wave refers to except surface is more in the present invention
All ripples outside subwave, namely primary wave here contain interbed multiple (it is unrelated with surface).
Represent a section of the frequency domain of the surface-related multiple of estimation, namely dataPropagated downwards as secondary focus, same to underground
Impulse response Green's functionFrequency domain product (time domain convolution) predict the section of surface-related multiple frequency domain.By
The formula of equation 1 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, and it is based on one by introducing
Adaptive the subtracting each other of the minimum hypothesis of secondary wave energy, estimates source wavelet from real data.In order to avoid due to solving focus
Ripple and caused nonlinear problem, Berkhout and Verschuur proposed iteration SRME methods in 1997.Wave energy is most
This small hypotheses condition is not that can meet in all cases.In addition, in actual applications, SRME is also present
Adaptive prediction subtracts each other the problem of infringement useful signal.
The problem of infringement useful signal is subtracted each other in prediction, van Groenestijn and Verschuur be present to improve SRME
Sparse inversion primary wave estimation technique, namely EPSI were proposed in 2009.In order to improve problem present in EPSI, Lin and
Herrmann proposed sane sparse inversion primary wave estimation technique, namely REPSI in 2013.EPSI and REPSI be with
The more related prior art of the present invention, their technical scheme and the problem of exist, shortcoming will be explained in detail below
State.
Further, since during actual field earthquake data acquisition, the source wavelet excited per big gun is not quite similar, and current institute
The more subwave methods of estimation having assume that all source wavelets are identical, influence more subwave estimated accuracies, and the present invention can estimate
With the source wavelet of big gun collection spatial position change.
The computing resource (internal memory, processor) of each iteration consuming is more when being estimated due to more subwaves, it is longer to calculate the time,
In addition to improving the estimated accuracy of more subwaves, source wavelet, the efficiency, the rate of convergence of boosting algorithm, simplified calculation of method are improved
The complexity of method, the adjustment parameter for reducing method, and the development trend of more subwave estimation techniques.
Prior art related to the present invention is described below.
The technical scheme of prior art one
Van Groenestijn and Verschuur proposed sparse inversion one in 2009 in SRME theoretical foundation
Subwave method of estimation EPSI, this method, using fitting data residual error as driving, are handed over during iterative inversion by steepest descent method
For renewal 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 adaptive the step of subtracting each other, and preferably protects one
Subwave useful signal.
Van Groenestijn and Verschuur EPSI methods 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 ensure underground impulse response G in the openness of time domain, EPSI during impulse response G is updated, according to
Time is ascending, amplitude is by by force to weak sequential update pulse.Namely in order to promote underground impulse response G time domain sparse
Property, EPSI methods devise numerous artificial adjustment parameters, for example, EPSI it needs to be determined that：One time window, for judging Green's letter
Several more new ranges (namely being updated in time window, the zero setting outside time window)；Reflected during each iteration in Green's function
The number of lineups, etc..As pointed by Lin and Herrmann in paper in 2013：Even if it has been accurately known that underground is anti-
The number (we can not can know that in practice) in firing area face, EPSI methods can not also make full use of this prior information.EPSI
Underground impulse response renewal, easily influenceed by the parameter such as pulse sliced window and pickup quantity, and then influence more subwaves
Estimated accuracy.Intrinsically, the degree of rarefication of Green's function is controlled by artificially in EPSI methods, rather than by counting
Optimized algorithm is learned to control.This shortcoming limits application of the EPSI methods in actual industrial.
Prior art two related to the present invention
Theoretical (namely the formula of feedback model 1) based on EPSI identicals, Lin and Herrmann are in 2013 by EPSI undergrounds arteries and veins
The sparse constraint condition of punching response is revised as a norm constraint, and using the spectrum mapping norm algorithm of gradient oneTo solve
To sparse underground impulse response Green's function, the stability of original EPSI methods is improved to a certain extent.Lin and
Herrmann was named as Robust EPSI (REPSI) in proposed method in 2013.What Lin and Herrmann was proposed
REPSI methods, it is the same with EPSI, and alternately renewal underground impulse response G and two unknown parameters of source wavelet S, Jin Erchong
Structure primary waveAnd its corresponding surface-related multipleEPSI and REPSI maximum difference exists
In：The degree of rarefication majority of Green's function is controlled by artificial adjustment parameter in EPSI, and REPSI is then employedMathematics
Optimized algorithm solves to obtain sparse underground impulse response Green's function.
Lin and Herrmann REPSI methods also assume 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 algorithms
In generation, solves and obtains source wavelet, considerably increases the amount of calculation of method, reduces efficiency.And a step of the invention can be estimated to obtain
Source wavelet.
Another significant drawbacks of Lin and Herrmann REPSI methods are：It relies on an extremely complex spectrum mapping
The norm algorithm of gradient oneTo solve to obtain sparse underground impulse response Green's function.It is one to increase income
MATLAB programs, can be by accessing https://github.com/mpf/spgl1 is downloaded.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 simultaneously reduce the efficiency of REPSI methods, hinder the reality of REPSI methods in the industry should
With.
The content of the invention
The technical problem to be solved in the present invention is, for technology existing for above-mentioned EPSI in the prior art and REPSI
Defect, there is provided a kind of surface-related multiple and higher-order spectra method and system based on linear Bregman algorithms.
According to the wherein one side of the present invention, the present invention is its technical problem of solution, there is provided one kind is based on linear
The surface-related multiple and higher-order spectra method of Bregman algorithms, include following step：
S1, using following step S11 and S12 alternately update Green's function to be estimated and source wavelet to be estimated
Until reach preset times, wherein the source wavelet to be estimated for carrying out step S11 first has default initial value：
S11, the input data according to source wavelet to be estimated and collection, treated using the renewal of linear Bregman algorithms
The Green's function of estimation；
S12, according to the Green's function to be estimated after renewal and the input data, updated using following formula and wait to 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 the result after final updated, draw the estimation knot of at least one of source wavelet and surface-related multiple
Fruit；
Wherein, subscript est is represented to be estimated,The element on the diagonal of source wavelet matrix section is represented, j is represented
Numbered with big gun collection locus or big gun collection, vec () represents a matrix being changed into column vector, subscriptRepresentative takes a frequency
The jth row of section matrix, subscriptAll row for taking a frequency slice matrix are represented,For the square of the frequency domain of input data
The section of battle array expression formula,For the section of the matrix expression of the frequency domain of source wavelet, rsurfFor the contact surface of air and water
Reflectance factor,For the section of the matrix expression of the frequency domain of Green's function.
Further, in the surface-related multiple and higher-order spectra method based on linear Bregman algorithms of the present invention, root
According to source wavelet to be estimated and the input data of collection, Green's function to be estimated is updated using linear Bregman algorithms
Specifically include：
According to formula yl+1=yl-tlAHΠσ(Axl- b) and xl+1=shrink (yl+1, μ) it is iterated, obtain successively
y1、 x1、y2、x2、y3…、Wherein,For the renewal result of the Green's function to be estimated after this execution step S11
Time-domain vectorial expression-form, lmaxFor this execution step S11 when linear Bregman algorithms greatest iteration time
Number；
In formula, x0With y0The vectorial expression-form of the time-domain of the default initial value of Green's function to be estimated is equal to,
Observe the vectorial expression-form for the time-domain that data vector b is the input data, tlWhat is represented is step-length, and its expression formula is：
Mapping function Пσ(Axl- b) definition be：
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 takes maximum, | | function representation takes absolute value, sign
(x) symbol value function is represented；
Matrix operator A represents to carry out the operational set of operations described below：Vector by Green's function to be estimated by time-domain
Sheet form carries out the vectorial sheet form that Fourier transformation is changed into frequency domain, the vector table by Green's function to be estimated by frequency domain
Form is changed into the matrix sheet form of frequency domain, according to formulaCalculate respectively corresponding to input data
Each section of the expression matrix form of the frequency domain of analogue dataAccording to each sectionCarry out anti-Fourier's change
Get the vectorial expression-form of the time-domain of analogue data corresponding to input data in return.
Further, in the surface-related multiple and higher-order spectra method based on linear Bregman algorithms of the present invention, step
In rapid S2,
The estimated result of source wavelet obtains 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 vectorial expression-form of the time-domain of Green's function, and d is the vectorial expression-form of the time-domain of input data.
Further, in the surface-related multiple and higher-order spectra method based on linear Bregman algorithms of the present invention, if
Green's function to be estimated has been maxed out update times and source wavelet to be estimated is not up to maximum update times, then again
During secondary progress step S11, lmaxBe updated to source wavelet to be estimated maximum update times subtract updated number after
Difference, and before every time perform step S11 when lmaxIt is equal；If source wavelet to be estimated has been maxed out renewal time
Several and to be estimated Green's functions is not up to maximum update times, then follow-up only to carry out a step S12, and performs step every time
L during S11maxIt is equal.
The present invention for solve its technical problem, additionally provide a kind of surface-related multiple based on linear Bregman algorithms and
Higher-order spectra system, comprising：
Alternately updating block, for using following Green's function updating blocks and source wavelet updating block alternating more
New Green's function to be estimated is with source wavelet to be estimated until reaching preset times, wherein calling Green's function renewal first
The source wavelet to be estimated of unit has default initial value：
Green's function updating block, according to source wavelet to be estimated and the input data of collection, using linear
Bregman algorithms update Green's function to be estimated；
Source wavelet updating block, according to the Green's function to be estimated after renewal and the input data, under
State formula and update 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 the result after final updated, drawing source wavelet to be estimated and treating
The estimated result of at least one of the surface-related multiple of estimation；
Wherein, subscript est is represented to be estimated,The element on the diagonal of source wavelet matrix section is represented, j is represented
Numbered with big gun collection locus or big gun collection, vec () represents a matrix being changed into column vector, subscriptRepresentative takes a frequency
The jth row of section matrix, subscriptAll row for taking a frequency slice matrix are represented,For the square of the frequency domain of input data
The section of battle array expression formula,For the section of the matrix expression of the frequency domain of source wavelet, rsurfFor the contact surface of air and water
Reflectance factor,For the section of the matrix expression of the frequency domain of Green's function.
Further, in the surface-related multiple and higher-order spectra system based on linear Bregman algorithms of the present invention, root
According to source wavelet to be estimated and the input data of collection, Green's function to be estimated is updated using linear Bregman algorithms
Specifically include：
According to formula yl+1=yl-tlAHПσ(Axl- b) and xl+1=shrink (yl+1, μ) it is iterated, obtain successively
y1、 x1、y2、x2、y3…、Wherein,For the renewal knot of the Green's function to be estimated after this execution step S11
The vectorial expression-form of the time-domain of fruit, lmaxFor this execution step S11 when linear Bregman algorithms greatest iteration time
Number；
In formula, x0With y0The vectorial expression-form of the time-domain of the default initial value of Green's function to be estimated is equal to,
Observe the vectorial expression-form for the time-domain that data vector b is the input data, tlWhat is represented is step-length, and its expression formula is：
Mapping function Πσ(Axl- b) definition be：
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 takes maximum, | | function representation takes absolute value, sign
(x) symbol value function is represented；
Matrix operator A represents to carry out the operational set of operations described below：Vector by Green's function to be estimated by time-domain
Sheet form carries out the vectorial sheet form that Fourier transformation is changed into frequency domain, the vector by Green's function to be estimated by frequency domain
Sheet form is changed into the matrix sheet form of frequency domain, according to formulaIt is corresponding that input data is calculated respectively
Analogue data frequency domain expression matrix form each sectionAccording to each sectionCarry out anti-Fourier
Conversion obtains the vectorial expression-form of the time-domain of analogue data corresponding to input data.
Further, in the surface-related multiple and higher-order spectra system based on linear Bregman algorithms of the present invention, estimate
Count in result acquiring unit,
The estimated result of source wavelet obtains 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 vectorial expression-form of the time-domain of Green's function, and d is the vectorial expression-form of the time-domain of input data.
Further, in the surface-related multiple and higher-order spectra system based on linear Bregman algorithms of the present invention, if
Green's function to be estimated has been maxed out update times and source wavelet to be estimated is not up to maximum update times, then
When carrying out step S11 again, lmaxThe maximum update times for being updated to source wavelet to be estimated subtract updated time
Difference after number, and before every time perform step S11 when lmaxIt is equal；If source wavelet to be estimated has been maxed out more
New number and Green's function to be estimated is not up to maximum update times, then it is follow-up only to carry out a step S12, and hold every time
L during row step S11maxIt is equal.
In the method for estimation and system based on linear Bregman algorithms for implementing the present invention, have the advantages that：This
The method of estimation and system of invention can estimate the source wavelet with big gun collection spatial position change, and method of estimation and system
Efficiency high, the precision height for taking less, estimating result.
Brief description of the drawings
Below in conjunction with drawings and Examples, the invention will be further described, in accompanying drawing：
Fig. 1 is linear Bregman Algorithm for Solving Ax=b of the invention to obtain the flow of the preferred embodiments of sparse solution x mono-
Figure；
Fig. 2 is the flow chart of the estimation source wavelet and the preferred embodiment of surface-related multiple one of the method for estimation of the present invention；
The schematic diagram data for including more subwaves of Fig. 3 (a) inputs；
Fig. 3 (b) is the primary wave schematic diagram data that the present invention estimates；
Fig. 3 (c) is the surface-related multiple schematic diagram data that the present invention estimates；
Fig. 4 (a) is that zero big gun of data in Fig. 3 (a) examines diagrammatic cross-section；
Fig. 4 (b) is the surface-related multiple schematic diagram data that the present invention estimates；
Fig. 4 (c) is the primary wave schematic diagram data of one embodiment of the invention estimation；
The primary wave schematic diagram data that the REPSI methods that Fig. 4 (d) is Lin and Herrmann are estimated；
Fig. 5 (a) is the schematic diagram of the source wavelet of estimation；
Fig. 5 (b) is the frequency domain amplitude spectrum schematic diagram of the source wavelet with the change of big gun collection of one embodiment of the invention estimation；
Fig. 6 (a) is the Pluto 1.5 for being used to test multiple ripple drawing method of the open release of international SMAART JV tissues
Schematic diagram data；
Fig. 6 (b) is the primary wave schematic diagram data that the present invention estimates；
Fig. 6 (c) is the surface-related multiple schematic diagram data that the present invention estimates；
Fig. 7 (a) is that zero big gun of data in Fig. 6 (a) examines diagrammatic cross-section；
Fig. 7 (b) is the primary wave schematic diagram data that the present invention estimates；
The surface-related multiple schematic diagram data that Fig. 7 (c) present invention estimates；
Fig. 8 (a) is the source wavelet schematic diagram changed with big gun collection estimated from the data of Pluto 1.5；
Fig. 8 (b) is that the frequency domain amplitude spectrum for the source wavelet with the change of big gun collection estimated from the data of Pluto 1.5 is illustrated
Figure；
Fig. 9 (a) is the data diagrammatic cross-section of the middle big gun of data in Fig. 7 (a)；
Fig. 9 (b) is the primary wave schematic diagram data that the present invention estimates；
Fig. 9 (c) is the single order surface-related multiple schematic diagram data that the present invention estimates；
Fig. 9 (d) is the multiple ripple schematic diagram data of second-order surface that the present invention estimates；
Fig. 9 (e) is the three rank surface-related multiple schematic diagram datas that the present invention estimates；
Fig. 9 (f) is the quadravalence surface-related multiple schematic diagram data that the present invention estimates；
Fig. 9 (g) is the five rank surface-related multiple schematic diagram datas that the present invention estimates；
Fig. 9 (h) is the schematic diagram of Fig. 9 (b) cumulative summations of data into Fig. 9 (g)；
Fig. 9 (i) is Fig. 9 (a) and Fig. 9 (h) poor schematic diagram；
Schematic diagrames of the Figure 10 (a) for big gun number in the data of Pluto 1.5 between 701-880 zero big gun inspection office section；
Figure 10 (b) is the primary wave schematic diagram data that the present invention estimates；
Figure 10 (c) is the surface-related multiple schematic diagram data that the present invention estimates.
Noun explanation
Primary wave (Primaries)：The seismic wave for referring to excite only occurs one during earth-layer propagation in subsurface interface
Secondary reflection is that the wave detector for being embedded in earth's surface is received, referred to as primary reflection, abbreviation primary wave.It is involved in the present invention
Primary wave refer to all ripples in addition to surface-related multiple, namely primary wave here contains interbed multiple (itself and table
Face is unrelated).
More subwaves (Multiples)：The seismic wave for referring to excite 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 received, referred to as multiple reflection is simple
Claim more subwaves.
Surface-related multiple (Surface-related multiples) or surface-related multiple (Surface
multiples)：Refer to the multiple reflection related 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 (between stratum) production
Raw multiple reflection.
Source wavelet：Refer to ripple caused by epicenter excitation.
Feedback model：Feedback model.
Surface-Related Multiple Elimination(SRME)：Surface-related multiple is suppressed.
Estimation of Primaries by Sparse Inversion(EPSI)：Sparse inversion primary wave is estimated.
Robust Estimation of Primaries by Sparse Inversion(EPSI)：Sane is sparse anti-
Drill primary wave estimation.
Spectral-Projected GradientThe spectrum mapping norm algorithm of gradient one.
LSQR：Least-squares QR, least square QR- factorization algorithm.
Linearized Bregman(LB)：Linear Bregman algorithms.
Accelerated Linearized Bregman(ALB)：The linear Bregman algorithms accelerated.
MATLAB：Be MathWorks companies of the U.S. produce business mathematics software, for algorithm development, data visualization,
Data analysis and the advanced techniques computational language and interactive environment of numerical computations.MATLAB is matrix＆laboratory
Two contaminations, mean matrix factory (matrix labotstory).
Signal-to-Noise Ratios(SNR)：Signal to noise ratio.
Embodiment
In order to which technical characteristic, purpose and the effect of the present invention is more clearly understood, now compares accompanying drawing and describe in detail
The embodiment of the present 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 vectorial expression-form for representing time-domain, and S, G, D, P, M are the square of time-domain
Battle array expression-form.For the section of the expression matrix form of S, G, D, P, M frequency domain, subscript est
Represent to be estimated,For the vectorial expression-form of s, g, d, p, m frequency domain, S, G, D, P, M are three-dimensional
Matrix, and its section is then two-dimensional matrix, each three-dimensional matrice can be split as multiple two-dimensional matrixs.Any one two dimension
Matrix Z=[Z1Z2Z3…Zn] vectorial expression-form can be expressed as：
Wherein Z1To ZnTwo-dimensional matrix Z row are represented, conversely, expression matrix can be obtained by vectorial expression-form
Form.
The main formulas that primary wave, surface-related multiple and source wavelet estimation utilize is the mathematical physics relation of feedback model
Formula, namely formula (1).The source wavelet s excited is passed to underground, and time domain convolution is carried out with underground impulse response Green's function g
(asterisk * represents convolution), produces primary wave p, and mathematical relationship expression formula is：
P=g*s (2)
Caused primary wave p runs into seawater and reflected with air contact surfaces, propagates downwards, and then is passed again as focus
Enter underground, carry out time domain convolution with underground impulse response g, produce the more subwave m of single order1
m1=g*rsurfP=rsurfg*g*s (3)
The caused more subwave m of single order1Run into Free Surface to reflect, and then underground is passed to again as focus, with ground
Lower impulse response g carries out time domain convolution, produces second-order multiples m2
m2=g*rsurfm1=(rsurf)2g*g*g*s (4)
The surface-related multiple of higher order time produces formula by that analogy.
Observe the summation that obtained total data d is primary wave p and all order surface-related multiples, its mathematical expression
Formula is：
And then have
Above-mentioned feedback model mathematical physics relational expression can be summarized as：
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 ease of contrast, by its transcription in this：
The main object of the present invention is：According to expression formula 7-8, known variables are solved：Source wavelet s and Green's function g, enters
And the primary wave p=g*s and surface-related multiple m=r estimatedsurfg*d。
In the case of the source wavelet of the given estimation of present invention explanation first, inverting Sparse Pulse responds the side of Green's function
Method.
In the source wavelet s of given estimationestIn the case of, according to expression formula (7)-(8), we write out inverting gestWhen institute
Formula, namely：
HereIt is a representational forward modeling operator of simplification, its core is in calculation formula 8
Specifically, matrix operatorRepresent the operational set of progress operations described below：By Green's function to be estimated from time-domain to
Scale form g carries out the vectorial sheet form that Fourier transformation is changed into frequency domainBy Green's function to be estimated by frequency domain
Vectorial sheet formIt is changed into the matrix sheet form of frequency domainAccording to formulaInput is calculated respectively
Each section of the expression matrix form of the frequency domain of analogue data corresponding to dataAccording to the input data calculated
Each section of the expression matrix form of the frequency domain of corresponding analogue dataInversefouriertransform is carried out to be inputted
The vectorial expression-form of the time-domain of analogue data corresponding to data
For the ease of being expanded on further how the present invention solves to obtain sparse Green's function g, 9 formulas are rewritten as to following mark
Quasi- form：
Ax=b (10)
Wherein,X=g, b=d.
In order to seek Green's function g sparse solution, the present invention utilizes the linear formula of Bregman Algorithm for Solving 10.Linearly
Bregman algorithms were proposed that Yin has done detailed analysis in 2010 to linear Bregman algorithms by Yin etc. in 2008.
Cai is studied in the paper of 2009 shows that linear Bregman algorithms are intended to solve：
Expression solves x and causes expression formula to obtain minimum value, and subject to are represented under qualifications.
Wherein, μ >=0 is degree of rarefication controlling elements, | | x | |1With | | x | |2Difference representation vectorNorm andNorm, σ
Represent the noise level factor in data.Following algorithm 1 gives the detailed iterative step of linear Bregman algorithms.
With reference to figure 1, Fig. 1 be the present invention linear Bregman Algorithm for Solving Ax=b to obtain sparse solution x flow chart,
Namely the flow chart of algorithm 1, algorithm 1 specifically include following step：
A) parameters described below is obtained：Forward modeling operator A, observe data vector b, initial value vector x0, threshold value μ, greatest iteration time
Number lmax。
B) initialize:Iteration variable l=0, by x0It is assigned to y0
C) following each d are repeated to f steps lmaxIt is secondary
D) y is calculatedl+1=yl-tlAHПσ(Axl-b)
E) x is calculatedl+1=shrink (yl+1,μ)
F) l is updated to l+1
G) sparse solution is drawnSparse solution x is that this carries out linear Bregman algorithms more
New gestResult afterwards.
In above-mentioned algorithm 1, step d tlWhat is represented is step-length, and its expression formula is：
Wherein, subscript H represents complex conjugate transposition.Mapping function Πσ(Axl- b) be for consider data present in noise, its
It is defined as：
And soft-threshold function shrink (x, μ) is defined as：
Shrink (x, μ)=max (| x |-μ, 0) sign (x) (14)
Wherein, max () function representation takes maximum, | | function representation takes absolute value, and sign (x) represents symbol value
Function.
The present invention is following will to be illustrated, in the case where inverting obtains Sparse Pulse response Green's function, to estimate source wavelet.
When inverting obtains gestAfterwards, according to formula 8, we estimate to solve the problems, such as during source wavelet be：
NowFor unknown number to be solved.We carry out source wavelet with simple derives to explain 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 i.e. according to the source wavelet derived to estimate with big gun collection spatial position change above.According to the He of formula 16
Formula 19, when the source wavelet of estimation is not with big gun collection spatial position change (namely normal higher-order spectra), have：
When the source wavelet of estimation is with big gun collection spatial position change (namely becoming higher-order spectra), have：
Wherein, S (j,j) element on source wavelet diagonal of a matrix is represented, j represents to compile with big gun collection locus or big gun collection
Number change, vec () represent a matrix being changed into column vector, and H represents complex conjugate transposition, subscriptRepresentative takes a frequency
The jth row of section matrix, subscriptRepresent all row for taking a frequency slice matrix.When carrying out source wavelet estimation, this
Invention does not need any assumed condition, namely need not do any hypothesis to the phase of wavelet.In addition, do not examined compared to tradition
Consider the higher-order spectra problem of more subwaves, more subwave items when the present invention carries out higher-order spectra help to overcome traditional higher-order spectra
Existing (amplitude and phase) nonuniqueness problem.
Fig. 2 is the flow chart of a preferred embodiment of the method for estimation of the present invention, and the method for estimation comprises the following steps：
A) following data are obtained：Data d, threshold value μ, gestGreatest iteration update times kmax, linear Bregman algorithms
Inner iterative number lmax, sestGreatest iteration update times mmax。
B) initialize:Iteration variable k=0, m=0, b=d, gestAnd sestIt is initialized as 0.
C) k ＜ k are judgedmaxWhether set up, if so, step c-i is then repeated, otherwise performs step 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) m ＜ m are judgedmaxWhether set up, if so, then perform step g-i.
G) according to the g after b and renewalest, call the formula of formula 20 or 21 to update to obtain sest.Calculated by formula 20 or 21
Go out each diagonal matrixIn diagonal on each element, it is all according to what is calculatedUpdate sest。
H) m=m is judgedmaxWhether -1 set up, by k if setting upmax- k is assigned to lmax。
I) m+1 is assigned to m.
J) according to the result after final updated, the estimated result of source wavelet and surface-related multiple is drawn.
There are two unknown numbers in said process：The section s of source waveletestWith the section g of Green's functionest, it is anti-solving
Traditional blind deconvolution method is similar in problem flow, 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, often
Need iteration several times because according to sparse Renew theory, generally require iteration could restrain to obtain several times one it is optimal sparse
Solution.Obtaining a sparse solution gestAfterwards, then source wavelet estimation is carried out.Estimate obtained sestIt is immediately used for inverting renewal
gest.Whole iteration major cycle is meeting that some stopping criterion backed off after random terminates.Alternately renewal gestWith sestPreset times determine
Due to kmaxAnd mmax, work as gestReach k after maximum update timesmaxAnd sestWhen being not reaching to maximum update times, performing
gestLast time renewal after, perform a sestRenewal i.e. complete gestWith sestAlternating renewal, wherein every time carry out
The l of linear Bregman algorithmsmaxIt is equal；Work as sestReach k after maximum update timesmaxAnd gestIt is not reaching to maximum update times
When, performing sestLast time renewal after, only carry out the linear Bregman algorithms of last time, last time performs linear
The l of Bregman algorithmsmaxIt is assigned gestMaximum update times and the difference of updated number, and before each
Perform the l of linear Bregman algorithmsmaxIt is equal.The l of linear Bregman algorithms is wherein carried out every timemaxEqual last time is more
New gestThat is g estimated result, last time update sestAs s estimated result, it is following to use gestRepresent g estimation knot
Fruit, sestRepresent s estimated result.
Obtaining gestAnd sestAfterwards, calculate and obtain 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
The section for the surface-related multiple time domain for obtaining different orders is calculated, such as, single order surface-related multiple m1Calculation formula be：
m1=rsurfgest* p=rsurfgest*gest*sest (22)
The section m of the time domain of the more subwaves of second-order surface2Calculation formula be：
m2=rsurfgest*m1=(rsurf)2gest*gest*gest*sest (23)
The calculation formula of the surface-related multiple of other orders is by that analogy.
Computationally intensive part when being estimated due to more subwaves is consumed is doing matrix-matrix product in each frequency slice, namely
8 formulas are calculated, so being taken to further reduce the computing of more subwave estimations, the present invention limits characteristic according to the band of geological data
(namely the frequency content of geological data is concentrated mainly in the frequency band range of a certain finite length), the present invention is in more subwaves, shake
Merely with the dominant frequency band that geological data confidence level and signal to noise ratio are high during the higher-order spectra of source, the calculating of method is greatly reduced
Time-consuming, this is different from conventional method and make use of all (including confidence level and signal to noise ratio low) frequency contents.
After primary wave or surface-related multiple is obtained, it is possible to draw subsurface structure figure respectively using them.
In A Salt-dome Model data example
The data of the more subwave research group simulation generations of DELPHI of TU Delft Polytechnics are shown in Fig. 4 a, the data
It is input by rate pattern and density model, using finite-difference forward modeling, time sampling interval is 4ms.During forward simulation
Source wavelet used is Ricker wavelets, the same per big gun wavelet, namely Chang Zibo.The shallow-layer of model is the water of about 200 meters of depths
Layer.Contrasted for same Lin and Herrmann REPSI methods, total iterations is set in processing procedure as 71.
The result of REPSI methods is provided by Lin and Herrmann.The present invention frequency band range used in processing is (former for 0-70Hz
The whole frequency band range that begins is 0-125Hz).
In from Fig. 3 to Fig. 4 from the point of view of result of the present invention, the present invention has obtained preferable effect.In original input data
In Fig. 4 a, more same subwaves of subwave are mixed in together, interfere, differentiate it is unclear.After the inventive method is handled, more subwaves
It is accurately separated with primary wave.The weak lineups for belonging to tomography in Fig. 4 a initial data are difficult identification, at the present invention
After reason, it is able to clearly identify in Fig. 4 c, shows.Due to it is contemplated that processing more subwaves related to surface, and then interlayer
More subwaves can be comprised in the primary wave obtained after processing, see Fig. 4 c.
By contrasting result (Fig. 4 c) of the present invention, Lin and Herrmann REPSI methods and resultses (Fig. 4 d), it is seen that,
The primary wave that present invention processing obtains is more accurate, and more subwaves are able to preferably estimation, compacting；And Lin and Herrmann
Some more subwave lineups still can be seen in REPSI methods and resultses (Fig. 4 d), as shown in Fig. 4 d arrows.
Result (Fig. 5 of source wavelet, Lin and Herrmann that the present invention estimates REPSI methods estimation in comparison diagram 5a
(a) what solid black lines represented is the source wavelet that the present invention estimates in, and what dotted line represented is Lin and Herrmann REPSI side
The source wavelet of method estimation), it can be seen that both results are consistent, and this demonstrates the correct of the inventive method 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
Caused result is more consistent, sees Fig. 5 b.
The data examples of Pluto 1.5
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 modelings in the world and discloses release, is more for open test in the world
One of signal data of subwave algorithm.Source wavelet used is Ricker wavelets during the data forward simulation, per big gun wavelet
Equally, namely Chang Zibo.Time sampling interval is 8 milliseconds.We directly make use of the original of the open release of SMAART JV tissues
Beginning data, in addition to direct wave is cut off, any pretreatment is not done to initial data.From being originally inputted in Fig. 6 a and Fig. 7 a
Substantial amounts of multiple wave energy from the point of view of data, in the data be present.Frequency band range used is 0- in processing procedure of the present invention
50Hz。
From Fig. 6 and Fig. 7 from the point of view of result, in original input data Fig. 6 a and Fig. 7 a, more same subwaves of subwave mix
It is miscellaneous together, interfere, differentiate it is unclear, and the present invention realize accurate primary wave, more subwaves estimation with separating (especially
It is the signified part of arrow in Fig. 7).Source wavelet used is not change with big gun collection when being generated due to this digital simulation, so
It is more consistent that the present invention carries out result caused by change higher-order spectra, sees Fig. 8.
Fig. 9 illustrates another scheme for evaluating more subwave estimation technique validity.Essence caused by Fig. 9 is that feedback changes
For model, i.e.,：The primary wave (Fig. 9 b) that the present invention estimates gets off to predict that single order surface is related repeatedly incomingly as secondary focus
Ripple (Fig. 9 c), single order surface-related multiple (Fig. 9 c) get off to predict that second-order surface is related repeatedly incomingly as secondary focus
Ripple (Fig. 9 d), second-order surface related multiple (Fig. 9 d) get off three rank surfaces of prediction correlation repeatedly incomingly as secondary focus
Ripple (Fig. 9 e), three rank surface-related multiples (Fig. 9 e) get off to predict that quadravalence surface is related repeatedly incomingly as secondary focus
Ripple (Fig. 9 f), quadravalence surface-related multiple (Fig. 9 f) get off five rank surfaces of prediction correlation repeatedly incomingly as secondary focus
Ripple (Fig. 9 g), etc. is by that analogy.For this data, surface-related multiple (Fig. 9 c to figure for the different orders predicted
It can 9g) be identified in original input data (Fig. 9 a) by identification easily.Primary wave estimated by the present invention, not same order
The same original input data (Fig. 9 a) of total data (Fig. 9 h) caused by the cumulative summation of secondary surface-related multiple is almost completely
Unanimously, difference (Fig. 9 i) very little between the two.
Figure 10 further illustrates the result that the present invention is applied to other big gun collection in the data of Pluto 1.5, can see
Go out primary wave and more subwaves are able to more accurately estimate with separating.
Embodiments of the invention are described above in conjunction with accompanying drawing, but the invention is not limited in above-mentioned specific
Embodiment, above-mentioned embodiment is only schematical, rather than restricted, the ordinary skill people of this area
Member in the case of present inventive concept and scope of the claimed protection is not departed from, can also make very under the enlightenment of the present invention
Multi-form, these are belonged within the protection of the present invention.
Claims (8)
1. a kind of surface-related multiple and higher-order spectra method based on linear Bregman algorithms, it is characterised in that include following steps
Suddenly：
S1, using following step S11 and S12 alternately update Green's function to be estimated with source wavelet to be estimated until
Reach preset times, wherein the source wavelet to be estimated for carrying out step S11 first has default initial value：
S11, the input data according to source wavelet to be estimated and collection, updated using linear Bregman algorithms to be estimated
Green's function；
S12, according to the Green's function to be estimated after renewal 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 the result after final updated, draw the estimated result of at least one of source wavelet and surface-related multiple；
Wherein, subscript est is represented to be estimated,The element on the diagonal of source wavelet matrix section is represented, j is represented with big gun
Collect locus or big gun collection numbering, vec () represents a matrix being changed into column vector, subscriptRepresentative takes a frequency slice
The jth row of matrix, subscriptAll row for taking a frequency slice matrix are represented,For the matrix table of the frequency domain of input data
Up to the section of formula,For the section 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 section of the matrix expression of the frequency domain of Green's function.
2. surface-related multiple according to claim 1 and higher-order spectra method, it is characterised in that according to focus to be estimated
Wavelet and the input data of collection, updating Green's function to be estimated using linear Bregman algorithms specifically includes：
According to formula yl+1=yl-tlAHΠσ(Axl- b) and xl+1=shrink (yl+1, μ) it is iterated, y is obtained successively1、x1、
y2、x2、y3…、Wherein,For the time of the renewal result of the Green's function to be estimated after this execution step S11
The vectorial expression-form in domain, lmaxFor this execution step S11 when linear Bregman algorithms maximum iteration；
In formula, x0With y0The vectorial expression-form of the time-domain of the default initial value of Green's function to be estimated is equal to, is observed
Data vector b be the input data time-domain vectorial expression-form, tlWhat is represented is step-length, and its expression formula is：
<mrow>
<msub>
<mi>t</mi>
<mi>l</mi>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
</mrow>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msup>
<mi>A</mi>
<mi>H</mi>
</msup>
<msub>
<mi>&Pi;</mi>
<mi>&sigma;</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
</mrow>
</mfrac>
<mo>,</mo>
</mrow>
Mapping function Πσ(Axl- b) definition be：
<mrow>
<msub>
<mi>&Pi;</mi>
<mi>&sigma;</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<mi>&sigma;</mi>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>2</mn>
</msub>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
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 takes maximum, | | function representation takes absolute value, sign (x) tables
Show symbol value function；
Matrix operator A represents to carry out the operational set of operations described below：Vector table shape by Green's function to be estimated by time-domain
Formula carries out the vectorial sheet form that Fourier transformation is changed into frequency domain, the vectorial sheet form by Green's function to be estimated by frequency domain
It is changed into the matrix sheet form of frequency domain, according to formulaCalculate corresponding to input data and simulate respectively
Each section of the expression matrix form of the frequency domain of dataAccording to each sectionInversefouriertransform is carried out to obtain
The vectorial expression-form of the time-domain of analogue data corresponding to input data.
3. surface-related multiple according to claim 1 and higher-order spectra method, it is characterised in that in the step S2, shake
The estimated result of source wavelet obtains 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 table of the time-domain of surface-related multiple reaches shape
Formula, g are the vectorial expression-form of the time-domain of Green's function, and d is the vectorial expression-form of the time-domain of input data.
4. surface-related multiple according to claim 2 and higher-order spectra method, it is characterised in that if Green's letter to be estimated
Number has been maxed out 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 updated number, and before
L during step S11 is performed every timemaxIt is equal；If source wavelet to be estimated has been maxed out update times and lattice to be estimated
Woods function is not up to maximum update times, then follow-up only to carry out a step S12, and l during each execution step S11maxPhase
Deng.
5. a kind of surface-related multiple and higher-order spectra system based on linear Bregman algorithms, it is characterised in that include：
Alternately updating block, for being treated using following Green's function updating blocks and the alternating renewal of source wavelet updating block
The Green's function of estimation is with source wavelet to be estimated until reaching preset times, wherein calling Green's function updating block first
Source wavelet to be estimated there is default initial value：
Green's function updating block, according to source wavelet to be estimated and the input data of collection, calculated using linear Bregman
Method updates Green's function to be estimated；
Source wavelet updating block, according to the Green's function to be estimated after renewal 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 the result after final updated, drawing source wavelet to be estimated and to be estimated
At least one of surface-related multiple estimated result；
Wherein, subscript est is represented to be estimated,The element on the diagonal of source wavelet matrix section is represented, j is represented with big gun
Collect locus or big gun collection numbering, vec () represents a matrix being changed into column vector, subscriptRepresentative takes a frequency slice
The jth row of matrix, subscriptAll row for taking a frequency slice matrix are represented,For the matrix table of the frequency domain of input data
Up to the section of formula,For the section 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 section of the matrix expression of the frequency domain of Green's function.
6. surface-related multiple according to claim 5 and higher-order spectra system, it is characterised in that according to focus to be estimated
Wavelet and the input data of collection, updating Green's function to be estimated using linear Bregman algorithms specifically includes：
According to formula yl+1=yl-tlAHΠσ(Axl- b) and xl+1=shrink (yl+1, μ) it is iterated, y is obtained successively1、x1、
y2、x2、y3…、Wherein,For the time of the renewal result of the Green's function to be estimated after this execution step S11
The vectorial expression-form in domain, lmaxFor this execution step S11 when linear Bregman algorithms maximum iteration；
In formula, x0With y0The vectorial expression-form of the time-domain of the default initial value of Green's function to be estimated is equal to, is observed
Data vector b be the input data time-domain vectorial expression-form, tlWhat is represented is step-length, and its expression formula is：
<mrow>
<msub>
<mi>t</mi>
<mi>l</mi>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
</mrow>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msup>
<mi>A</mi>
<mi>H</mi>
</msup>
<msub>
<mi>&Pi;</mi>
<mi>&sigma;</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
</mrow>
</mfrac>
<mo>,</mo>
</mrow>
Mapping function Πσ(Axl- b) definition be：
<mrow>
<msub>
<mi>&Pi;</mi>
<mi>&sigma;</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<mi>&sigma;</mi>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>2</mn>
</msub>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>Ax</mi>
<mi>l</mi>
</msub>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
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 takes maximum, | | function representation takes absolute value, sign (x) tables
Show symbol value function；
Matrix operator A represents to carry out the operational set of operations described below：Vector table shape by Green's function to be estimated by time-domain
Formula carries out the vectorial sheet form that Fourier transformation is changed into frequency domain, the vectorial sheet form by Green's function to be estimated by frequency domain
It is changed into the matrix sheet form of frequency domain, according to formulaCalculate corresponding to input data and simulate respectively
Each section of the expression matrix form of the frequency domain of dataAccording to each sectionInversefouriertransform is carried out to obtain
To the vectorial expression-form of the time-domain of analogue data corresponding to input data.
7. surface-related multiple according to claim 5 and higher-order spectra system, it is characterised in that the estimated result obtains
In unit,
The estimated result of source wavelet obtains 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 table of the time-domain of surface-related multiple reaches shape
Formula, g are the vectorial expression-form of the time-domain of Green's function, and d is the vectorial expression-form of the time-domain of input data.
8. surface-related multiple according to claim 6 and higher-order spectra system, it is characterised in that if Green's letter to be estimated
Number has been maxed out 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 updated number, and before
L during step S11 is performed every timemaxIt is equal；If source wavelet to be estimated has been maxed out update times and lattice to be estimated
Woods function is not up to maximum update times, then follow-up only to carry out a step S12, and l during each execution step S11maxPhase
Deng.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710524014.7A CN107390261B (en) | 2017-06-27 | 2017-06-27 | Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710524014.7A CN107390261B (en) | 2017-06-27 | 2017-06-27 | Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithm |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107390261A true CN107390261A (en) | 2017-11-24 |
| CN107390261B CN107390261B (en) | 2019-02-12 |
Family
ID=60334679
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710524014.7A Expired - Fee Related CN107390261B (en) | 2017-06-27 | 2017-06-27 | Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107390261B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111045084A (en) * | 2020-01-06 | 2020-04-21 | 中国石油化工股份有限公司 | Multi-wave self-adaptive subtraction method based on prediction feature extraction |
| CN111308556A (en) * | 2020-03-17 | 2020-06-19 | 清华大学 | Fast robust curvelet domain multiple subtraction technology based on frequency division constraint |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120051179A1 (en) * | 2010-08-26 | 2012-03-01 | Snu R&Db Foundation | Method and apparatus for time-domain reverse-time migration with source estimation |
| EP2755053A2 (en) * | 2013-01-11 | 2014-07-16 | CGG Services SA | A system and method for the removal of shallow water multiples using a hybrid multi-channel prediction method |
| CN104199089A (en) * | 2014-08-22 | 2014-12-10 | 电子科技大学 | AVO inversion method based on information geometry |
| CN105259575A (en) * | 2015-10-12 | 2016-01-20 | 中国石油大学(华东) | Method for fast predicting 3D surface-related multiples |
| CN105334537A (en) * | 2015-10-26 | 2016-02-17 | 中国石油大学(华东) | Primary wave and multiple wave separation method based on alternative splitting Bregman iterative algorithm |
-
2017
- 2017-06-27 CN CN201710524014.7A patent/CN107390261B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120051179A1 (en) * | 2010-08-26 | 2012-03-01 | Snu R&Db Foundation | Method and apparatus for time-domain reverse-time migration with source estimation |
| EP2755053A2 (en) * | 2013-01-11 | 2014-07-16 | CGG Services SA | A system and method for the removal of shallow water multiples using a hybrid multi-channel prediction method |
| CN104199089A (en) * | 2014-08-22 | 2014-12-10 | 电子科技大学 | AVO inversion method based on information geometry |
| CN105259575A (en) * | 2015-10-12 | 2016-01-20 | 中国石油大学(华东) | Method for fast predicting 3D surface-related multiples |
| 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 (3)
| Title |
|---|
| XINTAO CHAI 等: ""A linearized Bregman method for compressive waveform inversion"", 《SEG INTERNATIONAL EXPOSITION AND 86TH ANNUAL MEETING》 * |
| 张慧 等: ""A-线性Bregman迭代算法"", 《计算数学》 * |
| 彭才 等: ""动态地震子波估计"", 《大庆石油地质与开发》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111045084A (en) * | 2020-01-06 | 2020-04-21 | 中国石油化工股份有限公司 | Multi-wave self-adaptive subtraction method based on prediction feature extraction |
| CN111308556A (en) * | 2020-03-17 | 2020-06-19 | 清华大学 | Fast robust curvelet domain multiple subtraction technology based on frequency division constraint |
| CN111308556B (en) * | 2020-03-17 | 2021-04-13 | 清华大学 | Fast robust curvelet domain multiple subtraction method based on frequency division constraint |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107390261B (en) | 2019-02-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gholami | Sparse time–frequency decomposition and some applications | |
| Pinnegar et al. | Application of the S transform to prestack noise attenuation filtering | |
| Liu et al. | Seismic data interpolation beyond aliasing using regularized nonstationary autoregression | |
| Lark et al. | Changes in variance and correlation of soil properties with scale and location: analysis using an adapted maximal overlap discrete wavelet transform | |
| CN107367760B (en) | Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm | |
| Zhou et al. | Robust noise attenuation based on nuclear norm minimization and a trace prediction strategy | |
| Wang et al. | Self-training and learning the waveform features of microseismic data using an adaptive dictionary | |
| Zhou et al. | Sparse dictionary learning for seismic noise attenuation using a fast orthogonal matching pursuit algorithm | |
| Yang et al. | Denoising of distributed acoustic sensing data using supervised deep learning | |
| de Oliveira Lyrio et al. | Efficient automatic denoising of gravity gradiometry data | |
| Nascimento et al. | High-resolution acoustic impedance inversion to characterize turbidites at Marlim Field, Campos Basin, Brazil | |
| US20220342103A1 (en) | Noise Attenuation Methods Applied During Simultaneous Source Deblending and Separation | |
| Saad et al. | Uncovering the microseismic signals from noisy data for high-fidelity 3D source-location imaging using deep learning | |
| CN107390261B (en) | Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithm | |
| 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 | |
| CN114091538B (en) | Intelligent noise reduction method for discrimination loss convolutional neural network based on signal characteristics | |
| Oboué et al. | Erratic and random noise attenuation using adaptive local orthogonalization | |
| Gentilhomme et al. | Ensemble-based multi-scale history-matching using second-generation wavelet transform | |
| Xie et al. | Using wavelet‐domain adaptive filtering to improve signal‐to‐noise ratio of nuclear magnetic resonance log data from tight gas sands | |
| Heidary et al. | Wavelet analysis in determination of reservoir fluid contacts | |
| Yuan et al. | Ground roll attenuation based on an empirical curvelet transform | |
| Huang et al. | Shannon entropy-based seismic local correlation measure and enhancement | |
| CN110008633A (en) | Road noise drawing method and system based on artificial intelligence deep neural network | |
| Huang et al. | Data-assimilated time-lapse visco-acoustic full-waveform inversion: Theory and application for injected CO2 plume monitoring |
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 |
Granted publication date: 20190212 Termination date: 20200627 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |