CN107367760B - Based on the surface-related multiple and higher-order spectra method and system for accelerating linear Bregman algorithm - Google Patents

Abstract

The invention discloses a kind of based on the surface-related multiple and higher-order spectra method and system that accelerate linear Bregman algorithm, feedback model based on description primary wave and surface-related multiple relationship, this method carries out sparse inversion first with the linear Bregman algorithm of acceleration, to obtain the impulse response unrelated with surface, then the source wavelet with big gun collection spatial position change is updated using a step linear spectral estimation formulas, and the impulse response and source wavelet unrelated with surface are updated by cycle alternation, final estimation obtains primary wave and surface-related multiple.Compared to traditional multi-parameter, the surface-related multiple estimation method based on Chang Zibo hypothesis and complicated algorithm, implement estimation method and system of the invention, the source wavelet with big gun collection spatial position change can be estimated, and high-efficient, the time-consuming precision for less, estimating result of estimation method and system is high.

Description

Surface multiple and wavelet estimation method based on accelerated linear Bregman algorithm and system

technical field

The invention relates to the field of seismic exploration, in particular to seismic signal processing in seismic exploration technology, and more particularly to a method and system for estimating surface multiples and wavelets based on accelerated linear Bregman algorithm.

Background technique

The problem related to multiple reflections (multiples) is one of the most prominent problems commonly existing in reflection seismic exploration. Because the sea surface is a strong reflection interface, the problem of surface-related multiples is particularly serious in marine seismic exploration. Today, the exploration and development of offshore oil and natural gas combustible ice resources has become the top priority of major oil and energy companies at home and abroad. The South my country Sea is rich in oil and natural gas combustible ice resources, which is the hot spot and focus of current ocean exploration. In the current oil industry production, a basic assumption of conventional seismic data processing technology is that the input reflection data volume is composed of only one reflected wave (primary wave). Multiples can be mistakenly identified as primary waves or as part of a mixture of primary waves. The multiple waves thus cause strong interference to the processing and imaging of the primary waves and affect the resource exploration. Therefore, the multiple wave estimation technology has excellent practical industrial application prospects.

At present, the multiple wave estimation or suppression methods used in the industry are mainly divided into the following categories: one type uses the characteristic difference between multiple waves and primary waves to suppress multiple waves, which is called filtering method; one type is first predicted from seismic data. Multiples are then adaptively subtracted from the original data, known as predictive subtraction; there is also an inverse scattering series method to suppress multiples.

The filtering method assumes that multiples and primarys have obvious characteristic differences in a specific transform domain, and designs various mathematical algorithms to obtain such differences and eliminate multiples. The main characteristic differences exploited by filtering methods include the periodicity and separability of the multiples. When the filtering method is applied to the attenuation of some complex structures of multiples, the characteristic difference between the multiple and the primary wave is small or no difference. To satisfy, the multiple wave suppression effect is not obvious.

There are also some problems in the prediction of multiples by inverse scattering series: large amount of calculation and high cost; it is difficult to predict multiples at long offsets; the prediction ability of multiples in complex media also has certain limitations. Among them, the high computational cost makes this method difficult to apply to actual data processing.

Berkhout proposed a feedback model theory to describe complex multiples in 1982, which laid the mathematical and physical foundation of the feedback iterative multiple suppression method, but to eliminate multiples requires known source wavelets. In order to illustrate this problem and facilitate subsequent descriptions, we first introduce the feedback model formula that is closely related to the present invention, that is,

Among them, D is the matrix expression form of the collected input data in the time domain, is a slice of the matrix expression in the frequency domain after the forward Fourier transformation of D, S is the matrix expression in the time domain of the source wavelet to be estimated, is a slice of the matrix expression in the frequency domain after the forward Fourier transformation of S, G is the matrix expression in the time domain of the Green's function (that is, the subsurface impulse response without multiples), It is a slice of the matrix expression in the frequency domain after the forward Fourier transformation of G. We use the ^ symbol to represent the data or variables in the frequency domain, r _surf is the reflection coefficient of the contact surface between air and water, assuming r _surf =-1 ,in and The frequency-domain product of (also called the time-domain convolution) represents a slice of the estimated frequency of the primary wave, namely in is a slice of the frequency domain of the primary wave. The primary waves involved in the present invention refer to all waves except the surface multiples, that is, the primary waves here include the interlayer multiples (which have nothing to do with the surface). represents a slice of the frequency domain of the estimated surface multiples, i.e. the data Propagating downward as a secondary source, the same subsurface impulse responds to the Green's function The frequency-domain product of (time-domain convolution) is used to predict slices in the frequency domain of surface-correlated multiples. It can be seen from Equation 1 that accurate surface multiple estimation depends on the accuracy of the estimated source wavelets.

Verschuur et al. proposed the surface correlated multiple suppression technique SRME in 1992, which estimated the source wavelet from the actual data by introducing an adaptive subtraction based on the assumption of minimum primary wave energy. In order to avoid the nonlinear problem caused by solving the source wavelet, Berkhout and Verschuur proposed the iterative SRME method in 1997. The premise of the minimum primary wave energy is not satisfied in all cases. In addition, in practical applications, SRME also has the problem that the effective signal is damaged by adaptive prediction and subtraction.

In order to improve the problem that SRME has the problem of predicting subtractive damage to the effective signal, van Groenestijn and Verschuur proposed the sparse inversion primary wave estimation technique, or EPSI, in 2009. In order to improve the existing problems in EPSI, Lin and Herrmann proposed a robust sparse inversion primary wave estimation technique, namely REPSI, in 2013. EPSI and REPSI are the prior art relatively related to the present invention, and their technical solutions, existing problems and shortcomings will be described in detail below.

In addition, when the actual field seismic data is collected, the source wavelets excited by each shot are not the same, and all the current multiple wave estimation methods assume that all the source wavelets are the same, which affects the multiple wave estimation accuracy. The source wavelet can be estimated as a function of the spatial position of the shot collection.

Since each iteration consumes a lot of computing resources (memory, processor) and takes a long time to estimate multiples, in addition to improving the estimation accuracy of multiples and source wavelets, the efficiency of the method and the convergence of the algorithm are improved. The rate, the complexity of the simplified algorithm, and the adjustment parameters of the reduction method are also the development trend of multiple estimation technology.

The prior art related to the present invention is described below.

The technical solution of the prior art

In 2009, van Groenestijn and Verschuur proposed the sparse inversion primary wave estimation method EPSI based on the theoretical basis of SRME. In the iterative inversion process, the method is driven by the residual of the fitted data and alternately updates the subsurface pulse through the steepest descent method. Response to two unknown parameters G and source wavelet S, and then reconstruct the primary wave and its corresponding surface-correlated multiples Compared with SRME, EPSI avoids the step of adaptive subtraction and better protects the effective signal of primary wave.

The EPSI method of van Groenestijn and Verschuur assumes that all source wavelets are the same when estimating source wavelets, and does not provide an estimation scheme for source wavelets that vary with the spatial position of the shot collection (ie, constant wavelet estimation). The invention can estimate the source wavelet which changes with the spatial position of the shot collection.

In order to ensure the sparsity of the subsurface impulse response G in the time domain, EPSI updates the impulses in the order of time from small to large and amplitude from strong to weak in the process of updating the impulse response G. That is to say, in order to promote the time-domain sparsity of the underground impulse response G, the EPSI method has designed many artificial adjustment parameters. For example, EPSI needs to determine: a time window to determine the update range of the Green's function (that is, the update within the time window, zeroing outside the time window); the number of reflection events in the Green's function per iteration, etc. As Lin and Herrmann pointed out in their 2013 paper: Even if the number of subsurface reflecting interfaces is known precisely (which we cannot actually know), the EPSI method cannot fully utilize this prior information. The update of the subsurface impulse response of EPSI is easily affected by parameters such as the pulse selection window and the number of pickups, which in turn affects the estimation accuracy of multiples. Essentially, the sparsity of the Green's function in the EPSI method is controlled by humans rather than by mathematical optimization algorithms. This shortcoming limits the application of EPSI method in practical industry.

Prior art two related to the present invention

Based on the same theory of EPSI (that is, feedback model 1), Lin and Herrmann in 2013 modified the sparse constraint of EPSI subsurface impulse response to a one-norm constraint, and adopted the spectral mapping gradient one-norm algorithm To solve the sparse underground impulse response Green's function, 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, also alternately updates the two unknown parameters of the subsurface impulse response G and the source wavelet S, and then reconstructs the primary wave and its corresponding surface-correlated multiples The biggest difference between EPSI and REPSI is that the sparsity of Green's function in EPSI is mostly controlled by artificial adjustment parameters, while REPSI uses Mathematical optimization algorithm to solve the Green's function to obtain the sparse subsurface impulse response.

Lin and Herrmann's REPSI method also assumes that all source wavelets are the same when estimating source wavelets, and does not provide an estimation scheme for source wavelets (ie, constant wavelet estimation) that varies with the spatial position of the shot collection. The invention can estimate the source wavelet which changes with the spatial position of the shot collection. In addition, REPSI uses the LSQR algorithm to iteratively solve the source wavelet when estimating the constant wavelet, which greatly increases the computational complexity and reduces the efficiency of the method. In the present invention, the source wavelet can be estimated and obtained in one step.

Another important disadvantage of Lin and Herrmann's REPSI method is that it relies on an extremely complex gradient-norm algorithm for spectral mapping to solve the Green's function for the sparse subsurface impulse response. is an open source MATLAB program that can be downloaded by visiting https://github.com/mpf/spgl1. The complexity of the algorithm is mainly reflected in: many adjustment parameters, many sub-functions, and thousands of lines of code. Schlumberger oilfield services tried to adapt apply it in the actual industry, but due to the The complexity of the algorithm was ultimately unsuccessful. The complexity of the algorithm also reduces the efficiency of the REPSI method and hinders the practical application of the REPSI method in the industry.

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is to provide a method and system for estimating surface multiples and wavelets based on the accelerated linear Bregman algorithm, aiming at the above-mentioned technical defects of EPSI and REPSI in the prior art.

According to one aspect of the present invention, in order to solve its technical problem, the present invention provides a method for estimating surface multiples and wavelets based on accelerated linear Bregman algorithm, comprising the following steps:

S1, use the following steps S11 and S12 to alternately update the Green's function to be estimated and the source wavelet to be estimated until reaching a preset number of times, wherein the source wavelet to be estimated that performs step S11 for the first time has a preset initial value:

S11. According to the source wavelet to be estimated and the collected input data, use the accelerated linear Bregman algorithm to update the Green's function to be estimated;

S12, according to the updated Green's function to be estimated and the input data, adopt the following formula to update the source wavelet to be estimated;

When the source wavelet to be estimated does not change with the spatial position of the shot collection:

When the source wavelet to be estimated changes with the spatial position of the shot collection:

S2. According to the final updated result, obtain the estimation result of at least one of the source wavelet and the surface multiple waves;

Among them, the superscript est means to be estimated, Represents the elements on the diagonal of the source wavelet matrix slice, j represents the spatial position of the shot set or the shot set number, vec() represents the transformation of a matrix into a column vector, and the subscript Represents taking the jth column of a frequency slice matrix, subscript represents taking all columns of a frequency slice matrix, is a slice of a matrix expression in the frequency domain of the input data, is the slice of the matrix expression in the frequency domain of the source wavelet, r _surf is the reflection coefficient of the contact surface between air and water, A slice of the matrix expression in the frequency domain for the Green's function.

Further, in the method for estimating surface multiples and wavelets based on the accelerated linear Bregman algorithm of the present invention, according to the source wavelet to be estimated and the collected input data, the accelerated linear Bregman algorithm is used to update the Green's function to be estimated, which specifically includes: :

According to the formula and Iterate to find y ₁ , x ₁ , y ₂ , x ₂ , y ₃ …, in, is the vector representation in the time domain of the update result of the Green's function to be estimated after step S11 is executed this time, and ^lmax is the maximum number of iterations of the accelerated linear Bregman algorithm when step S11 is executed this time;

In the formula, x ₀ , and y ₀ are both equal to the vector representation in the time domain of the preset initial value of the Green's function to be estimated, and the observation data vector b is the vector representation in the time domain of the input data, t represents the step size, and its expression is:

The mapping function ∏ _σ (Ax _l -b) is defined as:

The soft threshold function shrink(x, μ) is defined as:

shrink(x,μ)=max(|x|-μ,0)sign(x),

Among them, σ represents the noise level factor, the max() function represents the maximum value, the || function represents the absolute value, and sign(x) represents the sign value function;

The matrix operator A represents a set of operations to perform the following operations: Fourier transform the Green's function to be estimated from a vector table in the time domain to a vector table in the frequency domain, and convert the Green's function to be estimated from the frequency domain. The vector table form becomes the matrix table form in the frequency domain, according to the formula Calculate each slice of the matrix representation of the frequency domain of the analog data corresponding to the input data separately According to each slice Inverse Fourier transform is performed to obtain the vector representation in the time domain of the analog data corresponding to the input data.

Further, in the method for estimating surface multiples and wavelets based on the accelerated linear Bregman algorithm of the present invention, 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 multiples is obtained according to m=r _surf g*d, where m is the vector representation in the time domain of the surface multiples, g is the vector representation in the time domain of the Green’s function, and d is the time domain of the input data vector representation of .

Further, in the method for estimating surface multiples and wavelets based on the accelerated linear Bregman algorithm of the present invention, 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 carry out again. In step S11, ^lmax is updated to the difference between the maximum update times of the source wavelet to be estimated minus the number of times that have been updated, and the ^lmax when step S11 is executed before each time is equal; if the source wavelet to be estimated If the maximum number of updates has been reached but the Green's function to be estimated has not reached the maximum number of updates, then step S12 is performed only once, and ^lmax is the same each time step S11 is performed.

In order to solve the technical problem, the present invention also provides a surface multiple wave and wavelet estimation system based on accelerated linear Bregman algorithm, including:

The alternate update unit is used to alternately update the Green's function to be estimated and the source wavelet to be estimated by using the following Green's function update unit and the source wavelet update unit until a preset number of times is reached, wherein the Green's function update unit to be estimated is called for the first time. The hypocenter wavelet has preset initial values:

The Green's function updating unit adopts the accelerated linear Bregman algorithm to update the Green's function to be estimated according to the source wavelet to be estimated and the collected input data;

The source wavelet updating unit adopts the following formula to update the source wavelet to be estimated according to the updated Green's function to be estimated and the input data;

When the source wavelet to be estimated does not change with the spatial position of the shot collection:

When the source wavelet to be estimated changes with the spatial position of the shot collection:

an estimation result obtaining unit, configured to obtain an estimation result of at least one of the source wavelet to be estimated and the surface multiple to be estimated according to the final updated result;

Among them, the superscript est means to be estimated, Represents the elements on the diagonal of the source wavelet matrix slice, j represents the spatial position of the shot set or the shot set number, vec() represents the transformation of a matrix into a column vector, and the subscript Represents taking the jth column of a frequency slice matrix, subscript represents taking all columns of a frequency slice matrix, is a slice of a matrix expression in the frequency domain of the input data, is the slice of the matrix expression in the frequency domain of the source wavelet, r _surf is the reflection coefficient of the contact surface between air and water, A slice of the matrix expression in the frequency domain for the Green's function.

Further, in the surface multiple wave and wavelet estimation system based on the accelerated linear Bregman algorithm of the present invention, according to the source wavelet to be estimated and the collected input data, the accelerated linear Bregman algorithm is used to update the Green's function to be estimated, which specifically includes: :

According to the formula and Iterate to find y ₁ , x ₁ , y ₂ , x ₂ , y ₃ …, in, is the vector representation in the time domain of the update result of the Green's function to be estimated after step S11 is executed this time, and ^lmax is the maximum number of iterations of the accelerated linear Bregman algorithm when step S11 is executed this time;

In the formula, x ₀ , and y ₀ are both equal to the vector representation in the time domain of the preset initial value of the Green's function to be estimated, and the observation data vector b is the vector representation in the time domain of the input data, t _l represents the step size, and its expression is:

The mapping function ∏ _σ (Ax _l -b) is defined as:

The soft threshold function shrink(x, μ) is defined as:

shrink(x,μ)=max(|x|-μ,0)sign(x),

Among them, σ represents the noise level factor, the max() function represents the maximum value, the || function represents the absolute value, and sign(x) represents the sign value function;

The matrix operator A represents a set of operations to perform the following operations: Fourier transform the Green's function to be estimated from a vector table in the time domain to a vector table in the frequency domain, and convert the Green's function to be estimated from the frequency domain. The vector table form becomes the matrix table form in the frequency domain, according to the formula Calculate each slice of the matrix representation of the frequency domain of the analog data corresponding to the input data separately According to each slice Inverse Fourier transform is performed to obtain the vector representation in the time domain of the analog data corresponding to the input data.

Further, in the surface multiple and wavelet estimation system based on the accelerated linear Bregman algorithm of the present invention, in the estimation result obtaining unit,

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 multiples is obtained according to m=r _surf g*d, where m is the vector representation in the time domain of the surface multiples, g is the vector representation in the time domain of the Green’s function, and d is the time domain of the input data vector representation of .

Further, in the surface multiple wave and wavelet estimation system based on the accelerated linear Bregman algorithm of the present invention, 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 carry out again. In step S11, ^lmax is updated to the difference between the maximum update times of the source wavelet to be estimated minus the number of times that have been updated, and the ^lmax when step S11 is executed before each time is equal; if the source wavelet to be estimated If the maximum number of updates has been reached but the Green's function to be estimated has not reached the maximum number of updates, then step S12 is performed only once, and ^lmax is the same each time step S11 is performed.

The estimation method and system based on the accelerated linear Bregman algorithm of the present invention have the following beneficial effects: the estimation method and system of the present invention can estimate the source wavelet that varies with the spatial position of the shot collection, and the estimation method and system have high efficiency , less time-consuming, and high accuracy of the estimated results.

Description of drawings

The present invention will be further described below in conjunction with the accompanying drawings and embodiments, in which:

1 is a flowchart of a preferred embodiment of the accelerated linear Bregman algorithm of the present invention to solve Ax=b to obtain a sparse solution x;

2 is a flowchart of a preferred embodiment of the estimation method of the present invention for estimating source wavelets and surface multiples;

3 is a schematic diagram of the true solution, the solution obtained by some algorithms, and the error between the obtained solution and the true solution;

Figure 4 is a comparison diagram of the variation curve of the data residual obtained by some algorithms with the number of iterations;

Figure 5(a) is a schematic diagram of the input data including multiples;

Figure 5(b) is a schematic diagram of primary wave data estimated by the present invention;

Figure 5(c) is a schematic diagram of surface multiple data estimated by the present invention;

Figure 6(a) is a schematic diagram of the zero-offset cross-section of the data in Figure 5(a);

Figure 6(b) is a schematic diagram of surface multiple data estimated by the present invention;

FIG. 6(c) is a schematic diagram of primary wave data estimated by an embodiment of the present invention;

Figure 6(d) is a schematic diagram of the primary wave data estimated by Lin and Herrmann's REPSI method;

Figure 7(a) is a schematic diagram of the estimated source wavelet;

Fig. 7(b) is a schematic diagram of the frequency domain amplitude spectrum of the source wavelet estimated with the shot collection according to an embodiment of the present invention;

Fig. 8 is a graph showing the inversion error of Lin and Herrmann's REPSI method, the original linear Bregman algorithm, and the accelerated linear Bregman algorithm with the number of iterations

Figure 9(a) is a schematic diagram of the Pluto 1.5 data released by the international SMAART JV organization for testing the multiple-wave suppression method;

Figure 9(b) is a schematic diagram of primary wave data estimated by the present invention;

Figure 9(c) is a schematic diagram of surface-related multiple data estimated by the present invention;

Fig. 10(a) is a schematic diagram of the zero-offset cross-section of the data in Fig. 9(a);

Figure 10(b) is a schematic diagram of primary wave data estimated by the present invention;

Figure 10(c) is a schematic diagram of surface-related multiple data estimated by the present invention;

Figure 11(a) is a schematic diagram of the source wavelet estimated from the Pluto 1.5 data as a function of shot collection;

Figure 11(b) is a schematic diagram of the frequency domain amplitude spectrum of the source wavelet estimated from the Pluto 1.5 data as a function of shot collection;

Fig. 12 (a) is the data section schematic diagram of the middle gun of the data in Fig. 9 (a);

Figure 12(b) is a schematic diagram of primary wave data estimated by the present invention;

Fig. 12(c) is a schematic diagram of first-order surface multiple data estimated by the present invention;

Figure 12(d) is a schematic diagram of second-order surface multiple data estimated by the present invention;

Figure 12(e) is a schematic diagram of third-order surface multiple data estimated by the present invention;

Figure 12(f) is a schematic diagram of the fourth-order surface multiple data estimated by the present invention;

Figure 12(g) is a schematic diagram of fifth-order surface multiple data estimated by the present invention;

Figure 12(h) is a schematic diagram of the accumulation and summation of the data in Figures 12(b) to 12(g);

Figure 12 (i) is a schematic diagram of the difference between Figure 12 (a) and Figure 12 (h);

Figure 13(a) is a schematic diagram of the cross-section of the zero-scanning bureau with shot numbers between 701-880 in the Pluto 1.5 data;

Figure 13(b) is a schematic diagram of primary wave data estimated by the present invention;

Fig. 13(c) is a schematic diagram of surface-related multiple data estimated by the present invention.

noun description

Primary wave (Primaries): It means that the excited seismic wave is only reflected once at the subsurface interface during the propagation of the stratum, that is, it is received by the geophone buried on the surface, which is called the primary reflected wave, or the primary wave for short. The primary waves involved in the present invention refer to all waves except the surface multiples, that is, the primary waves here include the interlayer multiples (which have nothing to do with the surface).

Multiples: It means that the excited seismic waves are received by the geophone buried on the surface after multiple reflections on the underground interface or free surface (contact surface between air and water) during the propagation process. Secondary reflected waves, referred to as multiple waves.

Surface-related multiples (Surface-related multiples) or surface multiples (Surfacemultiples): refers to the multiple reflection waves related to the free surface (air and water contact surface), which are briefly described as surface multiples in the present invention Wave.

Interlayer multiples (Internal multiples, IM): refers to the multiple reflection waves generated between the subsurface reflection interfaces (between formations).

Source wavelet: refers to the wave generated by the source excitation.

Feedback model: Feedback model.

Surface-Related Multiple Elimination (SRME): Surface-Related Multiple Elimination.

Estimation of Primaries by Sparse Inversion (EPSI): Sparse primary wave estimation for inversion.

Robust Estimation of Primaries by Sparse Inversion (EPSI): Robust sparse inversion primary wave estimation.

: Spectral-Projected Gradient-norm algorithm for spectral mapping.

LSQR: Least-squares QR, least squares QR decomposition algorithm.

Linearized Bregman (LB): Linearized Bregman algorithm.

Accelerated Linearized Bregman (ALB): Accelerated Linearized Bregman algorithm.

MATLAB: It is a commercial mathematical software produced by MathWorks in the United States. It is a high-level technical computing language and interactive environment for algorithm development, data visualization, data analysis, and numerical computing. MATLAB is a combination of the words matrix&laboratory, which means matrix factory (matrix laboratory).

Signal-to-Noise Ratios (SNR): Signal-to-Noise Ratio.

Detailed ways

In order to have a clearer understanding of the technical features, objects and effects of the present invention, the specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings. In order to facilitate the understanding of the present invention, the above and the following expression letters are now described as follows, the lowercase letters s, g, d, p, m are the vector expressions representing the time domain, S, G, D, P, M is the matrix representation in the time domain. is a slice of the matrix representation of the frequency domain of S, G, D, P, and M, and the superscript est indicates to be estimated, It is the vector expression form of the frequency domain of s, g, d, p, and m. S, G, D, P, and M are three-dimensional matrices, and their slices are two-dimensional matrices. Each three-dimensional matrix can be divided into multiple 2D matrix. The vector representation of any two-dimensional matrix Z=[Z ₁ Z ₂ Z ₃ ... Z _n ] can be expressed as:

Among them, Z ₁ to Z _n represent the columns of the two-dimensional matrix Z. On the contrary, the matrix expression form can be obtained from the vector expression form.

The main formula used in the estimation of primary waves, surface multiples and source wavelets is the mathematical-physical relational formula of the feedback model, that is, formula (1). The excited source wavelet s is introduced into the subsurface and undergoes time domain convolution with the subsurface impulse response Green's function g (asterisk * indicates convolution) to generate a primary wave p. The mathematical relationship is expressed as:

p=g*s (2)

The generated primary wave p is reflected at the contact surface between seawater and air, propagates downward, and then re-introduces into the ground as a seismic source, where it undergoes time-domain convolution with the underground impulse response g to generate first-order multiple waves m ₁

m ₁ =g*r _surf p=r _surf g*g*s (3)

The generated first-order multiples m ₁ are reflected from the free surface, and then re-introduced into the ground as the source, where they are convolved with the subsurface impulse response g in the time domain to generate second-order multiples m ₂

m ₂ =g*r _surf m ₁ =(r _surf ) ² g*g*g*s (4)

The higher-order surface-related multiple generation formulas are analogous.

The total observed data d is the sum of the primary wave p and the surface-related multiples of all orders, and its mathematical expression is:

and then have

The mathematical-physical relationship expression of the above feedback model can be summarized as:

Expression (1) is the frequency domain matrix expression form of Expression 7 in the case of two-dimensional or three-dimensional, and the parameters in formula 7 are modified to the corresponding matrix expression form in the frequency domain. It is transcribed here:

The main purpose of the present invention is: according to expressions 7-8, solve the unknown variables: the source wavelet s and the Green's function g, and then obtain the estimated primary wave p=g*s and surface multiples m=r _surf g*d .

The present invention first describes a method for inverting a Green's function of a sparse impulse response given an estimated source wavelet.

Given the estimated source wavelet s ^est , according to expressions (7)-(8), we write the formula used for inversion g ^est , namely:

here is a simplified and representative forward calculus operator, the core of which is in the calculation formula 8 Specifically, the matrix operator Represents a set of operations that perform the following operations: Fourier transform the Green's function to be estimated from the vector table form g in the time domain into a vector table form in the frequency domain Convert the Green's function to be estimated in the form of a vector table in the frequency domain into a matrix table form in the frequency domain According to the formula Calculate each slice of the matrix representation of the frequency domain of the analog data corresponding to the input data separately Each slice of the matrix representation of the frequency domain of the analog data corresponding to the calculated input data Perform the inverse Fourier transform to obtain the vector representation in the time domain of the analog data corresponding to the input data

In order to further explain how the present invention solves and obtains the sparse Green's function g, formula 9 is rewritten into the following standard form:

Ax=b (10)

in, x=g, b=d.

In order to seek the sparse solution of the Green's function g, the present invention uses the accelerated linear Bregman algorithm to solve Equation 10. The accelerated linear Bregman algorithm is designed to solve:

Indicates that x is solved so that the expression achieves the minimum value, and subject to is under limited conditions.

Among them, μ≥0 is the sparsity control factor, ||x|| ₁ and ||x|| ₂ represent the vector norm sum norm, σ represents the noise level factor in the data. Algorithm 1 below gives the detailed iterative steps to speed up the linear Bregman algorithm.

Referring to FIG. 1, FIG. 1 is a flowchart of the accelerated linear Bregman algorithm of the present invention to solve Ax=b to obtain a sparse solution x, that is, a flowchart of Algorithm 1, and Algorithm 1 specifically includes the following steps:

a) Obtain the following parameters: forward operator A, observation data vector b, initial value vector x ₀ , threshold μ, maximum number of iterations l ^max .

b) Initialization: iterative variable l=0, assign x ₀ to y ₀

c) Repeat steps d to g below for 1 ^max times

d) Calculate

e) calculation

f) calculation

g) update l to l+1

h) get a sparse solution The sparse solution x is the result after the update of the accelerated linear ^Bregman algorithm this time.

In the above algorithm 1, t _l of step d represents the step size, and its expression is:

where the superscript H stands for complex conjugate transpose. The mapping function ∏ _σ (Ax _l -b) is to consider the noise existing in the data, which is defined as:

The soft threshold function shrink(x, μ) is defined as:

shrink(x,μ)=max(|x|-μ,0)sign(x) (14)

Among them, the max(·) function represents the maximum value, the |·| function represents the absolute value, and the sign(x) represents the sign value function.

In the above Algorithm 1, the definition expression of α _l is:

This method accelerates the convergence rate of the algorithm by weighting the system α _l .

The following description of the present invention will describe the estimation of the source wavelet under the condition that the Green's function of the sparse impulse response is obtained by inversion.

When g ^est is obtained by inversion, according to Equation 8, the problem to be solved when estimating the source wavelet is:

at this time is the unknown to be solved. We illustrate the scheme of the present invention for source wavelet estimation with a simple derivation. To estimate the source wavelet, that is, to estimate the Fourier spectrum of the source wavelet in the frequency domain. Suppose two complex-valued vectors a and b satisfy the following relationship:

ac=b (17)

where c is the complex spectral coefficient of the Fourier transform of the source wavelet to be obtained. Then we have:

a ^H ac = a ^H b (18)

to get,

According to the above derivation, the present invention estimates the source wavelet which changes with the spatial position of the shot collection. According to Equation 16 and Equation 19, when the estimated source wavelet does not change with the spatial position of the shot collection (that is, constant wavelet estimation), there are:

When the estimated source wavelet changes with the spatial position of the shot collection (that is, the variable wavelet estimation), there are:

Among them, S _{(j, j)} represents the element on the diagonal of the source wavelet matrix, j represents the change with the spatial position of the shot set or the shot set number, vec( ) represents the transformation of a matrix into a column vector, and H represents the complex common yoke transpose, subscript Represents taking the jth column of a frequency slice matrix, subscript Represents taking all columns of a frequency slice matrix. When estimating the source wavelet, the present invention does not require any assumptions, that is, it does not need to make any assumptions about the phase of the wavelet. In addition, compared with the traditional wavelet estimation problem that does not consider multiples, the multiple-wave term in the present invention helps to overcome the (amplitude and phase) non-uniqueness problems existing in the traditional wavelet estimation.

FIG. 2 is a flowchart of a preferred embodiment of the estimation method of the present invention, and the estimation method includes the following steps:

a) Obtain the following data: data d, threshold μ, maximum iteration update times k ^max of g ^est , internal iteration times l ^max of the accelerated linear Bregman algorithm, and maximum iteration update times m ^max of s ^est .

b) Initialization: the iteration variables k=0, m=0, b=d, ^{gest and sest are} initialized to 0.

c) Determine whether k<k ^max is established, if so, repeat step ci, otherwise execute step j.

d) According to x ₀ =g ^est , s ^est and b=d, call Algorithm 1 to update to obtain g ^est .

e) Update k to k+l ^max .

f) Determine whether m<m ^max is established, if so, execute step gi.

g) According to b and the updated g ^est , call formula 20 or 21 to update to obtain s ^est . Calculate each diagonal matrix by Equation 20 or 21 For each element on the diagonal in , according to the calculation of all Update s ^est .

h) Determine whether m=m ^max -1 is established, and if so, assign k ^max -k to l ^max .

i) Assign m+1 to m.

j) According to the final updated results, the estimation results of source wavelets and surface multiples are obtained.

There are two unknowns in the above process: the slice s ^est of the source wavelet and the slice g ^est of the Green's function. The process of solving the inverse problem is similar to the traditional blind deconvolution method, that is, fixing s ^est to invert g ^est ; Then fix g ^est to estimate s ^est . When inverting g ^est , each iteration update uses the previous inversion results. When inverting the sparse ^gest solution, it often needs several iterations, because according to the sparse recovery theory, it often takes several iterations to converge to obtain an optimal sparse solution. After a sparse solution g ^est is obtained, the source wavelet estimation is performed. The estimated s ^est is immediately used for inversion to update g ^est . The entire iterative main loop exits after satisfying a certain stopping criterion. The preset times of alternately updating g ^est and s ^est are determined by k ^max and m ^max . When g ^est reaches the maximum number of updates and k ^max but s ^est does not reach the maximum number of updates, after the last update of g ^est is performed, Executing the update of s ^est once completes the alternate update of g ^est and s ^est , in which the l ^max of the accelerated linear Bregman algorithm is equal each time; when s ^est reaches the maximum number of updates, k ^max and g ^est does not reach the maximum number of updates , after the last update of s ^est , only the last accelerated linear Bregman algorithm is executed, and the l ^max of the last accelerated linear Bregman algorithm is assigned as the difference between the maximum number of updates of g ^est and the number of times that have been updated, and before The ^lmax of each execution of the accelerated linear Bregman algorithm is equal. Among them, the lmax of the accelerated linear ^Bregman algorithm is equal each time. The last update g ^est is the estimated result of g, and the last update of s ^est is the estimated result of s. The following uses g ^est to represent the estimated result of g, and s ^est to represent s estimated results.

After the g ^est and s ^est are obtained, the primary wave p=g ^est *s ^est is obtained by calculation.

by g ^est and input data Obtain surface correlated multiples m=r _surf g ^est *d.

According to the principle of the secondary source, after obtaining the Green function g ^est and the obtained primary wave p, the time domain slices of the surface-related multiples of different orders are calculated by the feedback iterative model, for example, the first-order surface multiples m ₁ The calculation formula is:

m ₁ =r _surf g ^est *p=r _surf g ^est *g ^est *u ^est (22)

The _formula for calculating the slice m2 in the time domain of the second-order surface multiple is:

m ₂ =r _surf g ^est *m ₁ =(r _surf ) ² g ^est *g ^est *g ^est *s ^est (23)

The calculation formulas of surface-related multiples of other orders are analogous.

Since most of the calculation amount in multiple wave estimation is spent in each frequency slice doing matrix-matrix product, that is, calculating Equation 8, in order to further reduce the computational time consuming of multiple wave estimation, the present invention uses seismic data according to the (that is, the frequency components of seismic data are mainly concentrated in a frequency band of a limited length), the present invention only utilizes the advantages of high reliability and high signal-to-noise ratio of seismic data when estimating multiples and source wavelets frequency bands, which greatly reduces the computational time of the method, which is different from the traditional method which utilizes all (including low reliability and signal-to-noise ratio) frequency components.

After the primary wave or the surface multiple waves are obtained, the underground structure map can be obtained by using them respectively.

We provide the following specific examples to illustrate the feasibility, effectiveness and superiority of the technical solution of the present invention.

(1) Numerical example of sparse restoration to illustrate the superiority of the accelerated linear Bregman algorithm of the present invention.

(2) Example of salt dome model data. The data was simulated and generated by the DELPHI multiples research group at Delft University of Technology in the Netherlands, and it is one of the landmark data used for public testing of multiples algorithms in the world. It has appeared in van Groenestijn and Verschuur's 2009 paper, and Lin and Herrmann's 2013 paper and used to test multiple estimation methods.

(3) Pluto 1.5 data example. This data is simulated and released by the international Subsalt Multiples Attenuation and Reduction Team Joint Venture (SMAART JV) organization, and it is also one of the landmark data used for public testing of multiple wave algorithms in the world.

Numerical example of sparse recovery

The real model and simulated data in this numerical example are generated randomly. We first randomly generated a sparse vector of length 512 with 30 non-zero values in 512 points. For the 512 unknowns, we randomly generated 450 observational data, that is, for the 512 unknowns, 450 equations were established. This system of equations is underdetermined because the number of equations is less than the number of unknowns. For all participating methods, we set the maximum number of iterations to 20. Since the simulated data is noise free, the effects of noise on all algorithms can be isolated. right For the algorithm, we adopted the default parameters recommended by its developers. Here we use the signal-to-noise ratio (SNR) and the running time to evaluate the algorithms involved in the comparison. The calculation formula of the signal-to-noise ratio is:

Here x represents the true solution, represents the solution obtained.

From Figure 3 (the first row represents the real sparse model, that is, the true solution; the second row represents the result obtained by the least squares QR decomposition method; the third row represents the The result obtained by the algorithm solution; the fourth row represents the The error between the result obtained by the algorithm and the real solution; the fifth row represents the result obtained by the linear Bregman algorithm; the sixth row represents the result obtained by the accelerated linear Bregman algorithm of the present invention; the seventh row represents the acceleration of the present invention From the error between the result obtained by the linear Bregman algorithm and the real solution), we can see that there is a lot of random noise in the solution obtained by the least squares QR decomposition algorithm. The accuracy of the solution obtained by the algorithm is far inferior to that of the linear Bregman algorithm and the accelerated linear Bregman algorithm of the present invention. The error of the solution obtained by the algorithm is much larger than that of the accelerated linear Bregman algorithm of the present invention. The solution obtained by the accelerated linear Bregman algorithm of the present invention has the highest signal-to-noise ratio SNR. Fig. 3 also shows that the running time of the accelerated linear Bregman algorithm of the present invention is higher than The running time of the algorithm is much less (about 90% reduction).

From Figure 4 (the black thin solid line represents the change curve of the data residual of the least squares QR decomposition algorithm with the number of iterations; the black thin dashed line represents the The change curve of the data residual of the algorithm with the number of iterations; the black thick dashed line represents the change curve of the data residual of the linear Bregman algorithm with the number of iterations; the black thick solid line represents the data residual of the accelerated linear Bregman algorithm of the present invention. The variation curve of the number of iterations) shows that the data matching error changes with the number of iterations in the solution process, the convergence speed of the linear Bregman algorithm is better than that of the least squares QR decomposition algorithm and algorithm. However, the convergence speed of the accelerated linear Bregman algorithm of the present invention is further superior to that of the linear Bregman algorithm.

The numerical examples in Fig. 3 and Fig. 4 demonstrate the effectiveness and superiority of the accelerated linear Bregman algorithm in the present invention, that is, the algorithm is simple and the convergence rate is high.

Salt Dome Model Data Example

Figure 5a shows the data generated by the DELPHI multiple wave research group at Delft University of Technology in the Netherlands. The data is input by the velocity model and the density model, and is simulated by finite difference forward modeling with a time sampling interval of 4ms. The source wavelet used in the forward modeling is the Ricker wavelet, which is the same for every shot, that is, the constant wavelet. The shallow layer of the model is about 200 meters deep in water. In order to compare with Lin and Herrmann's REPSI method, the total number of iterations is set to 71 during processing. The results of the REPSI method are provided by Lin and Herrmann. The frequency band range used in the present invention is 0-70 Hz (the original whole frequency band range is 0-125 Hz).

From the processing results of the present invention in FIGS. 5 to 6 , the present invention has obtained better effects. In the original input data in Figure 6a, the multiple waves and the primary waves are mixed together, interfere with each other, and cannot be distinguished clearly. After being processed by the method of the present invention, the multiples and the primary are precisely separated. In the original data of Fig. 6a, the weak event axis belonging to the fault is difficult to identify, but after being processed by the present invention, it can be clearly identified and displayed in Fig. 6c. Since the present invention aims to deal with the multiples related to the surface, the interlayer multiples will be included in the primary wave obtained after processing, see Fig. 6c.

By comparing the results of the present invention (Fig. 6c) and the results of the REPSI method of Lin and Herrmann (Fig. 6d), it is not difficult to find that the primary wave processed by the present invention is more accurate, and the multiple waves can be better estimated and suppressed; Some multiple events are still visible in Herrmann's REPSI method results (Fig. 6d), as indicated by the arrows in Fig. 6d.

Comparing the source wavelet estimated by the present invention in Fig. 7a, the results estimated by Lin and Herrmann's REPSI method (the black solid line in Fig. 7(a) represents the source wavelet estimated by the present invention, and the dashed line represents the source wavelet estimated by Lin and Herrmann. The source wavelet estimated by the REPSI method), 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 constant wavelet that does not change with the shot collection is used in the simulation and generation of this data, the results generated when the variable wavelet is estimated in the present invention are relatively consistent, as shown in Fig. 7b.

From Figure 8 (the black solid line represents the inversion error of Lin and Herrmann's REPSI method with the number of iterations; the gray solid line represents the inversion error when the original linear Bregman algorithm is used (in the case of constant wavelet estimation) The change of the number of iterations; the black dotted line represents the change of the inversion error with the iteration number when the original linear Bregman algorithm is used (in the case of variable wavelet estimation); the gray dotted line represents the acceleration of the linear Bregman algorithm (in the case of variable wavelet estimation) ) The variation of inversion error with the number of iterations), the convergence speed of the method of the present invention is faster than that of the REPSI method of Lin and Herrmann from the point of view of the variation of the data matching error in the iterative process. The convergence speed of constant wavelet estimation is basically the same as that of variable wavelet estimation. The convergence speed of the accelerated linear Bregman algorithm of the present invention is better than the convergence speed of the original linear Bregman algorithm.

The running time data in the following table shows that: under the same computing environment and computing resources, the running time of the method of the present invention is far less than the running time of the REPSI method of Lin and Herrmann. The calculation time of variable wavelet estimation is slightly more than that of constant wavelet estimation. The main reason why the running time of the method of the present invention is far less than the running time of the REPSI method of Lin and Herrmann is as follows:

(1) The present invention accelerates the simplicity of the steps of the linear Bregman algorithm, which is the same as that adopted by Lin and Herrmann. The complexity of the algorithms contrasts sharply;

(2) The present invention only needs to calculate the simple and efficient step of formula 20 or 21 when estimating the source wavelet; while the REPSI method of Lin and Herrmann calls the least squares QR decomposition algorithm when estimating the source wavelet, which requires multiple iterations;

(3) According to the band-limited characteristics of seismic data, the present invention only utilizes the dominant frequency band 0-70 Hz (see Fig. 6b), which has high reliability of seismic data and high signal-to-noise ratio, in the estimation of multiples and source wavelets (see Fig. 6b); And Herrmann's method uses the entire frequency band 0-125Hz, resulting in a lot of unnecessary computational time.

In summary, we use the same data and the same computing resource environment to compare with the REPSI method proposed by Lin and Herrmann in 2013, and demonstrate the advantages of the present invention.

method operation hours Lin and Herrmann's REPSI method 21 minutes 6 seconds The method of the present invention (constant wavelet estimation case) 5 minutes 57 seconds Method of the present invention (variable wavelet estimation case) 6 minutes 4 seconds

Pluto 1.5 data example

To further demonstrate the effectiveness of the present invention, we apply it to another convincing Pluto1.5 data. This data is simulated and released by the international SMAART JV organization, and it is one of the landmark data used for public testing of multiple wave algorithms in the world. The source wavelet used in the forward modeling of the data is the Ricker wavelet, which is the same for each shot, that is, the constant wavelet. The time sampling interval is 8 ms. We directly used the raw data publicly released by the SMAART JV tissue, and did not do any preprocessing of the raw data except for the removal of direct waves. From the raw input data in Figures 6a and 7a, there is a large amount of multiple energy in this data. The frequency band used in the process of the present invention is 0-50 Hz.

From the processing results in Fig. 9 and Fig. 10, in the original input data Fig. 9a and Fig. 10a, the multiple waves and the same primary wave are mixed together, interfere with each other, and are indistinguishable, while the present invention realizes accurate primary wave, Multiple estimation and separation (especially where the arrows in Figure 7 point to). Since the source wavelet used in the generation of the data simulation does not change with the shot collection, the results generated when the variable wavelet estimation is performed in the present invention are relatively consistent, as shown in Fig. 11 .

Figure 12 shows another scheme for evaluating the effectiveness of multiple estimation techniques. The essence of Fig. 12 is a feedback iterative model, that is: the primary wave (Fig. 12b) estimated by the present invention is introduced into the ground as a secondary source to predict the first-order surface-related multiples (Fig. 12c), and the first-order surface-related multiples (Fig. 12c). (Fig. 12c) As a secondary source, the second-order surface correlated multiples (Fig. 12d) were introduced underground to predict the second-order surface-related multiples (Fig. 12d), and the second-order surface-related multiples (Fig. 12d) were introduced underground as a secondary source to predict the third-order surface correlated multiples. Secondary waves (Fig. 12e), third-order surface-related multiples (Fig. 12e) are transmitted underground as secondary sources to predict fourth-order surface-related multiples (Fig. 12f), and fourth-order surface-related multiples (Fig. 12f) The fifth-order surface-related multiples (Fig. 12g) are predicted to be transmitted underground as a secondary source, and so on. For this data, the predicted surface-correlated multiples of different orders (Figs. 12c to 12g) can be easily identified in the raw input data (Fig. 12a). The total data (Fig. 12h) generated by the cumulative summation of the primary wave and surface-related multiples of different orders estimated by the present invention is almost identical to the original input data (Fig. 12a), and the difference between the two (Fig. 12a) 12i) is very small.

Fig. 13 further shows the processing results of the present invention applied to other shot sets in the Pluto 1.5 data. It can be seen that the primary and multiple waves are more accurately estimated and separated.

The embodiments of the present invention have been described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific embodiments, which are merely illustrative rather than restrictive. Under the inspiration of the present invention, without departing from the scope of protection of the present invention and the claims, many forms can be made, which all belong to the protection of the present invention.

Claims (6)

1. a surface multiple and wavelet estimation method based on accelerated linear Bregman algorithm, is characterized in that, comprises following steps: S1, use the following steps S11 and S12 to alternately update the Green's function to be estimated and the source wavelet to be estimated until reaching a preset number of times, wherein the source wavelet to be estimated that performs step S11 for the first time has a preset initial value: S11. According to the source wavelet to be estimated and the collected input data, use the accelerated linear Bregman algorithm to update the Green's function to be estimated; S12, according to the updated Green's function to be estimated and the input data, adopt the following formula to update the source wavelet to be estimated; When the source wavelet to be estimated does not change with the spatial position of the shot collection: When the source wavelet to be estimated changes with the spatial position of the shot collection: S2. According to the final updated result, obtain the estimation result of at least one of the source wavelet and the surface multiple waves; Among them, the superscript est means to be estimated, Represents the elements on the diagonal of the source wavelet matrix slice, j represents the spatial position of the shot set or the shot set number, vec() represents the transformation of a matrix into a column vector, and the subscript Represents taking the jth column of a frequency slice matrix, subscript represents taking all columns of a frequency slice matrix, is a slice of a matrix expression in the frequency domain of the input data, is the slice of the matrix expression in the frequency domain of the source wavelet, r _surf is the reflection coefficient of the contact surface between air and water, is a slice of the matrix expression in the frequency domain of the Green's function; According to the source wavelet to be estimated and the collected input data, using the accelerated linear Bregman algorithm to update the Green's function to be estimated specifically includes: According to the formula and Iterate to find y ₁ , x ₁ , y ₂ , x ₂ , y ₃ …, in, is the vector representation in the time domain of the update result of the Green's function to be estimated after step S11 is executed this time, and ^lmax is the maximum number of iterations of the accelerated linear Bregman algorithm when step S11 is executed this time; In the formula, x ₀ , and y ₀ are both equal to the vector representation in the time domain of the preset initial value of the Green's function to be estimated, and the observation data vector b is the vector representation in the time domain of the input data, t _l represents the step size, and its expression is: The mapping function ∏ _σ (Ax _l -b) is defined as: The soft threshold function shrink(x, μ) is defined as: shrink(x,μ)=max(|x|-μ,0)sign(x), Among them, σ represents the noise level factor, the max() function represents the maximum value, the || function represents the absolute value, sign(x) represents the sign value function, and μ represents the threshold; The matrix operator A represents a set of operations to perform the following operations: Fourier transform the Green's function to be estimated from a vector table in the time domain to a vector table in the frequency domain, and convert the Green's function to be estimated from the frequency domain. The vector table form becomes the matrix table form in the frequency domain, according to the formula Calculate each slice of the matrix representation of the frequency domain of the analog data corresponding to the input data separately According to each slice Inverse Fourier transform is performed to obtain the vector representation in the time domain of the analog data corresponding to the input data. 2. surface multiples and wavelet estimation method according to claim 1, is characterized in that, in described 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 multiples is obtained according to m=r _surf g*d, where m is the vector representation in the time domain of the surface multiples, g is the vector representation in the time domain of the Green’s function, and d is the time domain of the input data vector representation of . 3. surface multiples and wavelet estimation method according to claim 1, is characterized in that, if the Green's function to be estimated has reached the maximum number of updates and the source wavelet to be estimated does not reach the maximum number of updates, then step again In S11, ^lmax is updated to the difference between the maximum update times of the source wavelet to be estimated minus the number of times that have been updated, and ^lmax is the same each time step S11 is performed before; if the source wavelet to be estimated has been updated When the maximum number of updates is reached but the Green's function to be estimated does not reach the maximum number of updates, then step S12 is performed only once, and ^lmax is the same each time step S11 is performed. 4. a surface multiple and wavelet estimation system based on accelerated linear Bregman algorithm, is characterized in that, comprises: The alternate update unit is used to alternately update the Green's function to be estimated and the source wavelet to be estimated by using the following Green's function update unit and the source wavelet update unit until a preset number of times is reached, wherein the Green's function update unit to be estimated is called for the first time. The hypocenter wavelet has preset initial values: The Green's function updating unit adopts the accelerated linear Bregman algorithm to update the Green's function to be estimated according to the source wavelet to be estimated and the collected input data; The source wavelet updating unit adopts the following formula to update the source wavelet to be estimated according to the updated Green's function to be estimated and the input data; When the source wavelet to be estimated does not change with the spatial position of the shot collection: When the source wavelet to be estimated changes with the spatial position of the shot collection: an estimation result obtaining unit, configured to obtain an estimation result of at least one of the source wavelet to be estimated and the surface multiple to be estimated according to the final updated result; Among them, the superscript est means to be estimated, Represents the elements on the diagonal of the source wavelet matrix slice, j represents the spatial position of the shot set or the shot set number, vec() represents the transformation of a matrix into a column vector, and the subscript Represents taking the jth column of a frequency slice matrix, subscript represents taking all columns of a frequency slice matrix, is a slice of a matrix expression in the frequency domain of the input data, is the slice of the matrix expression in the frequency domain of the source wavelet, r _surf is the reflection coefficient of the contact surface between air and water, is a slice of the matrix expression in the frequency domain of the Green's function; According to the source wavelet to be estimated and the collected input data, using the accelerated linear Bregman algorithm to update the Green's function to be estimated specifically includes: According to the formula and Iterate to find y ₁ , x ₁ , y ₂ , x ₂ , y ₃ …, in, is the vector representation in the time domain of the update result of the Green's function to be estimated after the Green's function update unit is executed this time, and ^lmax is the maximum number of iterations of the accelerated linear Bregman algorithm when the Green's function update unit is executed this time; In the formula, x ₀ , and y ₀ are both equal to the vector representation in the time domain of the preset initial value of the Green's function to be estimated, and the observation data vector b is the vector representation in the time domain of the input data, t _l represents the step size, and its expression is: The mapping function ∏ _σ (Ax _l -b) is defined as: The soft threshold function shrink(x, μ) is defined as: shrink(x,μ)=max(|x|-μ,0)sign(x), Among them, σ represents the noise level factor, the max() function represents the maximum value, the || function represents the absolute value, sign(x) represents the sign value function, and μ represents the threshold; The matrix operator A represents a set of operations to perform the following operations: Fourier transform the Green's function to be estimated from a vector table in the time domain to a vector table in the frequency domain, and convert the Green's function to be estimated from the frequency domain. The vector table form becomes the matrix table form in the frequency domain, according to the formula Calculate each slice of the matrix representation of the frequency domain of the analog data corresponding to the input data separately According to each slice Inverse Fourier transform is performed to obtain the vector representation in the time domain of the analog data corresponding to the input data. 5. The surface multiple and wavelet estimation system according to claim 4, wherein, in the estimation result obtaining unit, 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 multiples is obtained according to m=r _surf g*d, where m is the vector representation in the time domain of the surface multiples, g is the vector representation in the time domain of the Green’s function, and d is the time domain of the input data vector representation of . 6. The surface multiples and wavelet estimation system according to claim 4, is characterized in that, 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 Green's function to be estimated has not reached the maximum number of updates. When the function updates the unit, ^lmax is updated to the difference between the maximum update times of the source wavelet to be estimated minus the number of times that have been updated, and the ^lmax when the Green’s function update unit is executed before is the same; If the source wavelet has reached the maximum number of updates but the Green's function to be estimated has not reached the maximum number of updates, then the source wavelet update unit is performed only once, and the ^lmax is equal each time the Green's function update unit is executed.