Marchenko Multiple Wave Suppression Method for Terrestrial Measured Data Based on Compressive Sensing

Abstract

The present invention belongs to the technical field of geophysical exploration, and relates to a Marchenko multiple suppression method for onshore measured data based on compressive sensing, including a data preprocessing link: original data encryption and regularization; Shearlet denoising and reconstruction; sparse constraint deconvolution; an iterative link for suppressing multiples: initializing the coda for the first iteration; obtaining the coda for odd and even iterations; and performing signal-to-noise separation based on Bayes' principle in even iterations to obtain the sum of the coda terms for even iterations that suppress noise. By means of an optimized preprocessing process, the present invention obtains higher-quality input data and performs a more noise-resistant Marchenko iterative process, effectively suppressing the multiple wave signals in onshore measured data without the need for any model information and predictive subtraction, laying a good foundation for subsequent imaging and interpretation, improving the accuracy of oil and gas resource exploration, and saving the costs consumed by more manual judgment of false horizons.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical exploration, and particularly relates to suppression of interlayer multiple noise applicable to onshore measured seismic data. Background Art

[0002] In energy exploration, the quality of seismic data directly affects the accuracy of reservoir identification and the efficiency of oil and gas reservoir development. Seismic data is obtained by artificially exciting seismic waves and recording their propagation characteristics (such as reflection, refraction, and attenuation) in underground media. The primary wave in seismic data is the effective signal that the seismic wave starts from the source, reflects once from the underground interface, and then returns to the ground. It can truly reflect information such as the depth and structural form of the formation interface. Different from it, the multiple wave is an interference wave formed by multiple reflections between underground wave impedance interfaces. Usually, it will cause false structural information to appear in the seismic imaging result and interfere with resource positioning. Therefore, it is also regarded as a kind of coherent noise in the seismic wave field. Multiple waves can be divided into surface multiples and interlayer multiples according to their formation mechanisms. Compared with the surface multiples with strong energy and periodicity, the interlayer multiples have non-periodic characteristics, and the aliasing phenomenon with the primary wave is more serious. Therefore, the suppression of interlayer multiples has become a key technical problem to be solved urgently. An effective method for suppressing interlayer multiples is of great significance for improving the quality of seismic data, improving the subsequent imaging effect, realizing accurate reservoir identification and interpretation, and accurately positioning oil and gas resources.

[0003] Conventional multiple suppression methods, such as the surface-related multiple elimination method (SRME), have been maturely applied. However, this method is usually applicable to the suppression of surface multiples in marine data. With the increasing demand for resource exploration, the present invention finds that the development of interlayer multiples is very serious in many land areas, and the suppression of interlayer multiples has become a key problem to be solved. Conventional methods for suppressing interlayer multiples, such as the inverse scattering series interlayer multiple suppression method (ISS) and the interlayer multiple suppression method based on the focusing transform (CFP), can successfully suppress interlayer multiples to a certain extent. However, these conventional interlayer multiple suppression methods and the above-mentioned SRME method all rely on the strategy of predicting multiples and matching and subtracting. Such a prediction and subtraction strategy is based on some characteristic differences between the primary wave and the multiple wave to predict the multiple wave (such as amplitude, waveform, phase, etc.). However, in actual data, especially in mid-deep land data, the interlayer multiples are usually implicitly developed and the aliasing phenomenon between the interlayer multiples and the primary wave is very serious. If following such a strategy, subtracting after predicting the multiple wave is very likely to damage the effective signal of the primary wave.

[0004] The Marchenko method directly separates the primary wave and the multiple wave by constructing the Green's function. Therefore, prediction subtraction is not required. In the Marchenko multiple wave elimination method (MME) in the Marchenko method, on the basis of the traditional Marchenko method that relies on the velocity model to obtain the initial extrapolation operator to perform the extrapolation to the subsurface on one side of the shot point or the receiver point, a convolution of the extrapolation operator is applied to both sides of the extrapolated Green's function expression, and the shot point and the receiver point are projected back to the surface. In this forward process and inverse process, the effects of the extrapolation operator estimated by the velocity model cancel each other out, and finally it becomes a delta function that only depends on time. To sum up, the Marchenko multiple wave elimination method neither requires the input of the velocity model nor prediction subtraction, and has a more solid theoretical support in suppressing the interlayer multiple waves.

[0005] However, the Marchenko multiple wave elimination method needs to follow three assumptions:

[0006] (1) There is no noise interference in the data;

[0007] (2) The data needs to meet the ultra-high density sampling condition;

[0008] (3) The wavelet is approximately a pulse sequence, and the reflection response converges to a sharp pulse signal.

[0009] However, the noise development mechanism of the measured onshore data is very complex, and the surface wave noise coverage is usually relatively serious. Considering the amplitude preservation of the effective signal, there will still be some residual noise in the data after conventional preprocessing. And due to actual terrain and other factors, even under high-density sampling conditions, it is very difficult for the actual onshore data obtained to meet the ultra-high density sampling assumption of the Marchenko multiple wave elimination, and the wavelet in the actual onshore data often is not a pulse sequence either.

[0010] Based on this, there is an urgent need to develop a new Marchenko interlayer multiple wave suppression method for the measured onshore seismic data to effectively solve the above problems. Summary of the Invention

[0011] The object of the present invention is to provide a compressed sensing improved Marchenko multiple wave suppression method applicable to onshore actual seismic data, so as to solve the problems of improving the suppression accuracy of the Marchenko method for the interlayer multiple waves in onshore actual data and the accuracy of the final seismic data interpretation, and provide a good data basis for subsequent reservoir identification and resource exploration.

[0012] The object of the present invention is achieved by the following technical solutions:

[0013] A Marchenko multiple suppression method for land measured data based on compressive sensing, comprising the following steps:

[0014] Step 1, extract three-dimensional seismic data and perform regularization processing to obtain original data ;

[0015] Step 2, according to the compressive sensing theory, perform Shearlet transform on the original data obtained in Step 1 :

[0016] (1)

[0017] Wherein, represents the Shearlet transform, is the Shearlet coefficient after transformation;

[0018] First, construct the Shearlet basis function :

[0019] (2)

[0020] Wherein, is the scale matrix, is the shear matrix, is the scale parameter, is the shear parameter, is the translation parameter, and the scale matrix and the shear matrix are defined in the following form:

[0021] (3)

[0022] Then, perform Shearlet transform on the original data obtained in Step 1 :

[0023] (4)

[0024] Step 3, perform denoising processing in the Shearlet domain, and obtain the denoised data from formula (5):

[0025] (5)

[0026] Wherein, is the original seismic data to be denoised observed, is the sparse representation matrix, represents the first norm of, represents the error threshold, is the solution to solve the above denoising problem, represents the solution estimate value, represents the denoised data;

[0027] Step 4: Reconstruct in the Shearlet domain to obtain the input seismic record for the multiple removal section. Specifically, interpolate the missing parts of the denoised data to obtain the high-density sampled seismic data without missing parts, and perform sparse reconstruction of the seismic data in the Shearlet domain, so as to transform the problem of interpolative reconstruction of the missing seismic traces into the following L1 norm minimization problem, and obtain the reconstructed data from Equation (7): wherein, the missing parts of the denoised data are interpolated to obtain the high-density sampled seismic data without missing parts, and sparse reconstruction of the seismic data is performed in the Shearlet domain, so as to transform the problem of interpolative reconstruction of the missing seismic traces into the following L1 norm minimization problem, and obtain the reconstructed data from Equation (7):

[0028] (7)

[0029] wherein, represents the inverse Shearlet operator, is the synthesis interpolation operator, is the Lagrange multiplier, represents the target eigenvector, represents the regularization solution, represents the reconstructed data;

[0030] Step 5: Perform sparse constraint deconvolution on the reconstructed data obtained in Step 4 to obtain the input reflection response. Specifically, optimize the wavelet estimation to meet the assumptions of Marchenko multiple elimination while obtaining the high-resolution impulse response for the input of Marchenko multiple removal , and the sparse constraint deconvolution model is expressed as:

[0031] (8)

[0032] wherein, represents the seismic trace, is the wavelet convolution matrix, is the reflection coefficient sequence, represents the noise term;

[0033] Solve Equation (9) using the iterative soft threshold algorithm in the Shearlet domain, and obtain the final reflection coefficient in the Shearlet domain obtained by the action of the inverse Shearlet transform operator The corresponding reflection coefficient in the time domain is:

[0034] (9)

[0035] wherein, represents the obtained reflection coefficient in the Shearlet domain, is the wavelet operator, represents the noise level;

[0036] Step 6: Use the input seismic records and input reflection responses obtained in Steps 4 and 5 to perform Marchenko multiple wave removal, and then solve for the wavelet. The final multiple wave suppression form after solving the wavelet is expressed as:

[0037] (10)

[0038] Where, is the original reflection response containing multiple waves, is the sum of the even terms of the solved wavelet, representing an inverse scattering term with the same energy as the multiple waves but opposite in polarity. After and are added together, the multiple waves in are suppressed to obtain the reflection response without multiple waves ; ;

[0039] Step 7: Solve for the coda wave through iterative calculation , where the subscript represents the number of iterative updates. Initialize the first iteration of the coda wave:

[0040] (17)

[0041] Where, represents time, represents the shortest two-way travel time of the first arrival wave, represents a very small smooth window function quantity, represents the current receiver position, represents the current shot point position;

[0042] Step 8: Odd-numbered iterations. The odd-numbered coda wave iterations follow the following formula:

[0043] (18)

[0044] Where, represents the receiver coordinate when the source is excited at the position, is the integration time variable;

[0045] Step 9: Even-numbered iterations. The coda wave calculation follows the following formula. In the iterative solution process of the even-numbered terms, a signal-to-noise separation algorithm based on the Bayesian principle is introduced to obtain the sum of the even-numbered coda waves ;

[0046] (19)

[0047] Where, represents the integration time variable;

[0048] Step 10: The sum of the even-numbered coda waves obtained in Step 9 Substitute into Step 6 to finally obtain data without multiple reflections between layers .

[0049] Further, Step 1 is specifically: extract the originally stored data of time×number of channels into 3D seismic data of time×number of channels×number of shots, and regularize the original seismic data. Regularize the data with originally relatively large shot intervals and geophone intervals into high sampling density, that is, sparse data with relatively small shot intervals and geophone intervals. The originally existent channels with valid values are regularized to corresponding positions, and positions without data are assigned zero values to obtain the original data .

[0050] Further, in Step 3, solve the seismic data without noise through the hard threshold algorithm, specifically:[[]]

[0051] (6)[[]]

[0052] Among them,[[]] represents the independent variable of the function substituted into the hard threshold algorithm,[[]] represents the decision threshold,[[]] represents the hard threshold function.[[]]

[0053] Further, in Step 4, adopt the iterative threshold cooling strategy to determine the Lagrange multiplier[[]] .

[0054] Further, in Step 6, the Marchenko multiple reflection removal step is: convolve both sides of the Green's function expressions based on the autofocus theory in Formula (11) and Formula (12) with the continuation operator[[]] , that is, the downward Green's function of the direct wave, to obtain Formula (13) and Formula (14):[[]]

[0055] (11)[[]]

[0056] (12)[[]]

[0057] (13)[[]]

[0058] (14)[[]]

[0059] Among them,[[]] represents the spatial position coordinates, and its subscript[[]] and[[]] respectively represent the depths of the spatial positions,[[]] represents the spatial position at a depth of 0,[[]] represents the spatial position at a depth of[[]] , and the subscript[[]] It is used to distinguish the spatial positions in the marking formula (11) and formula (12) separately and convert them into formula (13) and formula (14) after convolution. spatial location; Represents a plane with a depth of 0;

[0060] 5 to obtain the input reflection response, where Indicates the location of the earthquake source, Indicates the receiver location, both parts are in superior;

[0061] and The earthquake source is located at , the detection point is located at The descending and ascending Green's functions, and represents the downlink and uplink focusing functions;

[0062] represents the Green's function and and The convolution result of , where the superscripts “+” and “-” represent the downgoing wavefield and the upgoing wavefield respectively; for The convolution form of the coda component, for The convolution form of

[0063] Representation space is limited to two dimensions function, and Representing the time domain function; definition for First arrival time, i.e., the shortest two-way travel time from the surface to the focal plane and back; and The time interval is ; In this range , It can be determined by solving formula (13) and formula (14);

[0064] Will Substituting into formula (13) we can get The final expression of:

[0065] (15)

[0066] Will Defined as , collect each Corresponding values, and represent the final data without multiples in the form of scattered coda waves as:

[0067] (16)

[0068] where is the original input reflection response, is the reflection response after removing multiples, represents the even coda waves, represents a signal with the same amplitude as the interlayer multiples but opposite polarity.

[0069] Furthermore, step 9 is specifically: when solving for the even coda waves, first predict the effective terms and noise terms in the input coda waves in the Shearlet domain, and then obtain the Shearlet coefficients and that maximize the Bayesian conditional probability, and then perform Bayesian subtraction on the noise terms:

[0070] (20)

[0071] The Shearlet coefficients and that maximize the Bayesian conditional probability are obtained by solving formula (21):

[0072] (21).

[0073] Compared with the prior art, the beneficial effects of the present invention are:

[0074] 1. The present invention provides an improved Marchenko interlayer multiple suppression method for compressive sensing of onshore actual data, which specifically improves the limitations of the conventional Marchenko multiple elimination method that leads to poor final results in onshore seismic data and improves the conventional Marchenko iterative algorithm, and finally forms a complete multiple suppression system suitable for onshore seismic data. The improved multiple suppression system includes a data and processing link and a Marchenko multiple suppression link. In the data preprocessing link, the present invention improves the commonly used industrial filtering method to Shearlet denoising with multi-angle, multi-direction, and multi-scale recognition, which can denoise more accurately.

[0075] 2. The present invention improves the commonly used simple spline interpolation method in industry to Shearlet reconstruction to obtain a more continuous and smooth reconstruction result. The present invention upgrades the conventional pulse deconvolution method to sparse constraint deconvolution, which can obtain a more high-resolution and noise-resistant reflection response, making the suppression of multiple waves more accurate. Aiming at the noise that may be amplified during the iteration process and remains due to the consideration of amplitude preservation, the present invention improves the traditional Marchenko iteration to an iterative process with high noise resistance based on Bayesian principle for Shearlet noise separation, and effectively suppresses the residual noise and the possible new noise generated during each iteration process.

[0076] 3. The present invention is applicable to the suppression of interlayer multiple waves in onshore measured data, and successfully breaks through the application limitations that are difficult to overcome by the Marchenko multiple wave removal method in actual onshore data. By optimizing the preprocessing process to obtain higher-quality input data that conforms to the assumptions and performing a more noise-resistant Marchenko iterative suppression process, the interlayer multiple wave signals in the measured onshore data can finally be effectively and accurately suppressed, and no model information and predictive subtraction are required, laying a good data foundation for subsequent migration imaging and interpretation, improving the accuracy of oil and gas resource exploration, and saving the cost consumed by more manual judgment of false horizons. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following will briefly introduce the drawings required in the embodiments. It should be understood that the following drawings only show some embodiments of the present invention, and therefore should not be regarded as limiting the scope. For those of ordinary skill in the art, without creative efforts, other related drawings can also be obtained based on these drawings.

[0078] Figure 1 Flowchart of the steps of the Marchenko multiple wave suppression method for onshore measured data based on compressive sensing;

[0079] Figure 2 is the original measured onshore data;

[0080] Figure 3 is the display of the Shearlet denoising effect; among them, (a) is the original data input to the denoising process, (b) is the data after conventional filtering denoising, (c) is the data after Shearlet denoising, and (d) is the noise signal removed by Shearlet denoising compared to the conventional method;

[0081] Figure 4It shows the Shearlet reconstruction effect after denoising. Among them, (a) is the original data before denoising and reconstruction and its local enlargement, (b) is the data after conventional filtering denoising and conventional spline interpolation and its local enlargement, and (c) is the data after improved Shearlet denoising and reconstruction and its local enlargement;

[0082] Figure 5 It shows the effect of sparse constraint deconvolution. Among them, (a) is the reflection response obtained by the conventional pulse deconvolution method and its local enlargement, and (b) is the reflection response obtained by the sparse constraint deconvolution method and its local enlargement;

[0083] Figure 6 It shows the effect after processing by the improved compressed sensing Marchenko interlayer multiple suppression system of the present invention; (a) is the original data before suppression and its local enlargement, (b) is the result of conventional Marchenko multiple suppression, and (c) is the effect after the improved compressed sensing Marchenko interlayer multiple suppression of the present invention;

[0084] Figure 7 It is a common offset gather. Among them, (a) is the common offset gather before suppression, and (b) is the common offset gather after suppression by the system of the present invention;

[0085] Figure 8 It is an autocorrelation spectrum. Among them, (a) is the autocorrelation spectrum before suppression, (b) is the long-period multiple of the autocorrelation spectrum after suppression by the system of the present invention, (c) is before noise suppression, and (d) is after noise suppression. Detailed implementation manners

[0086] The present invention will be further described below in conjunction with embodiments:

[0087] The present invention will be further described in detail below in conjunction with the drawings and embodiments. It can be understood that the specific embodiments described herein are only used to explain the present invention, rather than limiting the present invention. In addition, it should be noted that for the sake of description, only the parts related to the present invention are shown in the drawings, rather than all the structures.

[0088] It should be noted that: similar reference numerals and letters represent similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are only used for distinguishing descriptions, and cannot be understood as indicating or implying relative importance.

[0089] The existing Marchenko multiple elimination methods have problems such as complex noise in land measured data, undersampling in onshore actual data, and the wavelets in land measured data usually not being the same for the entire shot gather and usually not being close to a sharp pulse signal. Therefore, the present invention provides an improved Marchenko multiple suppression method based on compressive sensing applicable to onshore actual seismic data, including an improved data preprocessing and an improved Marchenko multiple suppression link. Different from the conventional Marchenko multiple elimination methods restricted by the complex conditions of onshore actual data, the present invention designs specific solutions for each restrictive problem existing in onshore actual data, and finally forms a complete Marchenko multiple suppression method for onshore actual data.

[0090] Specifically, in the preprocessing link: for the complex noise in land measured data, the present invention improves the conventional industrial filtering process, performs Shearlet transform on the data in the sparse domain, and performs Shearlet denoising in the Shearlet domain with multi-angle, multi-scale, and multi-direction recognition, so that the complexly developed noise can be more accurately identified and removed, and the amplitude of the primary wave can be better protected. For the undersampling problem in onshore actual data, the present invention upgrades the conventional industrial spline interpolation method to Shearlet wavefield reconstruction. Compared with the conventional method, Shearlet wavefield reconstruction can better improve the continuity of the in-phase axis after reconstructing the undersampled data. Moreover, Shearlet wavefield reconstruction under the compressive sensing system can represent the sparse sampling signal in a certain sparse transform domain, and the ideal signal can be reconstructed from the non-uniform undersampled data by solving the optimization problem, which means that the sampling cost can be greatly saved. For the problem that the wavelets in land measured data usually are not the same for the entire shot gather and usually not close to a sharp pulse signal, the present invention replaces the commonly used pulse deconvolution method in industry with a deconvolution method with sparse constraints. The reflection response after convolving the wavelet obtained by this processing with the observed data converges more to a sharp pulse sequence, and such an input reflection response can better meet the iterative convergence conditions of Marchenko.

[0091] In the link of Marchenko multiple elimination: The present invention improves the iterative algorithm for conventional Marchenko multiple elimination. In the conventional Marchenko multiple removal process, considering the amplitude preservation of the primary wave, there are usually some residual noises after preprocessing. These residual noises will be enhanced in the subsequent iterative process, ultimately affecting the suppression effect of the interlayer multiples. To solve this problem, the present invention improves the original iterative process into an iterative process with high anti-noise ability. Specifically, in the process of updating the even-numbered coda waves in each iteration, the Shearlet signal-to-noise separation algorithm based on the Bayesian principle is added. First, the noise term and the effective signal term in the even-numbered coda waves are predicted, and then the noise term is subtracted by Bayesian, successfully suppressing the residual noise and the possible new noise effectively in each iterative process. Through the above preprocessing link and iterative link, corresponding solutions are designed for each problem existing in the land data, and finally a complete Marchenko interlayer multiple suppression method for land data based on the principle of compressive sensing is formed, improving the suppression accuracy of the Marchenko method for the interlayer multiples in the onshore actual data and enhancing the accuracy of the final seismic data interpretation, providing a good data basis for subsequent reservoir identification and resource exploration.

[0092] Example 1:

[0093] A Marchenko multiple suppression method for onshore measured data based on compressive sensing has been applied to a piece of onshore measured data for effect analysis. The original data sampling interval is 8 ms, the total sampling time is 6 s, the original shot interval is 25 m, and the original trace interval is 25 m. In the preprocessing link, both the shot interval and the trace interval are encrypted to 12.5 m.

[0094] As Figure 1 shown, it includes the following steps:

[0095] Data preprocessing link:

[0096] Step 1, establishment of the acquisition system. Extract the data stored in the original time * number of traces into 3D seismic data of time * number of traces * number of shots. At the same time, regularize the original seismic data, regularize the data with relatively large shot intervals and geophone intervals into sparse data with high sampling density (relatively small shot intervals and geophone intervals). The traces with original valid values are regularized to the corresponding positions, and the positions without data are assigned zero values to obtain the original data , and are encrypted into high-sampling-density data in the subsequent interpolation step. The original data processed by this step is shown in Figure 2, it can be seen that there is relatively strong interfering surface wave noise in the original data after preliminary processing. Since the data has not been reconstructed and encrypted yet, it presents an undersampled sparse form. The signal is not close to a sharp pulse and there is some high-energy noise interference before 0.4 s, and these interferences are subsequently removed.

[0097] Step 2: According to the compressive sensing theory, perform the following operations on the original data obtained in Step 1 for Shearlet transform:

[0098] (1)

[0099] where represents the Shearlet transform, is the Shearlet coefficient after transformation, which is used for subsequent denoising and other processing.

[0100] Specifically: First, construct the Shearlet basis function :

[0101] (2)

[0102] where is the scale matrix, is the shear matrix, is the scale parameter, is the shear parameter, is the translation parameter. The scale matrix and the shear matrix are defined in the following form:

[0103] (3)

[0104] Then perform the Shearlet transform on the original data obtained in Step 1 :

[0105] (4)

[0106] Step 3: Perform denoising processing in the Shearlet domain. Since the distribution of noise in the Shearlet domain is different, transform the denoising problem into solving the noise-free data:

[0107] (5)

[0108] where is the observed original seismic data to be denoised, is the sparse representation matrix, denotes the l1-norm of, denotes the error threshold, is the solution to solve the above denoising problem, denotes the solution estimate value, denotes the data after denoising. The ultimate goal of this denoising problem is to solve the noise-free data. Specifically, the computationally expensive iterative solution method can be replaced by the hard thresholding algorithm. The development of the thresholding algorithm is relatively mature. The specific formula can be seen as follows:

[0109] (6)

[0110] Among them, denotes the independent variable of the function substituted into the hard thresholding algorithm, denotes the decision threshold, denotes the hard thresholding function.

[0111] In this experiment, the results of Shearlet denoising can be specifically seen in Figure 3 , and through comparative analysis Figure 3 (c), it can be seen that the data processed through this step has successfully removed most of the surface wave noise interference compared to the original input data Figure 3 (a). Compared with the results after conventional filtering denoising Figure 3 (b), it has removed more steeply dipping noise represented by the white arrows, Figure 3 (d) denotes the noise signal removed by Shearlet denoising compared to the conventional filtering denoising method.

[0112] Step 4: Reconstruct in the Shearlet domain to obtain the input seismic record for the multiple removal link; specifically, interpolate the missing part of the denoised data . The high-density sampled seismic data without missing parts obtained through this step is used as the input seismic record (for constructing the Green's function) that meets the assumptions of the Marchenko multiple removal link. The problem of interpolating and reconstructing missing seismic traces by performing seismic data sparse reconstruction in the Shearlet domain can be transformed into the following L1 norm minimization problem:

[0113] (7)

[0114] Among them, denotes the Shearlet inverse operator, is the synthesis interpolation operator, is the Lagrange multiplier, denotes the target eigenvector, denotes the regularized solution, denotes the data after reconstruction. The purpose of this denoising problem is to minimize the Shearlet coefficients. Therefore, the present invention adopts an iterative threshold cooling strategy to determine the Lagrange multiplier .

[0115] In the present invention, the Shearlet reconstruction is performed on the basis of the previous Shearlet denoising. Figure 4 The comparison of the effects after denoising and reconstruction is shown. Figure 4 (a) is the original data before denoising and reconstruction. Figure 4 (b) is the data after conventional filtering denoising and simple spline interpolation processing. Figure 4 (c) is the data after improved Shearlet denoising and reconstruction processing. Comparing Figure 4 (a), Figure 4 (b) and Figure 4 (c), it can be seen that Figure 4 the steep dip angle noise within the black dotted parallelogram frame in (b) is Figure 4 eliminated in (c). Compared with Figure 4 (b), Figure 4 the event axes at many positions indicated by the white arrows in (c) are significantly more natural in transition, and Figure 4 many discontinuous weak event axes in (b) are Figure 4 restored in (c).

[0116] Step 5: Perform sparse constraint deconvolution on the reconstructed data obtained in Step 4 to obtain the input reflection response. Optimize the wavelet estimation to meet the above assumptions for Marchenko multiple elimination, and at the same time obtain the high-resolution impulse response of the input for Marchenko multiple removal , which improves the noise resistance of the data. The sparse constraint deconvolution model is expressed as:

[0117] (8)

[0118] where, represents the seismic trace, is the wavelet convolution matrix, is the reflection coefficient sequence, represents the noise term. Solve the following optimization problem using the iterative soft threshold algorithm in the Shearlet domain:

[0119] (9)

[0120] where, represents the reflection coefficient in the Shearlet domain obtained by solving, is the wavelet operator, represents the noise level.

[0121] Figure 5 The effect after sparse constraint deconvolution processing is shown. Compared with Figure 5 the processing result of the conventional pulse deconvolution in (a), the result of the improved method of the present invention Figure 5The event axis of the reflection response in (b) is more converged and close to the sharp pulse sequence, with higher resolution, which provides preliminary guarantee for the noise resistance of the data in the subsequent processing flow.

[0122] Marchenko multiple wave removal process:

[0123] Step 6: Input the input seismic records and input reflection responses obtained in steps 4 and 5 into the improved Marchenko multiple wave removal process to prepare for wavelet solution. The final suppressed multiple wave form after wavelet solution is expressed as:

[0124] (10)

[0125] in, is the original reflection response including multiple waves, The sum of the even-numbered terms of the wavelets obtained represents an inverse scattering term with the same energy but opposite polarity to the energy of the multiple waves. and Add and general The multiple waves in the reflection response are suppressed to obtain the reflection response without multiple waves. .

[0126] In the traditional Marchenko method, the two sides of the Green function expression based on the self-focusing theory of formula (11) and formula (12) are combined with the continuation operator (the downward Green's function of the direct wave) is convolved to obtain formulas (13) and (14):

[0127] (11)

[0128] (12)

[0129] (13)

[0130] (14)

[0131] in, represents the integrated time variable; Indicates the spatial position coordinates, its subscript and Represent the depth of the spatial position, represents a spatial position with a depth of 0, Indicates the depth The spatial position of It is used to distinguish the spatial positions in the marking formula (11) and formula (12) separately and convert them into formula (13) and formula (14) after convolution. Spatial position; Indicates the plane at a depth of 0.

[0132] The input reflection response obtained at 5, where Indicates the position of the seismic source, Indicates the receiver position, and both of these parts are on .

[0133] By processing the observed data, an expression for the relationship between the seismic event and the receiver position is obtained. Through the conversion in the time domain, the above Green's function expression is derived therefrom, where and respectively represent the downgoing and upgoing Green's functions when the seismic source is located at and the geophone is located at (it can be seen from the reciprocity theorem that the positions of the seismic source and the geophone can be swapped). and represent the downgoing and upgoing focusing functions.

[0134] The Marchenko multiple removal method convolves both sides of equations (10) and (11) with the continuation operator and then projects the observation system back to the surface again. represents the Green's function and convolved with (the superscripts "+" and "-" represent the downgoing wavefield and the upgoing wavefield respectively); is the convolution form of the coda component. Similarly, [[ID=!47]]is the convolution form of.

[0135] represents a spatially band-limited two-dimensional [[ID=!functions]], while represents the time-domain function. The focusing function and the Green's function have different time distribution characteristics. The present invention defines as the arrival time of the first-arrival wave (i.e., the shortest two-way travel time from the surface to the focusing plane and back). and The time interval of is . Within this interval , can be determined by solving equations (13) and (14). Substituting into equation (13) gives the final expression of [[ID=7!6]] :

[0136] It should be noted that there seems to be an error in the original text where "!functions" and "!47" and "!76" are likely incorrect notations. This translation is based on the best understanding of the provided text with these possible errors considered. (15)

[0137] Define as , collect each corresponding value, and represent the final data without multiples in the form of scattered coda waves as:

[0138] (16)

[0139] where is the original input reflection response, is the reflection response after removing multiples, represents the even coda waves, represents a signal with the same amplitude as the interlayer multiples but opposite polarity. By adding these, the seismic data without interlayer multiples on the left side of the equation can be obtained. In this step is known. To find , it needs to be iteratively solved in the subsequent process .

[0140] Step 7, in the final form of Step 6, solve the coda wave by iterative calculation, that is, expand the improved anti-noise Marchenko multiple removal iterative process of the present invention. Its subscript represents the number of iterative updates. Initialize the first iteration of the coda wave:

[0141] (17)

[0142] where represents time, represents the shortest two-way travel time of the first arrival wave, represents a very small smooth window function quantity, represents the current receiver position, represents the current shot point position.

[0143] Step 8, odd-numbered iteration. The odd-numbered coda wave iteration follows the following formula:

[0144] (18)

[0145] where represents the receiver coordinates when the source is excited at the position, is the integration time variable.

[0146] Step 9, for even-numbered iterations, the coda calculation follows the following formula. Compared with the traditional Marchenko multiple elimination algorithm, the present invention introduces a signal-to-noise separation algorithm based on the Bayesian principle in the iterative solution process of the even-numbered terms to obtain the sum of the coda of the even-numbered terms. ;

[0147] (19)

[0148] In the final form of Step 7, it can be seen from the present invention that the final data without interlayer multiples is only related to the coda of the even-numbered terms. Therefore, in order to ensure the accuracy of the coda of the even-numbered terms during the iteration process, in this embodiment, a signal-to-noise separation algorithm based on the Bayesian principle is introduced in the iterative solution process of the even-numbered terms. After this processing, the calculated coda of the even-numbered terms can get rid of the noise that may be amplified during the iteration process and remains in the preprocessing step due to considering amplitude preservation, and can more accurately suppress the multiple signal during the iteration process. The specific explanation is as follows:

[0149] When solving the coda of the even-numbered terms, first predict the effective terms and noise terms in the input coda in the Shearlet domain, and then obtain the Shearlet coefficients and that maximize the Bayesian conditional probability, and then perform Bayesian subtraction on the noise terms:

[0150] (20)

[0151] The Shearlet coefficients and that maximize the Bayesian conditional probability are obtained by solving the following optimization problem: [[ID=CO33]]

[0152] (21)

[0153] Step 10, substitute the accurate sum of the coda of the even-numbered terms obtained in Step 9 into formula (16) of Step 6 to finally obtain the data without interlayer multiples . In this experiment, the single-shot suppression result of the interlayer multiples can be seen Figure 6 . Compared with the input original data Figure 6 (a), Figure 6 the conventional Marchenko multiple elimination method represented by (b) can only barely eliminate the in-phase axis corresponding to the interlayer multiple indicated by a white arrow, while Figure 6(c) In the result of the improved compressive sensing Marchenko multiple suppression system of the present invention represented, it can be clearly seen that the multiple event axes indicated by several white arrows are cleanly eliminated. Moreover, the event axis indicated by the black arrow is not new energy introduced relative to the original data. Instead, before the multiple suppression, the event axes of the multiple waves and the primary waves coincide and cancel each other out. After the multiple suppression, the primary wave signal is restored, and the event axis indicated by the black arrow is the restored primary wave information. In addition, Figure 7 and Figure 8 respectively show the common offset gather and the autocorrelation spectrum before and after the suppression. It can be seen that the multiple wave signals in the rectangular frame before the suppression are attenuated after the suppression, which proves the processing ability of the system proposed by the present invention in the multi-shot domain, and can successfully suppress the multiple wave interference in the land field measured data, greatly improving the accuracy of subsequent data interpretation.

[0154] Note that the above is only the preferred embodiment of the present invention and the applied technical principle. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described here. Various obvious changes, re-adjustments and substitutions can be made by those skilled in the art without departing from the protection scope of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments. Without departing from the concept of the present invention, more other equivalent embodiments can be included, and the scope of the present invention is determined by the scope of the appended claims.

Claims

1. A Marchenko multiple suppression method for land measured data based on compressive sensing, characterized in that Including the following steps: Step 1, extract three-dimensional seismic data and perform regularization processing to obtain the original data ; Step 2, according to the theory of compressive sensing, perform Shearlet transform on the original data obtained in Step 1 : （1） Among them, represents the Shearlet transform, is the Shearlet coefficient after the transform; First, construct the Shearlet basis function : （2） Among them, is the scale matrix, is the shear matrix, is the scale parameter, is the shear parameter, is the translation parameter, and the scale matrix and the shear matrix are defined in the following form: （3） Then, perform Shearlet transform on the original data obtained in step 1 : （4） Step 3, perform denoising processing in the Shearlet domain, and obtain the denoised data from formula (5): （5） Among them, is the original seismic data to be denoised that has been observed, is the sparse representation matrix, denotes the l1 norm of, denotes the error threshold, is the solution to solve the above denoising problem, denotes the solution estimate, denotes the data after denoising; Step 4, perform reconstruction in the Shearlet domain to obtain the input seismic record for the multiple wave removal step; specifically, interpolate the missing parts of the denoised data to obtain high-density sampled seismic data without missing parts, and perform sparse reconstruction of the seismic data in the Shearlet domain, so as to transform the problem of interpolative reconstruction of missing seismic traces into the following L1 norm minimization problem, and obtain the reconstructed data from Equation (7): （7） Among them, represents the inverse Shearlet operator, is the synthesis interpolation operator, is the Lagrange multiplier, represents the target eigenvector, represents the regularized solution, represents the data after reconstruction; Step 5: Perform sparse-constrained deconvolution on the reconstructed data obtained in Step 4 to obtain the input reflection response. Specifically, optimize the wavelet estimation to meet the assumptions of Marchenko multiple elimination while obtaining the high-resolution impulse response of the input for Marchenko multiple removal. , and the sparse-constrained deconvolution model is expressed as: （8） Among them, represents a seismic trace, is a wavelet convolution matrix, is a reflection coefficient sequence, represents the noise term; Solve formula (9) using the iterative soft threshold algorithm in the Shearlet domain, and obtain the final reflection coefficient in the Shearlet domain through the action of the inverse Shearlet transform operator The corresponding reflection coefficient in the time domain corresponding to the reflection coefficient in the Shearlet domain : （9） Among them, represents the Shearlet domain reflection coefficient obtained by solving, is the wavelet operator, represents the noise level; Step 6, use the input seismic record and input reflection response obtained in Step 4 and Step 5 to perform Marchenko multiple wave removal, and then perform wavelet solving. The final multiple wave suppression form after wavelet solving is expressed as: （10） Among them, is the original reflection response containing multiple waves, is the sum of the even terms of the obtained wavelet, representing an inverse scattering term with the same energy as the multiple wave but opposite in polarity. After and are added together, the multiple waves in are suppressed to obtain the reflection response without multiple waves ; Step 7, solve for the coda wave by iterative calculation , where the subscript represents the number of iterative updates, and initialize the first iteration of the coda wave: （17） wherein, represents time, represents the shortest two-way travel time of the first arrival wave, represents a very small smooth window function quantity, represents the current receiver position, represents the current shot point position; Step 8, odd-numbered iteration. The odd-numbered coda wave iteration follows the following formula: （18） Among them, represents the receiver coordinates when the seismic source is excited at the position, is the integration time variable; Step 9, even-numbered iteration. The aftershock calculation follows the following formula. During the iterative solution of even-numbered terms, a signal-to-noise separation algorithm based on the Bayesian principle is introduced to obtain the sum of aftershocks of even-numbered terms ; （19） Among them, represents the integration time variable; Step 10: Substitute the sum of the even-term tails obtained in Step 9 into Step 6 to finally obtain data without internal multiples .

2. The Marchenko multiple suppression method for land measured data based on compressive sensing according to claim 1, characterized in that Step 1 is specifically as follows: Extract the data stored by time * number of channels as three-dimensional seismic data of time * number of channels * number of shots, and regularize the original seismic data. Regularize the data with relatively large shot intervals and geophone intervals into high sampling density, that is, sparse data with relatively small shot intervals and geophone intervals. The channels with original valid values are regularized to the corresponding positions, and the positions without data are assigned zero values to obtain the original data 。 3. A Marchenko multiple suppression method for land measured data based on compressive sensing according to claim 1, characterized in that, Step 3, solve the seismic data without noise through the hard threshold algorithm, specifically: （6） Among them, represents the independent variable of the function substituted into the hard threshold algorithm, represents the determination threshold, represents the hard threshold function.

4. A Marchenko multiple suppression method for land measured data based on compressive sensing according to claim 1, characterized in that: Step 4, an iterative threshold cooling strategy is adopted to determine the Lagrange multiplier .

5. A Marchenko multiple suppression method for land measured data based on compressive sensing according to claim 1, characterized in that Step 6, the Marchenko multiple removal step is as follows: Convolve both sides of the Green's function expressions based on the autofocus theory in Equation (11) and Equation (12) with the continuation operator , which is the downgoing Green's function of the direct wave, to obtain Equation (13) and Equation (14): （11） （12） （13） （14） Among them, represents the spatial position coordinates, and its subscript and respectively represent the depth of the spatial position, represents the spatial position with a depth of 0, represents the spatial position with a depth of , and the subscript is used to separately distinguish the spatial positions in formula (11) and formula (12) after convolution and converted into the spatial positions in formula (13) and formula (14); represents the plane when the depth is 0; 5 The obtained input reflection response, where, Indicates the position of the seismic source, Indicates the receiver position, and both of these two parts are on ; and respectively represent the downgoing and upgoing Green's functions with the seismic source located at and the geophone located at , and and represent the downgoing and upgoing focusing functions; denotes the Green's function and with the convolution result of, where the superscripts "+" and "-" represent the downgoing wavefield and the upgoing wavefield, respectively; is the convolution form of the coda wave component, is the convolution form of; Representation space is limited to two dimensions function, and Representing the time domain function; definition for First arrival time, i.e., the shortest two-way travel time from the surface to the focal plane and back; and The time interval is ; In this range , It can be determined by solving formula (13) and formula (14); Substitute into formula (13) to obtain 's final expression: （15） Define as , collect the corresponding values, and represent the final data without multiples in the form of scattered coda waves as: （16） Among them, is the original input reflection response, is the reflection response after multiple elimination, represents the even-numbered coda waves, represents a signal with the same amplitude as the interlayer multiple waves but opposite polarity.

6. A Marchenko multiple suppression method for land measured data based on compressive sensing according to claim 1, characterized in that Step 9 is specifically as follows: When solving the even-order coda wave, first predict the valid terms and the noise terms in the input coda wave in the Shearlet domain, and then obtain the Shearlet coefficients that maximize the Bayesian conditional probability and , and then perform Bayesian subtraction on the noise terms: （20） The Shearlet coefficients that maximize the Bayesian conditional probability are obtained by solving formula (21). and : （21）。

Images