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
- msub
- mrow
- function
- 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 80
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 66
- 238000001228 spectrum Methods 0.000 title abstract description 8
- 230000008859 change Effects 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 73
- 239000013598 vector Substances 0.000 claims description 64
- 230000014509 gene expression Effects 0.000 claims description 60
- 238000004613 tight binding model Methods 0.000 claims description 52
- 230000009466 transformation Effects 0.000 claims description 10
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 8
- 238000013507 mapping Methods 0.000 claims description 7
- 230000004044 response Effects 0.000 abstract description 19
- 230000003595 spectral effect Effects 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 28
- 238000010586 diagram Methods 0.000 description 24
- 230000008569 process Effects 0.000 description 10
- 238000012545 processing Methods 0.000 description 10
- 238000004364 calculation method Methods 0.000 description 9
- 238000001914 filtration Methods 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 230000001629 suppression Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000008520 organization Effects 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 239000003345 natural gas Substances 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 239000003209 petroleum derivative Substances 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 239000013535 sea water Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 241000143252 Idaea infirmaria Species 0.000 description 1
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 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
- 230000006978 adaptation Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 239000003054 catalyst Substances 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000001186 cumulative effect Effects 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
- 238000013461 design Methods 0.000 description 1
- 230000000694 effects Effects 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
- 238000005755 formation reaction Methods 0.000 description 1
- 235000021384 green leafy vegetables Nutrition 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000003825 pressing Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 150000003839 salts Chemical class 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000002945 steepest descent method Methods 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
The invention relates to the field of seismic exploration, in particular to the aspect of seismic signal processing in the seismic exploration technology, and more particularly relates to a surface multiple and wavelet estimation method and system based on a linear Bregman algorithm.
Background
Multiple reflection (multiple) related problems are one of the most prominent problems common in reflection seismic exploration. The surface-related multiple problem is particularly serious in marine seismic exploration because the seawater surface is a strong reflection interface. At present, exploration and development of offshore petroleum and natural gas combustible ice resources become an urgent priority of various large petroleum energy companies at home and abroad. The petroleum and natural gas combustible ice in south China sea is rich in resources and is a hotspot and focus of current marine exploration. In current oil industry production, one of the basic assumptions of conventional seismic data processing techniques is: the input reflection data volume consists of only primary reflected waves (primaries). Multiples can be mistakenly considered as primaries or mixed in part of primaries. The multiple wave causes strong interference to the processing and imaging of the primary wave, and influences resource exploration, so the multiple wave estimation technology has excellent practical industrial application prospect.
The multiple estimation or suppression methods currently used in the industry are mainly classified into the following categories: one class of filtering methods is called as filtering methods, which use the characteristic difference between multiples and primaries to suppress the multiples; firstly, predicting multiple waves from seismic data, and then adaptively subtracting the multiple waves from original data, namely prediction subtraction; there is also an inverse scattering series method to suppress multiples.
The filtering method assumes that the multiples and the primaries have obvious characteristic difference in a specific transform domain, and various mathematical algorithms are designed to acquire the difference and eliminate the multiples. The filtering method utilizes the main characteristic differences including the periodicity and separability of multiple waves. When the filtering method is applied to multiple attenuation of some complex structures, the condition that the characteristic difference between the multiple and the primary wave is small or no difference is met can be met, and the multiple suppression effect is not obvious.
There are also problems with predicting multiples using the backscattering order: the calculated amount is large and the cost is high; it is difficult to predict multiples at far offsets; the prediction capability of multiples in complex media also has certain limitations. Among other things, the high computational cost makes this method difficult to apply to the processing of real data.
Berkhout proposed a feedback model theory for describing complex multiples in 1982, and laid the mathematical and physical basis of a feedback iteration multiple pressing method, but known source wavelets are needed when multiple waves are eliminated. To illustrate this problem and to facilitate the following description, we first introduce the feedback model equations that are closely related to the present invention, i.e., the equation
Wherein D is a matrix expression form of a time domain of the acquired input data,is a slice of a matrix expression of a frequency domain after Fourier transform, S is a matrix expression of a time domain of a seismic source wavelet to be estimated,is a slice of a matrix expression of S in a frequency domain after Fourier transform, G is a matrix expression of a time domain of a Green function (i.e. a subsurface impulse response without multiples),for a slice of the matrix expression of the G Fourier transformed frequency domain we use ^ symbols to represent the data or variables of the frequency domain, rsurfFor the reflection coefficient of the air-water interface, let r besurfIs 1 ═ 1, whereinAndrepresents a slice of the frequency of the estimated primary wave (i.e. a time domain convolution),that is to sayWhereinIs a slice of the frequency domain of the primary wave. The term "primary" as used herein refers to all waves except surface multiples, i.e., where the primary contains interbed multiples (which are not related to the surface).One slice, i.e. data, representing the frequency domain of the estimated surface multiplesAs a quadratic source down-propagation with the underground impulse response green functionIs used to predict the slice of the surface-related multiples frequency domain. From equation 1, it can be seen that accurate surface multiple estimation depends on the accuracy of the estimated source wavelet.
Verschuuur et al in 1992 proposed a surface-related multiple suppression technique SRME that estimates the source wavelet from actual data by introducing adaptive subtraction based on the assumption that the primary energy is minimal. To avoid the non-linearity problem created by solving the source wavelet, Berkhout and Verschuur proposed an iterative SRME method in 1997. The precondition that the primary wave energy is minimal is not fulfilled in all cases. In addition, in practical application, SRME also has the problem that adaptive prediction phase subtraction damages the effective signal.
In order to improve the problem of the SRME that the prediction subtraction damages the effective signal, van groenerstijn and Verschuur proposed a sparse inversion primary wave estimation technique, namely EPSI, in 2009. To improve the problems in EPSI, Lin and Herrmann proposed a robust sparse inversion primary estimation technique, REPSI, in 2013. EPSI and REPSI are prior art related to the present invention, and their technical solutions and problems and disadvantages will be described in detail below.
In addition, when the actual field seismic data is collected, the source wavelets excited by each shot are different, and all the existing multiple estimation methods assume that all the source wavelets are the same, so that the multiple estimation precision is influenced, and the source wavelets which change along with the shot gather space position can be estimated.
Due to the fact that the number of computing resources (memory and processor) consumed by each iteration is large and the computing time is long during multiple estimation, the efficiency of the method is improved, the convergence rate of the algorithm is improved, the complexity of the algorithm is simplified, the adjusting parameters of the method are reduced besides the estimation accuracy of multiple and source wavelets, and the method is also the development trend of the multiple estimation technology.
The prior art related to the present invention is described as follows.
Technical scheme of prior art I
van Groenetijn and Verschuur put forward a sparse inversion primary wave estimation method EPSI on the basis of SRME theory in 2009, in the iterative inversion process, fit data residual is used as driving, two unknown parameters of an underground impulse response G and a seismic source wavelet S are alternately updated through a steepest descent method, and then the primary wave is reconstructedAnd its corresponding surface-related multiplesCompared with SRME, EPSI avoids the step of self-adaptive subtraction and better protects the primary effective signal.
When the EPSI method of van Groenetijn and Verschuur carries out the estimation of the source wavelet, all the source wavelets are assumed to be the same, and an estimation scheme (namely, ordinary wavelet estimation) of the source wavelet changing along with the shot gather space position is not provided. The method can estimate the seismic source wavelet changing with the shot gather space position.
In order to ensure the sparsity of the underground impulse response G in the time domain, in the process of updating the impulse response G, the EPSI updates the impulse according to the sequence of time from small to large and amplitude from strong to weak. That is, in order to promote the time-domain sparsity of the subsurface impulse response G, the EPSI method designs a number of artificial adjustment parameters, such as the EPSI needs to determine: a time window for determining the update range of the green function (i.e. performing an update within the time window and setting zero outside the time window); the number of reflection event axes in the green function at each iteration, etc. As indicated in the 2013 paper by Lin and Herrmann: even if the number of subsurface reflecting interfaces is known precisely (which is not possible in practice), the EPSI method does not make full use of this a priori information. The underground impulse response updating of the EPSI is easily influenced by parameters such as an impulse selecting window, a pickup quantity and the like, and further influences the multiple estimation precision. In essence, the sparsity of the green function in the EPSI method is controlled by human, rather than by a mathematical optimization algorithm. This drawback limits the applicability of the EPSI process in the actual industry.
Prior art relating to the invention
Based on the same theory of EPSI (namely feedback model 1 formula), Lin and Herrmann modify sparse constraint conditions of EPSI underground impulse response into one-norm constraint in 2013, and adopt a spectrum mapping gradient one-norm algorithmThe sparse underground impulse response Green function is obtained by solving, and the stability of the original EPSI method is improved to a certain extent. The method proposed by Lin and Herrmann in 2013 was named Robust EPSI (REPSI). The REPSI method proposed by Lin and Herrmann, like EPSI, alternately updates two unknown parameters of the underground impulse response G and the seismic source wavelet S, and then reconstructs the primary waveAnd its corresponding surface-related multiplesThe biggest difference between EPSI and REPSI is that: the sparsity of Green function in EPSI is mostly controlled by artificial regulation parameters, and REPSI adoptsAnd solving to obtain a sparse underground impulse response Green function by a mathematical optimization algorithm.
The REPSI method of Lin and Herrmann also assumes that all source wavelets are the same when estimating the source wavelet, and does not provide an estimation scheme (i.e., constant wave estimation) for the source wavelet that changes with the shot gather spatial position. The method can estimate the seismic source wavelet changing with the shot gather space position. In addition, when the REPSI carries out the ordinary wave estimation, the LSQR algorithm is adopted to carry out iterative solution to obtain the seismic source wavelet, thereby greatly increasing the calculated amount of the method and reducing the efficiency. The method can estimate the source wavelet in one step.
Another important disadvantage of the REPSI method by Lin and Herrmann is that: it relies on a very complex spectral mapping gradient-norm algorithmTo solve to obtain sparse underground impulse response green's function.Is an open source MATLAB program that can be downloaded by accessing https:// github. com/mpf/spgl 1.The complexity of the algorithm is mainly reflected in that: numerous tuning parameters, numerous subfunctions, and thousands of lines of code. Adaptation was attempted by Schlumberger oilfield servicesIt is applied in the actual industry, but is hinderedThe complexity of the algorithm has not been successful in the end.The complexity of the algorithm also reduces the efficiency of the REPSI method and hinders the practical application of the REPSI method in industry.
Disclosure of Invention
The technical problem to be solved by the present invention is to provide a method and a system for estimating surface multiples and wavelets based on a linear Bregman algorithm, aiming at the technical defects of the above-mentioned EPSI and REPSI in the prior art.
According to one aspect of the present invention, to solve the technical problems, the present invention provides a surface multiple and wavelet estimation method based on a linear Bregman algorithm, comprising the following steps:
s1, alternately updating the green function to be estimated and the source wavelet to be estimated using the following steps S11 and S12 until reaching a preset number of times, wherein the source wavelet to be estimated, which is first performed in step S11, has a preset initial value:
s11, updating the Green function to be estimated by adopting a linear Bregman algorithm according to the seismic source wavelet to be estimated and the acquired input data;
s12, updating the seismic source wavelet to be estimated by adopting the following formula according to the updated Green function to be estimated and the input data;
when the source wavelet to be estimated does not change along with the shot gather space position:
when the source wavelet to be estimated changes along with the shot gather space position:
s2, obtaining at least one estimation result of the seismic source wavelet and the surface multiple according to the final updated result;
wherein the superscript est represents the value to be estimated,representing elements on the diagonal of a slice of the source wavelet matrix, j representing the shot gather spatial position or shot gather number, vec () representing the transformation of a matrix into a column vector, subscriptRepresenting column j, subscript, taking a frequency slice matrixRepresenting taking all columns of a frequency slice matrix,for a slice of the matrix expression of the frequency domain of the input data,being a slice of a matrix expression of the frequency domain of the seismic source sub-waves, rsurfIs the reflection coefficient of the contact surface of air and water,a slice of a matrix expression of the frequency domain of the green's function.
Further, in the method for estimating surface multiples and wavelets based on the linear Bregman algorithm of the present invention, the specific steps of updating the green function to be estimated by using the linear Bregman algorithm according to the source wavelet to be estimated and the acquired input data include:
according to the formula yl+1＝yl-tlAHΠσ(Axl-b) and xl+1＝shrink(yl+1μ) are iterated to find y in turn1、 x1、y2、x2、y3…、Wherein,a vector expression form l in the time domain for the update result of the green function to be estimated after the current execution of step S11maxThe maximum iteration number of the linear Bregman algorithm when the step S11 is executed this time;
in the formula, x0And y0A vector representation form of the time domain which is equal to the preset initial value of the green function to be estimated, an observation data vector b is the vector representation form of the time domain of the input data, tlWhat is represented is the step size, which is expressed as:
mapping function Пσ(Axl-b) is defined as:
the soft threshold function shrnk (x, μ) is defined as:
shrink(x,μ)＝max(|x|-μ,0)sign(x)，
wherein, σ represents a noise level factor, the max () function represents taking a maximum value, the | | function represents taking an absolute value, and sign (x) represents a symbol value taking function;
the matrix operator a represents a set of operations that perform the following operations: fourier transformation is carried out on the Green's function to be estimated from a vector table form of a time domain to a vector table form of a frequency domain, the Green's function to be estimated is changed from the vector table form of the frequency domain to a matrix table form of the frequency domain, and the Green's function to be estimated is transformed from the vector table form of the frequency domain to a matrix table form of the frequency domain according to a formulaCalculating each slice of the matrix expression form of the frequency domain of the analog data corresponding to the input dataAccording to each sliceAnd performing inverse Fourier transform to obtain a vector expression form of a time domain of the analog data corresponding to the input data.
Further, in the surface multiples and wavelet estimation method based on the linear Bregman algorithm of the present invention, in step S2,
the estimation result of the seismic source wavelet is obtained according to the result of updating the seismic source wavelet to be estimated for the last time;
the estimation result of the surface multiple is based on m ═ rsurfg x d, m is a vector expression form of a time domain of the surface multiple, g is a vector expression form of a time domain of the green function, and d is a vector expression form of a time domain of the input data.
Further, in the surface multiples and wavelet estimation method based on the linear Bregman algorithm of the present invention, if the green function to be estimated has reached the maximum update times and the source wavelet to be estimated has not reached the maximum update times, then when step S11 is performed again, lmaxUpdated to the maximum number of updates of the source wavelet to be estimated minus the number of updates already, and l is previously set each time step S11 is executedmaxEqual; if it is to be estimatedIf the source wavelet has reached the maximum update times and the greens function to be estimated has not reached the maximum update times, then step S12 is performed only once, and step S11 is executed each timemaxAre equal.
The invention also provides a surface multiple and wavelet estimation system based on linear Bregman algorithm for solving the technical problem, comprising:
the alternative updating unit is used for alternatively updating the Green function to be estimated and the seismic source wavelet to be estimated by using the following Green function updating unit and the seismic source wavelet updating unit until the preset times are reached, wherein the seismic source wavelet to be estimated, which calls the Green function updating unit for the first time, has a preset initial value:
the green function updating unit is used for updating the green function to be estimated by adopting a linear Bregman algorithm according to the seismic source wavelet to be estimated and the acquired input data;
the seismic source wavelet updating unit is used for updating the seismic source wavelet to be estimated by adopting the following formula according to the updated green function to be estimated and the input data;
when the source wavelet to be estimated does not change along with the shot gather space position:
when the source wavelet to be estimated changes along with the shot gather space position:
the estimation result acquisition unit is used for obtaining an estimation result of at least one of the seismic source wavelet to be estimated and the surface multiple to be estimated according to the final updated result;
wherein the superscript est represents the value to be estimated,representing elements on the diagonal of a slice of the source wavelet matrix, j representing the shot gather spatial position or shot gather number, vec () representing the transformation of a matrix into a column vector, subscriptRepresenting column j, subscript, taking a frequency slice matrixRepresenting taking all columns of a frequency slice matrix,for a slice of the matrix expression of the frequency domain of the input data,being a slice of a matrix expression of the frequency domain of the seismic source sub-waves, rsurfIs the reflection coefficient of the contact surface of air and water,a slice of a matrix expression of the frequency domain of the green's function.
Further, in the surface multiple and wavelet estimation system based on the linear Bregman algorithm of the present invention, the updating of the green function to be estimated by using the linear Bregman algorithm according to the source wavelet to be estimated and the acquired input data specifically includes:
according to the formula yl+1＝yl-tlAHПσ(Axl-b) and xl+1＝shrink(yl+1μ) are iterated to find y in turn1、 x1、y2、x2、y3…、Wherein,a vector expression form l in the time domain for the update result of the green function to be estimated after the current execution of step S11maxThe maximum iteration number of the linear Bregman algorithm when the step S11 is executed this time;
in the formula, x0And y0A vector representation form of the time domain which is equal to the preset initial value of the green function to be estimated, an observation data vector b is the vector representation form of the time domain of the input data, tlWhat is represented is the step size, which is expressed as:
mapping function Πσ(Axl-b) is defined as:
the soft threshold function shrnk (x, μ) is defined as:
shrink(x,μ)＝max(|x|-μ,0)sign(x)，
wherein, σ represents a noise level factor, the max () function represents taking a maximum value, the | | function represents taking an absolute value, and sign (x) represents a symbol value taking function;
the matrix operator a represents a set of operations that perform the following operations: fourier transformation is carried out on the Green's function to be estimated from the vector table form of the time domain to be converted into the vector table form of the frequency domain, the vector table form of the frequency domain is converted into the matrix table form of the frequency domain, and the Green's function to be estimated is converted into the matrix table form of the frequency domain according to the formulaCalculating each slice of the matrix expression form of the frequency domain of the analog data corresponding to the input dataAccording to each sliceAnd performing inverse Fourier transform to obtain a vector expression form of a time domain of the analog data corresponding to the input data.
Furthermore, in the surface multiple and wavelet estimation system based on the linear Bregman algorithm, in the estimation result acquisition unit,
the estimation result of the seismic source wavelet is obtained according to the result of updating the seismic source wavelet to be estimated for the last time;
the estimation result of the surface multiple is based on m ═ rsurfg x d, m is a vector expression form of a time domain of the surface multiple, g is a vector expression form of a time domain of the green function, and d is a vector expression form of a time domain of the input data.
Further, in the surface multiples and wavelet estimation system based on the linear Bregman algorithm of the present invention, if the green function to be estimated has reached the maximum update times and the source wavelet to be estimated has not reached the maximum update times, then when step S11 is performed again, lmaxUpdated to the maximum number of updates of the source wavelet to be estimated minus the number of updates already, and l is previously set each time step S11 is executedmaxEqual; if the source wavelet to be estimated has reached the maximum update times and the Green function to be estimated has not reached the maximum update times, then step S12 is performed only once, and l is performed each time step S11 is performedmaxAre equal.
The estimation method and the estimation system based on the linear Bregman algorithm have the following beneficial effects: the estimation method and the system can estimate the seismic source wavelet changing along with the shot gather space position, and have high efficiency, less time consumption and high accuracy of the estimated result.
Drawings
The invention will be further described with reference to the accompanying drawings and examples, in which:
FIG. 1 is a flow chart of a preferred embodiment of the linear Bregman algorithm of the present invention for solving Ax ═ b to obtain a sparse solution x;
FIG. 2 is a flow chart of a preferred embodiment of the estimation method of the present invention for estimating source wavelets and surface multiples;
FIG. 3(a) is a schematic diagram of input data containing multiples;
FIG. 3(b) is a schematic diagram of primary data estimated according to the present invention;
FIG. 3(c) is a schematic representation of estimated surface multiple data according to the present invention;
FIG. 4(a) is a schematic diagram of a zero shot section of the data of FIG. 3 (a);
FIG. 4(b) is a schematic diagram of estimated surface multiple data according to the present invention;
FIG. 4(c) is a schematic diagram of estimated primary beam data according to an embodiment of the present invention;
FIG. 4(d) is a schematic diagram of primary wave data estimated by the REPSI method of Lin and Herrmann;
FIG. 5(a) is a schematic diagram of an estimated source wavelet;
FIG. 5(b) is a schematic diagram of a frequency domain amplitude spectrum of a seismic source sub-wave as a function of shot gather estimated according to an embodiment of the present invention;
FIG. 6(a) is a graphical representation of Pluto1.5 data publicly released by the International SMAART JV organization for testing the multiple suppression method;
FIG. 6(b) is a schematic diagram of primary data estimated according to the present invention;
FIG. 6(c) is a schematic diagram of estimated surface-related multiples data according to the present invention;
FIG. 7(a) is a schematic diagram of a zero shot section of the data of FIG. 6 (a);
FIG. 7(b) is a schematic diagram of primary data estimated according to the present invention;
FIG. 7(c) is a schematic representation of estimated surface-related multiples data of the present invention;
FIG. 8(a) is a schematic diagram of source wavelet as shot gather estimated from Pluto1.5 data;
FIG. 8(b) is a schematic diagram of the frequency domain amplitude spectrum of shot gather varying source sub-waves estimated from Pluto1.5 data;
FIG. 9(a) is a data cross-sectional view of the cannon of data of FIG. 7 (a);
FIG. 9(b) is a schematic diagram of primary data estimated according to the present invention;
FIG. 9(c) is a schematic diagram of the first order surface multiple data estimated according to the present invention;
FIG. 9(d) is a schematic diagram of second order surface multiple data estimated according to the present invention;
FIG. 9(e) is a schematic diagram of third order surface multiple data estimated according to the present invention;
FIG. 9(f) is a schematic diagram of fourth order surface multiple data estimated according to the present invention;
FIG. 9(g) is a schematic diagram of fifth order surface multiple data estimated according to the present invention;
FIG. 9(h) is a schematic diagram of the cumulative summation of the data in FIGS. 9(b) to 9 (g);
FIG. 9(i) is a schematic illustration of the difference of FIGS. 9(a) and 9 (h);
FIG. 10(a) is a schematic diagram of a cross-section of a zero-shot-check office with shot numbers between 701-880 in Pluto1.5 data;
FIG. 10(b) is a schematic diagram of primary data estimated according to the present invention;
FIG. 10(c) is a schematic diagram of the estimated surface-related multiples data of the present invention.
Noun description
Primary waves (Primaries): the seismic waves are received by detectors embedded in the earth surface after being reflected once only at the underground interface in the process of stratum propagation, and are called primary reflected waves, which are called primary waves for short. The term "primary" as used herein refers to all waves except surface multiples, i.e., where the primary contains interbed multiples (which are independent of the surface).
Multiples (multiplets): the seismic waves are received by detectors embedded on the earth surface after multiple reflections occur on an underground interface or a free surface (a contact surface of air and water) in the process of propagation, and the seismic waves are called multiple reflected waves, which are called multiple reflected waves for short.
Surface-related multiples (Surface-related multiples) or Surface multiples (Surface multiples): it refers to multiple reflections associated with the free surface (the interface between air and water), and is referred to herein simply as surface multiples.
Interbed multiples (lnternal multiples, IM): refers to multiple reflections that occur between subsurface reflective interfaces (between formations).
Source wavelet: refers to the waves generated by the excitation of the seismic source.
Feedback model: and (5) feedback model.
Surface-Related Multiple Elimination (SRME): surface-related multiple suppression.
Estimation of Primaries by spark interrogation (EPSI): sparse inversion primary wave estimation.
Robust Estimation of fingerprints by spark interrogation (EPSI): robust sparse inversion primary wave estimation.
Spectral-Projected GradientThe spectrum maps a gradient norm algorithm.
LSQR: least-squares QR, Least squares QR decomposition algorithm.
Linear Bregman (LB): the linear Bregman algorithm.
Accelerated Linear Bregman (ALB): accelerated linear Bregman algorithm.
MATLAB: is commercial mathematical software produced by MathWorks company in America, and is used for high-level technical computing language and interactive environment of algorithm development, data visualization, data analysis and numerical calculation. MATLAB is a combination of two words of matrix & laboratory, meaning a matrix factory (matrix laboratory).
Signal-to-Noise Ratios (SNR): signal to noise ratio.
Detailed Description
For a more clear understanding of the technical features, objects and effects of the present invention, embodiments of the present invention will now be described in detail with reference to the accompanying drawings. For the convenience of understanding the present invention, the letters of the above and below expressions are explained as follows, and the lower case letters s, g, d, p, m are vector expressions representing the time domain, and S, G, D, P, M is a matrix expression of the time domain.A slice in the form of a matrix representation of the frequency domain of S, G, D, P, M, the superscript est representing the value to be estimated,in the form of a vector representation of the frequency domain of s, g, d, p, m, S, G, D, P, M is a three-dimensional matrix, and its slice is a two-dimensional matrix,each three-dimensional matrix may be split into multiple two-dimensional matrices. Any two-dimensional matrix Z ═ Z1Z2Z3…Zn]The vector expression form of (a) can be expressed as:
wherein Z1To ZnThe columns of the two-dimensional matrix Z are represented, whereas the matrix representation is obtained from the vector representation.
The primary, surface multiples, and source wavelet estimates utilize a primary formula that is the mathematical-physical relationship of the feedback model, i.e., formula (1). The excited seismic source wavelet s is transmitted into the ground, time domain convolution (star represents convolution) is carried out on the excited seismic source wavelet s and an underground impulse response Green function g, a primary wave p is generated, and the mathematical relational expression is as follows:
p＝g*s (2)
the generated primary wave p is reflected when meeting the contact surface of seawater and air, propagates downwards and then is transmitted into the underground as a seismic source again, and carries out time domain convolution with the underground pulse response g to generate a first-order multiple m1
m1＝g*rsurfp＝rsurfg*g*s (3)
The generated first order multiple m1When the reflection occurs on the free surface, the reflection is transmitted into the underground as a seismic source again, and time domain convolution is carried out with the underground pulse response g to generate a second-order multiple m2
m2＝g*rsurfm1＝(rsurf)2g*g*g*s (4)
Higher order surface-related multiples generate formulas and so on.
The observed total data d is the sum of the primary wave p and all the order surface-related multiples, and the mathematical expression is as follows:
and then have
The mathematical-physical relational expression of the feedback model can be summarized as follows:
expression (1) is a frequency domain matrix expression form of expression 7 in two-dimensional or three-dimensional condition, and also holds for the matrix expression form in which the parameters in expression 7 are modified into their corresponding frequency domains, and for comparison, it is rewritten here:
the main purposes of the invention are: solving for unknown variables according to expressions 7-8: the source wavelet s and the Green function g, and then the estimated primary wave p ═ g ═ s and the surface multiple m ═ rsurfg*d。
The invention first illustrates a method for inverting a sparse impulse response green's function given an estimated source wavelet.
At a given estimated source wavelet sestIn the case of (2), we write the inversion g according to expressions (7) - (8)estThe formula used when, that is:
here, theIs a simplified representative forward operator, the core of which is in the calculation formula 8In particular, matrix operatorsRepresents a set of operations that perform the following operations: fourier transformation is carried out on the Green function to be estimated from a vector table form g of a time domain to be changed into a vector table form of a frequency domainThe Green function to be estimated is represented by a vector table form of a frequency domainInto a matrix table form of the frequency domainAccording to the formulaCalculating each slice of the matrix expression form of the frequency domain of the analog data corresponding to the input dataEach slice of the frequency domain matrix expression form of the analog data corresponding to the calculated input dataPerforming inverse Fourier transform to obtain a time-domain vector expression form of analog data corresponding to the input data
To further illustrate how the present invention can be solved to obtain sparse green function g, formula 9 is rewritten as the following standard form:
Ax＝b (10)
wherein,x＝g，b＝d。
in order to seek a sparse solution of the green's function g, the present invention solves equation 10 using a linear Bregman algorithm. The linear Bregman algorithm was proposed by Yin et al in 2008, who performed a detailed analysis of the linear Bregman algorithm in 2010. Cai in 2009's paper showed that the linear Bregman algorithm aims to solve:
expressing solving for x so that the expression takes the minimum value, subject to is expressed under the defined condition.
Wherein mu is more than or equal to 0 and is a sparsity control factor, | | x | | luminance1And | | x | | non-conducting phosphor2Respectively representing vectorsNorm sumNorm, σ, represents the noise level factor in the data. Algorithm 1 below gives a detailed iteration step of the linear Bregman algorithm.
Referring to fig. 1, fig. 1 is a flowchart of the linear Bregman algorithm of the present invention for solving Ax ═ b to obtain a sparse solution x, that is, a flowchart of the algorithm 1, where the algorithm 1 specifically includes the following steps:
a) obtaining the following ginsengNumber: positive operator A, observation data vector b, initial value vector x0Threshold μ, maximum number of iterations lmax。
b) Initialization, iteration variable l is 0, x0Assign y to0
c) Repeating the following steps d to fmaxNext time
d) Calculating yl+1＝yl-tlAHПσ(Axl-b)
e) Calculating xl+1＝shrink(yl+1,μ)
f) Update l to l +1
g) To obtain a sparse solutionThe sparse solution x is the update g of the linear Bregman algorithm performed at this timeestAnd (5) obtaining the final result.
In Algorithm 1 above, t of step dlWhat is represented is the step size, which is expressed as:
wherein the superscript H represents the complex conjugate transpose. Mapping function Πσ(Axl-b) is to take into account the noise present in the data, which is defined as:
and the soft threshold function shrink (x, μ) is defined as:
shrink(x,μ)＝max(|x|-μ,0)sign(x) (14)
wherein, max (·) function represents taking maximum value, | · | function represents taking absolute value, sign (x) represents symbol value function.
The invention will be described below in the context of estimating a source wavelet in the case of a sparse impulse response green's function obtained by inversion.
When the inversion yields gestThen, according to equation 8, we estimate the source wavelet as the problem to be solved:
at this timeIs the unknown to be solved. We use a simple derivation to illustrate the approach of the present invention to source wavelet estimation. The source wavelet is estimated, i.e. the fourier spectrum of the source wavelet frequency domain is estimated. Assume that the two complex-valued vectors a and b satisfy the following relationship:
ac＝b (17)
and c is a Fourier transform spectrum complex coefficient of the seismic source wavelet to be solved. Further we have:
aHac＝aHb (18)
so as to obtain the compound with the characteristics of,
the invention estimates the source wavelet which changes with the shot gather space position according to the derivation. According to equations 16 and 19, when the estimated source wavelet does not vary with shot gather spatial location (i.e., the ordinary wavelet estimate), there are:
when the estimated source wavelet changes with the shot gather space position (i.e. the modified wavelet estimation), there are:
wherein, S: (j,j) Representing elements on the diagonal of the source wavelet matrix, j representing variation with shot gather spatial position or shot gather number, vec (-) representing the transformation of a matrix into a column vector, H representing the complex conjugate transpose, subscriptRepresenting column j, subscript, taking a frequency slice matrixRepresenting taking all columns of a frequency slice matrix. In the process of estimating the source wavelet, the invention does not need any assumption condition, namely, does not need any assumption on the phase of the wavelet. In addition, compared with the traditional wavelet estimation problem without considering the multiple waves, the multiple wave term in the wavelet estimation process helps to overcome the non-uniqueness problem (amplitude and phase) existing in the traditional wavelet estimation.
FIG. 2 is a flow chart of a preferred embodiment of the estimation method of the present invention, which comprises the steps of:
a) the following data were obtained: data d, threshold μ, gestMaximum number of iterative updates kmaxNumber of internal iterations l of the linear Bregman algorithmmax，sestMaximum number of iterative update times mmax。
b) Initialization, iterative variables k is 0, m is 0, b is d, gestAnd sestThe initialization is 0.
c) Judging k is more than kmaxAnd if so, repeatedly executing the step c-i, otherwise, executing the step j.
d) According to x0＝gest、sestAnd b ═ d, calling Algorithm 1 to update to get gest。
e) Update k to k + lmax。
f) M is judged to be less than mmaxAnd if yes, executing the step g-i.
g) According to b and updated gestCalling the formula 20 or 21 to update to obtain sest. Calculating each diagonal matrix by equation 20 or 21According to all the calculated values of all the elements on the diagonal lineUpdating sest。
h) Judging m as mmax-1 if true, k if truemax-k is assigned to lmax。
i) M +1 is assigned to m.
j) And obtaining the estimation results of the source wavelet and the surface multiple according to the final updated result.
There are two unknowns in the process: slices s of source waveletsestAnd section g of green functionestSimilar to the traditional blind deconvolution method in solving the inverse problem, i.e. fixing sestInversion of gest(ii) a Followed by fixing gestTo estimate sest. In the inversion of gestIn time, the inversion result of the previous time is used for each iteration updating. In inverting sparse gestDuring solution, iteration is often needed for several times, and according to the sparse recovery theory, iteration is often needed for several times to converge to obtain an optimal sparse solution. In order to obtain a sparse solution gestAnd then, estimating the source wavelet. Estimated sestImmediately for inverse updating gest. The whole iteration main loop exits and ends after a certain stopping criterion is met. Alternately updating gestAnd sestIs determined by kmaxAnd mmaxWhen g isestK after reaching the maximum number of updatesmaxAnd sestWhen the maximum updating times are not reached, g is executedestAfter the last update, s is executed onceestIs updated to finish gestAnd sestIn which each time l of the linear Bregman algorithm is performedmaxEqual; when s isestK after reaching the maximum number of updatesmaxAnd gestWhen the maximum number of updates is not reached, s is executedestOnly the last linear Bregman algorithm is executed after the last update, and the last linear Bregman algorithm is executedmaxIs assigned a value of gestIs different from the number of updates already, while each previous execution of l of the linear Bregman algorithmmaxAre equal. In which the linear Bregman algorithm is performed each timemaxEqual last update gestI.e. the estimated result of g, last update sestI.e., the estimation of s, g is used belowestDenotes the estimation of g, sestThe estimation result of s is indicated.
In the presence of a catalyst to obtain gestAnd sestThen, the primary wave p ═ g is obtained by calculationest*sest。
From gestAnd inputting dataObtaining surface-related multiples m-rsurfgest*d。
According to the principle of a secondary seismic source, a green function g is obtainedestAnd obtaining the time domain slices of the surface-related multiples of different orders, such as the first-order surface multiples m, by calculating through a feedback iterative model after obtaining the primary wave p1The calculation formula of (2) is as follows:
m1＝rsurfgest*p＝rsurfgest*gest*sest(22)
time domain slice m of second order surface multiples2The calculation formula of (2) is as follows:
m2＝rsurfgest*m1＝(rsurf)2gest*gest*gest*sest(23)
the formula for calculating surface-related multiples of other orders and so on.
Because most of the calculation amount during the multiple estimation is consumed in matrix-matrix multiplication of each frequency slice, namely 8-type calculation, in order to further reduce the calculation time consumption of the multiple estimation, the method only utilizes the dominant frequency band with high reliability and signal-to-noise ratio of the seismic data during the multiple and seismic source wavelet estimation according to the band-limited characteristic of the seismic data (namely, the frequency components of the seismic data are mainly concentrated in a certain frequency band with limited length), thereby greatly reducing the calculation time consumption of the method, which is different from the traditional method which utilizes all frequency components (including low reliability and signal-to-noise ratio).
After the first wave or the surface multiple waves are obtained, the underground structure map can be obtained by using the first wave or the surface multiple waves respectively.
Salt dome model data example
The data generated by simulation of the DELPHI multiple research team of Dutch Delff's theory of technology is shown in FIG. 4a, which is input by a velocity model and a density model, and is simulated by finite difference forward modeling, wherein the time sampling interval is 4 ms. The source wavelet used in the forward modeling is a Ricker wavelet, and each source wavelet is the same, namely a constant wavelet. The shallow layer of the model is a layer of water about 200 meters deep. For comparison with the REPSI method from Lin and Herrmann, a total number of iterations is set to 71 during the process. The results of the REPSI method are provided by Lin and Herrmann. The frequency band used in the present invention is in the range of 0-70Hz (the original whole frequency band is in the range of 0-125 Hz).
The present invention achieves better results from the processing results of the present invention in fig. 3 to 4. In the raw input data of fig. 4a, multiples and the same secondary are mixed together, interfere with each other, and are not clearly resolved. After the treatment by the method, the multiples and the primaries can be accurately separated. The weak event belonging to the fault in the raw data of fig. 4a is difficult to identify and is clearly identified and visualized in fig. 4c after the processing of the present invention. Since the present invention is directed to processing surface-related multiples, then interbed multiples will be included in the primary resulting after processing, see fig. 4 c.
By comparing the results of the present invention (fig. 4c) with the results of the REPSI methods by Lin and Herrmann (fig. 4d), it can be easily found that the primary wave obtained by the processing of the present invention is more accurate, and the multiples can be better estimated and suppressed; some multiples in-phase axes are still visible in the results of the REPSI method from Lin and Herrmann (FIG. 4d), as indicated by the arrows in FIG. 4 d.
Comparing the results of the estimated source wavelet of the present invention in fig. 5a and the results of the REPSI method estimation of Lin and Herrmann (the solid black line in fig. 5(a) represents the estimated source wavelet of the present invention, and the dashed line represents the source wavelet estimated by the REPSI method of Lin and Herrmann), it can be seen that the results of the two are consistent, which verifies the correctness of the method of the present invention from the side. Since the data simulation is generated by using the constant wave which does not change with the shot gather, the result generated when the invention carries out the estimation of the variable wave is more consistent, and is shown in figure 5 b.
Example of Pluto1.5 data
To further demonstrate the effectiveness of the present invention, we applied it to another convincing pluto1.5 data. This data, generated and released by simulation by the international SMAART JV organization, is one of the international signature data for open test multiple algorithms. The source wavelet used in the forward simulation of the data is Ricker wavelet, and each source wavelet is the same, namely, a constant wavelet. The time sampling interval is 8 milliseconds. We directly utilized the raw data released publicly by the SMAART JV organization without any pre-processing of the raw data except for the cutting of the direct wave. From the raw input data in fig. 6a and 7a, there is a large amount of multiple energy in the data. The frequency band used in the process of the invention is in the range of 0-50 Hz.
From the processing results in fig. 6 and 7, in fig. 6a and 7a of the raw input data, multiples and the same secondary are mixed together, interfere with each other, and are not clearly resolved, whereas the present invention achieves accurate primary and multiple estimation and separation (especially indicated by arrows in fig. 7). Since the source wavelet used in the data simulation generation is not changed with the shot gather, the result generated when the invention carries out the estimation of the changed wavelet is more consistent, and is shown in figure 8.
Figure 9 shows another scheme for evaluating the effectiveness of the multiple estimation technique. The essence of fig. 9 is that the feedback iteration model, namely: the first order surface-related multiples (fig. 9b) estimated by the present invention are used as the secondary source to transmit downward to predict the first order surface-related multiples (fig. 9c), the first order surface-related multiples (fig. 9c) are used as the secondary source to transmit downward to predict the second order surface-related multiples (fig. 9d), the second order surface-related multiples (fig. 9d) are used as the secondary source to transmit downward to predict the third order surface-related multiples (fig. 9e), the third order surface-related multiples (fig. 9e) are used as the secondary source to transmit downward to predict the fourth order surface-related multiples (fig. 9f), the fourth order surface-related multiples (fig. 9f) are used as the secondary source to transmit downward to predict the fifth order surface-related multiples (fig. 9g), and so on. For this data, the predicted surface-related multiples of different orders (fig. 9c to 9g) can be easily identified in the original input data (fig. 9 a). The summation of the estimated primaries and surface related multiples of different orders in the present invention yields a total data (fig. 9h) that is nearly identical to the original input data (fig. 9a), with little difference (fig. 9 i).
FIG. 10 further illustrates the results of the invention applied to other shot gathers in the Pluto1.5 data, and it can be seen that the primary and multiples are more accurately estimated and separated.
While the present invention has been described with reference to the embodiments shown in the drawings, the present invention is not limited to the embodiments, which are illustrative and not restrictive, and various modifications can be made by those skilled in the art without departing from the spirit and scope of the present invention as defined in the appended claims.
Claims (8)
1. A surface multiple and wavelet estimation method based on a linear Bregman algorithm is characterized by comprising the following steps:
s1, alternately updating the green function to be estimated and the source wavelet to be estimated using the following steps S11 and S12 until reaching a preset number of times, wherein the source wavelet to be estimated, which is first performed in step S11, has a preset initial value:
s11, updating the Green function to be estimated by adopting a linear Bregman algorithm according to the seismic source wavelet to be estimated and the acquired input data;
s12, updating the seismic source wavelet to be estimated by adopting the following formula according to the updated Green function to be estimated and the input data;
when the source wavelet to be estimated does not change along with the shot gather space position:
when the source wavelet to be estimated changes along with the shot gather space position:
s2, obtaining at least one estimation result of the seismic source wavelet and the surface multiple according to the final updated result;
wherein the superscript est represents the value to be estimated,representing elements on the diagonal of a slice of the source wavelet matrix, j representing the shot gather spatial position or shot gather number, vec () representing the transformation of a matrix into a column vector, subscriptRepresenting column j, subscript, taking a frequency slice matrixRepresenting taking all columns of a frequency slice matrix,for a slice of the matrix expression of the frequency domain of the input data,being a slice of a matrix expression of the frequency domain of the seismic source sub-waves, rsurfIs the reflection coefficient of the contact surface of air and water,a slice of a matrix expression of the frequency domain of the green's function.
2. The method of estimating surface multiples and wavelets according to claim 1, wherein updating the green's function to be estimated using a linear Bregman algorithm based on the source wavelet to be estimated and the acquired input data specifically comprises:
according to the formula yl+1＝yl-tlAHΠσ(Axl-b) and xl+1＝shrink(yl+1μ) are iterated to find y in turn1、x1、y2、x2、y3…、Wherein,a vector expression form l in the time domain for the update result of the green function to be estimated after the current execution of step S11maxThe maximum iteration number of the linear Bregman algorithm when the step S11 is executed this time;
in the formula, x0And y0A vector representation form of the time domain which is equal to the preset initial value of the green function to be estimated, an observation data vector b is the vector representation form of the time domain of the input data, tlWhat is represented is the step size, which is expressed as:
<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) is defined as:
<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>
the soft threshold function shrnk (x, μ) is defined as:
shrink(x,μ)＝max(|x|-μ,0)sign(x)，
wherein, σ represents a noise level factor, the max () function represents taking a maximum value, the | | function represents taking an absolute value, and sign (x) represents a symbol value taking function;
the matrix operator a represents a set of operations that perform the following operations: fourier transformation is carried out on the Green's function to be estimated from a vector table form of a time domain to a vector table form of a frequency domain, the Green's function to be estimated is changed from the vector table form of the frequency domain to a matrix table form of the frequency domain, and the Green's function to be estimated is transformed from the vector table form of the frequency domain to a matrix table form of the frequency domain according to a formulaCalculating each slice of the matrix expression form of the frequency domain of the analog data corresponding to the input dataAccording to each sliceAnd performing inverse Fourier transform to obtain a vector expression form of a time domain of the analog data corresponding to the input data.
3. The method for estimating surface multiples and wavelets according to claim 1, wherein in step S2, the estimation result of the source wavelet is obtained according to the result of the last update of the source wavelet to be estimated;
the estimation result of the surface multiple is based on m ═ rsurfg x d, m is a vector expression form of a time domain of the surface multiple, g is a vector expression form of a time domain of the green function, and d is a vector expression form of a time domain of the input data.
4. The method of estimating surface multiples and wavelets according to claim 2, wherein if the green' S function to be estimated has reached the maximum update times and the source wavelet to be estimated has not reached the maximum update times, then when step S11 is performed again, lmaxUpdated to the difference of the maximum update times of the source wavelet to be estimated minus the times of update already, and executed each time beforeL at step S11maxEqual; if the source wavelet to be estimated has reached the maximum update times and the Green function to be estimated has not reached the maximum update times, then step S12 is performed only once, and l is performed each time step S11 is performedmaxAre equal.
5. A surface multiple and wavelet estimation system based on a linear Bregman algorithm, comprising:
the alternative updating unit is used for alternatively updating the Green function to be estimated and the seismic source wavelet to be estimated by using the following Green function updating unit and the seismic source wavelet updating unit until the preset times are reached, wherein the seismic source wavelet to be estimated, which calls the Green function updating unit for the first time, has a preset initial value:
the green function updating unit is used for updating the green function to be estimated by adopting a linear Bregman algorithm according to the seismic source wavelet to be estimated and the acquired input data;
the seismic source wavelet updating unit is used for updating the seismic source wavelet to be estimated by adopting the following formula according to the updated green function to be estimated and the input data;
when the source wavelet to be estimated does not change along with the shot gather space position:
when the source wavelet to be estimated changes along with the shot gather space position:
the estimation result acquisition unit is used for obtaining an estimation result of at least one of the seismic source wavelet to be estimated and the surface multiple to be estimated according to the final updated result;
wherein the superscript est represents the value to be estimated,representing elements on the diagonal of a slice of the source wavelet matrix, j representing the shot gather spatial position or shot gather number, vec () representing the transformation of a matrix into a column vector, subscriptRepresenting column j, subscript, taking a frequency slice matrixRepresenting taking all columns of a frequency slice matrix,for a slice of the matrix expression of the frequency domain of the input data,being a slice of a matrix expression of the frequency domain of the seismic source sub-waves, rsurfIs the reflection coefficient of the contact surface of air and water,a slice of a matrix expression of the frequency domain of the green's function.
6. The surface multiples and wavelet estimation system of claim 5, wherein updating the green's function to be estimated using the linear Bregman algorithm based on the source wavelet to be estimated and the acquired input data specifically comprises:
according to the formula yl+1＝yl-tlAHΠσ(Axl-b) and xl+1＝shrink(yl+1μ) are iterated to find y in turn1、x1、y2、x2、y3…、Wherein,a vector expression form l in the time domain for the update result of the green function to be estimated after the current execution of step S11maxThe maximum iteration number of the linear Bregman algorithm when the step S11 is executed this time;
in the formula, x0And y0A vector representation form of the time domain which is equal to the preset initial value of the green function to be estimated, an observation data vector b is the vector representation form of the time domain of the input data, tlWhat is represented is the step size, which is expressed as:
<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) is defined as:
<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>
the soft threshold function shrnk (x, μ) is defined as:
shrink(x,μ)＝max(|x|-μ,0)sign(x)，
wherein, σ represents a noise level factor, the max () function represents taking a maximum value, the | | function represents taking an absolute value, and sign (x) represents a symbol value taking function;
the matrix operator a represents a set of operations that perform the following operations: fourier transformation is carried out on the Green's function to be estimated from a vector table form of a time domain to a vector table form of a frequency domain, the Green's function to be estimated is changed from the vector table form of the frequency domain to a matrix table form of the frequency domain, and the Green's function to be estimated is transformed from the vector table form of the frequency domain to a matrix table form of the frequency domain according to a formulaCalculating each slice of the matrix expression form of the frequency domain of the analog data corresponding to the input dataAccording to each sliceObtaining a time-domain vector expression form of analog data corresponding to the input data by performing inverse Fourier transform。
7. The surface multiples and wavelets estimation system of claim 5 wherein, in said estimation result acquisition unit,
the estimation result of the seismic source wavelet is obtained according to the result of updating the seismic source wavelet to be estimated for the last time;
the estimation result of the surface multiple is based on m ═ rsurfg x d, m is a vector expression form of a time domain of the surface multiple, g is a vector expression form of a time domain of the green function, and d is a vector expression form of a time domain of the input data.
8. The surface multiples and wavelet estimation system of claim 6, wherein if the Green' S function to be estimated has reached the maximum number of updates and the source wavelet to be estimated has not reached the maximum number of updates, then when step S11 is performed again,/maxUpdated to the maximum number of updates of the source wavelet to be estimated minus the number of updates already, and l is previously set each time step S11 is executedmaxEqual; if the source wavelet to be estimated has reached the maximum update times and the Green function to be estimated has not reached the maximum update times, then step S12 is performed only once, and l is performed each time step S11 is performedmaxAre equal.
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 |
|---|---|---|
| CN107728211B (en) | Seismic signal recovery algorithm based on tensor nuclear norm regularization | |
| US11333783B2 (en) | Integrated method for estimation of seismic wavelets and synthesis of seismic records in depth domain | |
| Yuan et al. | Double-scale supervised inversion with a data-driven forward model for low-frequency impedance recovery | |
| US10429530B2 (en) | Deghosting with adaptive operators | |
| CN112946749B (en) | Method for suppressing seismic multiples based on data augmentation training deep neural network | |
| CN107894612B (en) | A kind of the sound impedance inversion method and system of Q attenuation by absorption compensation | |
| WO2014189679A1 (en) | Multi-parameter inversion through offset dependent elastic fwi | |
| CA2825395A1 (en) | Convergence rate of full wavefield inversion using spectral shaping | |
| CN102884447A (en) | Q tomography method | |
| CN105467442B (en) | The time-varying sparse deconvolution method and device of global optimization | |
| CN108535775A (en) | Non-stationary seismic data sound impedance inversion method and device | |
| CN107367760B (en) | Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm | |
| CN109143331B (en) | Seismic wavelet extraction method | |
| EP4031910B1 (en) | Noise attenuation methods applied during simultaneous source deblending and separation | |
| Donno et al. | Curvelet-based multiple prediction | |
| CN108508489B (en) | Seismic inversion method based on waveform micro-variation matching | |
| Ibrahim et al. | Wavefield reconstruction using a Stolt-based asymptote and apex shifted hyperbolic Radon transform | |
| CN107390261B (en) | Surface-related multiple and higher-order spectra method and system based on linear Bregman algorithm | |
| CN114966860A (en) | Seismic data denoising method based on convolutional neural network | |
| Aghamiry et al. | Complex-valued imaging with total variation regularization: an application to full-waveform inversion in visco-acoustic media | |
| CN117724154B (en) | VTI medium parameter prediction method considering strong anisotropy characteristics | |
| Cheng et al. | Meta-Processing: A robust framework for multi-tasks seismic processing | |
| CN110146923A (en) | A kind of efficient high accuracy depth domain methods of seismic wavelet extraction | |
| Chen et al. | Joint data and model-driven simultaneous inversion of velocity and density | |
| CN108680957B (en) | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting |
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 |