Land measured data Marchenko multiple suppression method based on compressed sensing

Abstract

The invention belongs to the technical field of geophysical exploration, and relates to a land measured data Marchenko multiple suppression method based on compressed sensing, comprising a data preprocessing link: encrypting and regularizing original data; carrying out Shearlet de-noising and reconstruction; carrying out sparse constraint deconvolution; a step of suppressing multiple waves in iteration: initializing the first iteration of a wake wave; carrying out odd number and even number iterations to solve a tail wave; signal-noise separation based on the Bayesian principle is carried out in even iteration, and the sum of even term tail waves of suppression noise is obtained. According to the method, higher-quality input data is obtained through an optimized preprocessing process, a more anti-noise Marchenko iteration process is carried out, multiple signals in actually measured land data are effectively suppressed under the condition that any model information and prediction do not need to be subtracted, and a good foundation is laid for subsequent imaging and interpretation; the accuracy of oil and gas resource exploration is improved, and more cost consumed by manual false horizon judgment is saved.

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 and returns to the ground after a single reflection from the underground interface, and it can truly reflect information such as the depth of the formation interface and the structural morphology. Different from this, the multiple wave is an interference wave formed by multiple reflections between underground impedance interfaces, which usually causes false structural information to appear in the seismic imaging result and interferes 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 surface multiples with stronger energy and periodicity, interlayer multiples have non-periodic characteristics and more serious aliasing with primary waves. 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, achieving accurate reservoir identification and interpretation, and accurately locating oil and gas resources.

[0003] Conventional multiple suppression methods, such as the surface-related multiple elimination (SRME) method, 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 subtraction. Such a prediction subtraction strategy is based on some characteristic differences between primary waves and multiples to predict multiples (such as amplitude, waveform, phase, etc.). However, in actual data, especially in mid-deep land data, interlayer multiples usually develop implicitly and the aliasing phenomenon between interlayer multiples and primary waves is very serious. If following such a strategy, subtracting after predicting multiples is very likely to damage the effective signal of primary waves.

[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, based on the process of obtaining the initial extrapolation operator depending on the velocity model in the traditional Marchenko method to extrapolate to the subsurface on one side of the shot point or the receiver point, a convolution of an extrapolation operator is performed on both sides of the expression of the extrapolated Green's function, and the shot point and the receiver point are projected back to the surface again. In this forward process and inverse process, the effects of the extrapolation operator estimated depending on 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 a 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: (1) There is no noise interference in the data; (2) The data needs to meet the ultra-high density sampling condition; (3) The wavelet is approximately a pulse sequence, and the reflection response converges to a sharp pulse signal.

[0006] 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 is often not a pulse sequence.

[0007] Based on this, it is urgent 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

[0008] The purpose of the present invention is to provide an improved Marchenko multiple wave suppression method based on compressive sensing for 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.

[0009] The purpose of the present invention is achieved through the following technical solutions: A Marchenko multiple wave suppression method for measured onshore data based on compressive sensing includes the following steps: Step 1, extract the three-dimensional seismic data and perform regularization processing to obtain the original data ; Step 2: According to the compressive sensing theory, perform Shearlet transform on the original data obtained in Step 1 : (1) where represents the Shearlet transform, is the Shearlet coefficient after transformation; First, construct the Shearlet basis function : (2) 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: (3) Then, perform Shearlet transform on the original data obtained in Step 1 : (4) Step 3: Perform denoising processing in the Shearlet domain. The denoised data is obtained from formula (5): (5) where is the observed original seismic data to be denoised, is the sparse representation matrix, represents the l1 norm of represents the error threshold, is the solution to solve the above denoising problem, represents the solution estimate value, represents the denoised data; Step 4: Perform reconstruction in the Shearlet domain to obtain the input seismic record for the multiple wave removal link; specifically, interpolate the missing part 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 to convert the problem of interpolating and reconstructing the missing seismic traces into the following L1 norm minimization problem. The reconstructed data is obtained from formula (7): (7) where represents the Shearlet inverse operator, is a synthetic interpolation operator, is a Lagrange multiplier, represents the target eigenvector, represents the regularized solution, represents the reconstructed data; 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 satisfy the assumptions of Marchenko multiple elimination while obtaining the high-resolution impulse response of the Marchenko multiple removal input. , and the sparse-constrained deconvolution model is expressed as: (8) where, represents a seismic trace, is the wavelet convolution matrix, is the reflection coefficient sequence, represents the noise term; Solve formula (9) using the iterative soft-threshold algorithm in the Shearlet domain. After the action of the inverse Shearlet transform operator the final Shearlet-domain reflection coefficient obtained The corresponding time-domain reflection coefficient : (9) where, represents the obtained Shearlet-domain reflection coefficient, is the wavelet operator, represents the noise level; Step 6: Use the input seismic record and the input reflection response obtained in Steps 4 and 5 to perform Marchenko multiple removal, and then perform wavelet solution. The final form of suppressing multiples after wavelet solution is expressed as: (10) where, is the original reflection response containing multiples, is the sum of the even terms of the obtained wavelet, representing an inverse scattering term with the same energy and opposite polarity to the multiple energy. After and summation, the multiples in are suppressed to obtain the reflection response without multiples ; Step 7: Solve the coda wave through iterative calculation , where the subscript represents the number of iterative updates. Initialize the first iteration of the coda wave: (17) 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; Step 8, odd-numbered iteration. The odd-numbered coda wave iteration follows the following formula: (18) where, represents the receiver coordinates when the source is excited at the position, is the integration time variable; Step 9, even-numbered iteration. 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 coda waves of the even-numbered terms ; (19) where, represents the integration time variable; Step 10, substitute the sum of the coda waves of the even-numbered terms obtained in Step 9 into Step 6 to finally obtain the data without multiple reflections between layers .

[0010] Furthermore, Step 1 is specifically: extract the originally time * number of traces stored data into three-dimensional seismic data of time * number of traces * number of shots, and regularize the original seismic data. Regularize the data with originally relatively large shot spacing and geophone spacing into high sampling density, that is, sparse data with relatively small shot spacing and geophone spacing. The originally existing valid traces are regularized to the corresponding positions, and the positions without data are assigned zero values to obtain the original data .

[0011] Furthermore, in Step 3, solve the seismic data without noise through the hard threshold algorithm. Specifically: (6) where, represents the independent variable of the function substituted into the hard threshold algorithm, represents the decision threshold, represents the hard threshold function.

[0012] Furthermore, in Step 4, an iterative threshold cooling strategy is adopted to determine the Lagrange multiplier .

[0013] Further, in step 6, the Marchenko multiple removal step is as follows: Multiply both sides of the Green's function expressions based on the autofocus theory in Equation (11) and Equation (12) by the continuation operator , that is, convolve with the downward Green's function of the direct wave, to obtain Equation (13) and Equation (14): (11) (12) (13) (14) where 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 is used to separately distinguish and mark the spatial positions in Equation (11) and Equation (12) that are converted to the spatial positions in Equation (13) and Equation (14) after convolution; represents the plane at a depth of 0; 5 The input reflection response obtained, where indicates the position of the source, indicates the receiver position, and both of these parts are on ; and respectively represent the downward and upward Green's functions when the source is located at and the geophone is located at , and represent the downward and upward focusing functions; represents the Green's function and and convolution result, where the superscripts "+" and "-" respectively represent the downward wavefield and the upward wavefield; is the convolution form of the coda wave component, is the convolution form; represents the spatially band-limited two-dimensional function, while represents the time-domain function; Define 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); Substitute into equation (13) to obtain The final expression of: (15) Define as , collect each corresponding value, and represent the final data without multiples in the form of scattered coda waves as: (16) where is the original input reflection response, is the reflection response after eliminating multiples, represents the even-numbered coda wave, represents a signal with the opposite polarity and equivalent amplitude to the interlayer multiples.

[0014] Furthermore, step 9 is specifically: when solving the even-numbered coda wave, first predict the effective term and the noise term in the input coda wave 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 term: (20) The Shearlet coefficients and that maximize the Bayesian conditional probability are obtained by solving equation (21): (21).

[0015] Compared with the prior art, the beneficial effects of the present invention are: 1. The present invention provides an improved Marchenko interlayer multiple suppression method for 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, ultimately forming a complete multiple suppression system applicable to onshore seismic data. The improved multiple suppression system includes a data preprocessing 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.

[0016] 2. The present invention improves the commonly used industrial simple spline interpolation method 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 multiples more accurate. For the noise that may be amplified during the iteration due to amplitude preservation considerations, the present invention improves the traditional Marchenko iteration to an iterative process with high noise resistance based on the Bayesian principle of Shearlet noise separation, successfully suppressing the residual noise and possible new noise effectively during each iteration process.

[0017] 3. The present invention is applicable to the suppression of interlayer multiples in onshore measured data, successfully breaking through the application limitations of the Marchenko multiple removal method that are difficult to overcome in actual onshore data. By optimizing the preprocessing process to obtain higher-quality input data that fits the assumptions and performing a more noise-resistant Marchenko iterative suppression process, the interlayer multiple signals in the measured onshore data can be effectively and accurately suppressed ultimately, and no model information or predictive subtraction is 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

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

[0019] Figure 1 Flowchart of the steps of the Marchenko multiple suppression method for onshore measured data based on compressive sensing; Figure 2 is the original measured onshore data; Figure 3 For the display of the denoising effect of Shearlet; among them, (a) is the original data input in 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 with the conventional method; Figure 4 For the display of the reconstruction effect of Shearlet after denoising; among them, (a) is the original data before denoising and reconstruction and the local magnification, (b) is the data and local magnification after conventional filtering denoising and conventional spline interpolation, and (c) is the data and local magnification after improved Shearlet denoising and reconstruction; Figure 5 For the display of the effect of sparse constraint deconvolution; among them, (a) is the reflection response obtained by the conventional impulse deconvolution method and the local magnification, and (b) is the reflection response obtained by the sparse constraint deconvolution method and the local magnification; Figure 6 For the display of the effect after processing by the improved compressive sensing Marchenko interlayer multiple suppression system of the present invention; (a) is the original data before suppression and the local magnification, (b) is the result of conventional Marchenko multiple suppression, and (c) is the effect after the improved compressive sensing Marchenko interlayer multiple suppression of the present invention; Figure 7 For the 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; Figure 8 For the 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

[0020] The present invention will be further described below in conjunction with embodiments: The present invention will be further described in detail below in conjunction with the accompanying 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 convenience of description, only the parts related to the present invention are shown in the drawings, rather than all the structures.

[0021] 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.

[0022] 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 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 improved data preprocessing and improved Marchenko multiple suppression steps. 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.

[0023] Specifically, in the preprocessing step: 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 not close to a sharp pulse signal, the present invention replaces the commonly used pulse deconvolution method in industry with a sparse-constrained deconvolution method. 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 condition of Marchenko.

[0024] In the link of Marchenko multiple elimination: The present invention has improved the iterative algorithm for conventional Marchenko multiple elimination. The conventional Marchenko multiple removal process usually has some residual noises after preprocessing considering the amplitude preservation of the primary wave. 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 has improved 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 set of Marchenko interlayer multiple suppression methods 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.

[0025] Example 1: A method for suppressing Marchenko multiples in onshore measured data based on compressive sensing has been applied to a piece of onshore measured data for effect analysis. The original data has a sampling interval of 8 ms, a total sampling time of 6 s, an original shot interval of 25 m, and an original trace interval of 25 m. In the preprocessing link, both the shot interval and the trace interval are encrypted to 12.5 m.

[0026] As Figure 1 shown, it includes the following steps: Data preprocessing link: Step 1, the establishment of the acquisition system. The data originally stored in time * number of traces is extracted into three-dimensional seismic data of time * number of traces * number of shots. At the same time, the original seismic data is regularized. The data with relatively large shot intervals and geophone intervals is regularized into sparse data with a 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 is encrypted into high-sampling-density data in the subsequent interpolation step. The original data after this step of processing is shown in Figure 2 . It can be seen that there are relatively strong interfering surface wave noises in the original data after preliminary processing. Since the data has not been reconstructed and encrypted yet, it shows an undersampled sparse form. The signal is not close to a sharp pulse and there are some high-energy noise interferences before 0.4 s, and these interferences are subsequently removed.

[0027] Step 2: According to the theory of compressive sensing, perform Shearlet transform on the original data obtained in Step 1 as follows: (1) wherein, denotes the Shearlet transform, is the Shearlet coefficient after transformation, which is used for subsequent denoising and other processing.

[0028] Specifically: First, construct the Shearlet basis function as follows: (2) wherein, is the scaling matrix, is the shear matrix, is the scaling parameter, is the shear parameter, is the translation parameter. The scaling matrix and the shear matrix are defined in the following form: (3) Then, perform Shearlet transform on the original data obtained in Step 1 as follows: (4) Step 3: Perform denoising processing in the Shearlet domain. Since the noise distribution in the Shearlet domain is different, transform the denoising problem into solving the noise-free data: (5) wherein, is the original seismic data to be denoised observed, is the sparse representation matrix, denotes the l1 norm of, denotes the error threshold, is the solution to 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: (6) wherein, denotes the independent variable of the function substituting the hard thresholding algorithm, denotes the decision threshold, denotes the hard thresholding function.

[0029] In this experiment, the specific results of Shearlet denoising can be 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), and has removed more steeply dipping noise represented by the white arrows compared to the results after conventional filtering denoising processing Figure 3 (b), Figure 3 (d) represents the noise signal that Shearlet denoising removes more than the conventional filtering denoising method.

[0030] Step 4: Reconstruct in the Shearlet domain to obtain the input seismic record for the multiple wave 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 wave 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: (7) where represents the Shearlet inverse operator, is the synthetic interpolation operator, is the Lagrange multiplier, represents the target eigenvector, represents the regularized solution, represents the reconstructed data. This denoising problem aims to minimize the Shearlet coefficients. To this end, the present invention adopts an iterative threshold cooling strategy to determine the Lagrange multiplier .

[0031] In the present invention, Shearlet reconstruction is performed on the basis of the previous step of Shearlet denoising. Figure 4 shows the effect comparison after denoising and reconstruction. 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. By comparing Figure 4 (a), Figure 4 (b) and Figure 4 (c), it can be seen that Figure 4 the steeply dipping noise within the black dashed parallelogram frame in Figure 4 (b) is eliminated in Figure 4 (c). Compared with Figure 4 (b),Figure 4 The in-phase axes at many positions indicated by the white arrows in (c) are significantly more naturally transitional, and Figure 4 many discontinuous weak in-phase axes in (b) Figure 4 are recovered in (c).

[0032] Step 5: Perform sparse-constrained deconvolution on the reconstructed data obtained in Step 4 to obtain the input reflection response. Optimize the wavelet estimation to meet the above assumptions of Marchenko multiple elimination, and simultaneously obtain the high-resolution impulse response of the input for Marchenko multiple removal, which improves the noise resistance of the data. The sparse-constrained deconvolution model is expressed as: (8) 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: (9) where, represents the reflection coefficient in the Shearlet domain obtained by solving, is the wavelet operator, represents the noise level.

[0033] Figure 5 shows the effect after sparse-constrained deconvolution processing. Compared with Figure 5 the processing result of conventional pulse deconvolution in (a), the result of the improved method of the present invention Figure 5 has in-phase axes of the reflection response in (b) that are more convergent and closer to a sharp pulse sequence, with higher resolution, which provides a preliminary guarantee for the noise resistance of the data in the subsequent processing flow.

[0034] Marchenko multiple removal section: Step 6: Input the input seismic record and the input reflection response obtained in Step 4 and Step 5 into the improved Marchenko multiple removal section to prepare for wavelet solution. The final form of suppressing multiples after wavelet solution is expressed as: (10) where, is the original reflection response containing multiples, is the sum of the even terms of the calculated wavelet, representing an inverse scattering term with the same energy as the multiple energy but opposite in polarity. After and are added together, Multiple reflection suppression in [it] results in a reflection response without multiple reflections .

[0035] In the Marchenko multiple reflection removal method, in the traditional Marchenko method, both sides of the Green's function expressions based on the autofocus theory in equations (11) and (12) are convolved with the continuation operator (the downward Green's function of the direct wave) to obtain equations (13) and (14): (11) (12) (13) (14) where represents the integrated time variable; 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 is used to separately distinguish and mark the spatial positions in equations (11) and (12) that are converted to the spatial positions in equations (13) and (14) after convolution; represents the plane at a depth of 0.

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

[0037] 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. Here and respectively represent the downward and upward 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 swapping the positions of the seismic source and the geophone is the same), and represent the downward and upward focusing functions.

[0038] The Marchenko multiple reflection removal method convolves both sides of equations (10) and (11) with the continuation operator After convolution, the acquisition system is projected back to the surface again. denotes the Green's function and with (the superscripts "+" and "-" represent the downgoing and upgoing wavefields respectively); is the convolution form of the coda component. Similarly, is the convolution form of

[0039] denotes the spatially band-limited two-dimensional function, while represents the function in the time domain. 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 are in the time interval . Within this interval , can be determined by solving equations (13) and (14). Substituting into equation (13) gives the final expression of : (15) Define as , collect the corresponding values for each , and represent the final data without multiples in the form of scattered coda waves as: (16) where is the original input reflection response, is the reflection response after removing multiples, represents the even-numbered coda waves, represents the signal with the opposite polarity and equivalent amplitude to the interlayer multiples. By this summation, the seismic data without interlayer multiples on the left side of the equation can be obtained. In this step is known. To find , needs to be iteratively solved in subsequent processes.

[0040] Step 7, in the final form of step 6, iteratively calculate to solve for the coda wave , that is, expand the improved anti-noise Marchenko multiple removal iterative process of the present invention. The subscript represents the number of iterative updates. 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.

[0041] Step 8, odd-numbered iteration, the odd-numbered coda wave iteration follows the following formula: (18) Wherein, represents the receiver coordinates when the seismic source is excited at the position, is the integration time variable.

[0042] Step 9, even-numbered iteration, the coda wave calculation follows the following formula. Compared with the traditional Marchenko multiple wave 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 waves of the even-numbered terms ; (19) In the final form of Step 7, it can be seen in the present invention that the final data without interlayer multiple waves is only related to the coda waves of the even-numbered terms. Therefore, in order to ensure the accuracy of the coda waves 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 waves of the even-numbered terms can get rid of the noise that may be amplified during the iteration process due to the consideration of amplitude preservation in the preprocessing step, and can more accurately suppress the multiple wave signals during the iteration process. The specific explanation is as follows: When solving the coda waves of the even-numbered terms, first predict the effective terms and the 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: (20) The Shearlet coefficients and that maximize the Bayesian conditional probability are obtained by solving the following optimization problem: (21) Step 10, the accurate sum of the coda waves of the even-numbered terms obtained in Step 9 Substitute it into the formula (16) in step 6, and finally obtain the data without multiple reflections between layers. , in this experiment, the single-shot suppression result of multiple reflections between layers can be seen Figure 6 . Compared with the original input data Figure 6 (a), Figure 6 the conventional Marchenko multiple reflection elimination method represented by (b) can only barely eliminate the event corresponding to the multiple reflection between layers indicated by one white arrow, while Figure 6 in the result of the improved compressive sensing Marchenko multiple reflection suppression system of the present invention represented by (c), it can be clearly seen that the multiple reflection events indicated by several white arrows are cleanly eliminated, and the event indicated by the black arrow is not new energy introduced relative to the original data, but is due to the coincidence and mutual cancellation of the events of the multiple reflection and the primary wave before the multiple reflection is suppressed. After the suppression of the multiple reflection, the primary wave signal is restored, and the event 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 autocorrelation spectrum before and after suppression. It can be seen that the multiple reflection signals in the rectangular frame before suppression are attenuated after suppression, which proves the processing ability of the system proposed in the present invention in the multi-shot domain, and can successfully suppress the multiple reflection interference in the land field data, greatly improving the accuracy of subsequent data interpretation.

[0043] Note that the above is only the preferred embodiment of the present invention and the technical principles applied. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described here, and 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 compressed sensing, characterized in that: The following steps are involved: Step 1: Extract 3D seismic data and perform regularization processing to obtain raw data ; Step 2: According to the theory of compressed sensing, the original data obtained in step 1 is Perform a Shearlet transformation: （1） in, represents the Shearlet transform, is the transformed Shearlet coefficient; First construct the Shearlet basis function : （2） in, is the scale matrix, is the shear matrix, is the scale parameter, is the shearing parameter, is the translation parameter, the scale matrix and the shear matrix are defined as follows: （3） Then the original data obtained in step 1 Perform a Shearlet transformation: （4） Step 3: Perform denoising in the Shearlet domain and obtain the denoised data using formula (5): （5） in, is the observed original seismic data to be denoised, is a sparse expression matrix, express The one-norm of represents the error threshold, To solve the above denoising problem, represents the estimated value of the solution, represents the denoised data; Step 4: Reconstruct in the Shearlet domain to obtain the input seismic records for the multiple wave removal step; specifically, the denoised data The missing parts are interpolated to obtain high-density sampled non-missing seismic data, and the seismic data is sparsely reconstructed in the Shearlet domain to transform the interpolation reconstruction problem of missing seismic traces into the following L1 norm minimization problem. The reconstructed data is obtained by formula (7): （7） in, represents the Shearlet inverse operator, is the synthetic interpolation operator, is the Lagrange multiplier, represents the target feature vector, Indicates regular solution, Represents the reconstructed data; Step 5: Perform sparse constrained deconvolution on the reconstructed data obtained in step 4 to obtain the input reflection response. Specifically, the wavelet estimation is optimized to meet the assumption of Marchenko multiple elimination and obtain the high-resolution impulse response of the input with Marchenko multiple elimination. , the sparse constrained deconvolution model is expressed as: （8） in, represents the earthquake channel, is the wavelet convolution matrix, is the reflection coefficient sequence, represents the noise term; In the Shearlet domain, an iterative soft threshold algorithm is used to solve formula (9), and after the Shearlet inverse transform operator The final Shearlet domain reflection coefficient obtained by action The corresponding time domain reflection coefficient : （9） in, represents the Shearlet domain reflection coefficient obtained by solution, is the wavelet operator, Indicates the noise level; Step 6, using the input seismic records and input reflection responses obtained in steps 4 and 5 to remove Marchenko multiple waves, and then performing wavelet solution. The final suppressed multiple wave form after solving the wavelet is expressed as: （10） in, is the original reflection response including multiple waves, The sum of the even-numbered terms of the wavelets obtained represents a backscattering term with the same energy but opposite polarity to the energy of the multiple waves. and plus general The multiple waves in the reflection response without multiple waves are suppressed. ; Step 7: Solve the coda by iterative calculation , whose subscript indicates the number of iterative updates, initializes the first iteration of the tail wave: （17） in, Indicates time, represents the shortest two-way travel time of the first arrival wave, represents a very small smooth window function, represents the current receiver position, Indicates the current gun point position; Step 8, odd number of iterations, odd number of tail wave iterations follow the following formula: （18） in, Indicates that the earthquake source is The coordinates of the receiving point when the position is triggered, is the integrated time variable; Step 9, even-numbered iterations, the coda calculation follows the following formula. In the iterative solution process of even-numbered terms, the signal-noise separation algorithm based on the Bayesian principle is introduced to obtain the sum of the codas of the even-numbered terms: ; （19） in, represents the integrated time variable; Step 10: sum the tail waves of the even-order terms obtained in step 9 Substitute it into step 6, and finally obtain data without interlayer multiple waves .

2. The Marchenko multiple suppression method for land measured data based on compressed sensing according to claim 1, characterized in that: Step 1, specifically: extract the data originally stored in the format of time*channel number into three-dimensional seismic data of time*channel number*shot number, and regularize the original seismic data, regularize the data originally with relatively large shot spacing and detector spacing into sparse data with high sampling density, i.e. relatively small shot spacing and detector spacing, the channels originally with valid values ​​are regularized to the corresponding positions, and the positions without data are assigned zero values ​​to obtain the original data .

3. The Marchenko multiple suppression method for land measured data based on compressed sensing according to claim 1, characterized in that: Step 3: Solve the noise-free seismic data using a hard threshold algorithm, specifically: （6） in, represents the function argument substituted into the hard threshold algorithm, represents the judgment threshold, represents the hard threshold function.

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

5. The Marchenko multiple suppression method for land measured data based on compressed sensing according to claim 1, characterized in that: Step 6, the Marchenko multiple wave removal step is: Substitute both sides of the Green function expression based on the self-focusing theory of formula (11) and formula (12) with the extension operator , that is, the convolution of the downlink Green's function of the direct wave, we get formula (13) and formula (14): （11） （12） （13） （14） in, Indicates spatial position coordinates, with the subscript and They represent the depth of the spatial position, represents the 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 transform them into formula (13) and formula (14) after convolution. Spatial location; Represents a plane with a depth of 0; 5 to obtain the input reflection response, where Indicates the location of the earthquake source, Indicates the receiver location, both parts are in superior; and The earthquake sources are located at , the detection point is located at The downward and upward Green's functions, and represents the downlink and uplink focusing functions; 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 Representation space is limited to two dimensions function, and Representing the time domain Function; Definition for The first arrival time is the shortest two-way travel time from the surface to the focal plane and back; and The time interval is ; In this range , By solving formula (13) and formula (14), it can be determined that; Will Substituting into formula (13) we can get The final expression of is: （15） Will Defined as , collect each Corresponding value, and the final data without multiple waves is expressed in the form of scattered coda waves as: （16） in, is the original input reflection response, To eliminate the reflection response after multiple waves, represents an even-numbered tail wave, Indicates a signal with a polarity quite opposite to that of the interlayer multiple wave amplitude.

6. The Marchenko multiple suppression method for land measured data based on compressed sensing according to claim 1, characterized in that: Step 9, specifically: when solving the coda wave of even-order terms, first predict the effective terms in the coda wave input in the Shearlet domain and the noise term , and then find the Shearlet coefficient that maximizes the Bayesian conditional probability and , and then perform Bayesian subtraction on the noise term: （20） The Shearlet coefficient that maximizes the Bayesian conditional probability is obtained by solving formula (21): and : （21）。

Images