Three-dimensional surface multiple and wavelet estimation method based on accelerated linear bregman and gpu

Abstract

The application belongs to the technical field of oil and gas seismic exploration data processing, and specifically discloses a three-dimensional surface multiple wave and wavelet estimation method based on accelerated linear Bregman and GPU, which comprises the following steps: obtaining original three-dimensional seismic input data, recombining and transforming to obtain a frequency domain three-dimensional data matrix and time domain vectorized data, and transmitting the data to GPU display memory; taking the data as input on the GPU, solving a sparse inversion problem based on a feedback model expression using an accelerated linear Bregman algorithm under the condition of fixed wavelet to obtain a sparse estimation of the Green function; under the condition of fixed sparse estimation of the Green function, directly estimating the wavelet according to the sparse estimation of the Green function using a wavelet analytical solution formula; repeating the above two steps to alternately update until the iteration number or error tolerance is reached, and then separating the surface-related multiple waves and primary waves from the original data. The application can effectively improve the precision and efficiency of three-dimensional seismic data surface multiple wave elimination.

Description

Technical Field

[0001] This application belongs to the field of oil and gas seismic exploration data processing technology, and more specifically, relates to a three-dimensional surface multiple wavelet and wavelet estimation method based on accelerated linear Bregman and GPU. Background Technology

[0002] In actual seismic exploration, seismic waves propagate through complex three-dimensional subsurface media. With the development of modern industrial-grade three-dimensional acquisition technology, both longitudinal (inline) and transverse (crossline) three-dimensional data are acquired simultaneously. Traditional surface multiple elimination methods (such as Surface-Related Multiple Elimination (SRME) and Estimation of Primaries by Sparse Inversion (EPSI)) were initially developed based on two-dimensional wavefield propagation assumptions. However, in geological environments with significant three-dimensional effects, such as drastically undulating surfaces, cable feather effects, complex seafloor topography, or steeply dipping structures, the actual propagation paths (e.g., out-of-plane reflections) and travel times of three-dimensional multiples deviate significantly from the two-dimensional assumptions.

[0003] Suppression and elimination of free-surface correlated multiples has always been an important research direction in the field of seismic data processing. SRME, as the most representative feedback model "prediction-subtraction" method, often fails to achieve ideal results when primary and multiples strongly overlap. To overcome this shortcoming, EPSI fits the observed data by jointly estimating the primary and free-surface correlated multiples, thus avoiding the adaptive subtraction process, and has been proven to be superior to the traditional SRME method. The subsequently proposed Robust Sparse Inversion Primary Estimation Method (REPSI) utilizes the Spectral Projected Gradient for L1 minimization (SPGL1) algorithm to enhance the robustness of the inversion process. As a further development of REPSI, the Estimation of Multiples and Sources (EMS) method employs a simpler and more efficient Accelerated Linearized Bregman (ALB) algorithm to recover the Green's function without a free surface, and completes the wavelet recovery through one-step scaling estimation. Studies have shown that EMS not only shows promise in both synthetic and real-world data, but may also outperform REPSI in terms of convergence speed, surface-correlated multiple suppression, algorithm simplicity, and computational efficiency. Under complex geological conditions, three-dimensional effects are not negligible, and two-dimensional methods are prone to amplitude and travel-time errors. Therefore, the prediction and elimination of free surface-correlated multiples must be performed within a complete three-dimensional framework. However, compared to the more mature two-dimensional methods, the research and application of three-dimensional methods remain relatively insufficient.

[0004] The shortcomings of existing technologies are mainly manifested in the following aspects: First, if two-dimensional algorithms are still used to process three-dimensional data line by line, it will lead to serious errors in amplitude and travel time prediction. The processing effect of surface multiples is "suppression" rather than complete "elimination", and may even damage the effective primary wave. Second, in terms of primary wave estimation methods, conventional SRME is based on "prediction-adaptive subtraction". When primary waves and multiples overlap severely, it will lead to primary wave distortion. EPSI directly inverts the primary wave by solving the feedback model. However, Robust EPSI technology is highly dependent on the complex SPGL1 solver, which is extremely difficult to debug in industrial-grade massive three-dimensional data. Third, three-dimensional seismic data processing involves extremely large data volumes and dense spatial frequency domain cross matrix convolution, which leads to low computational efficiency of conventional central processing unit (CPU) computing architectures and makes it difficult to implement in industrial applications.

[0005] Therefore, improving the accuracy and efficiency of surface multiples elimination in 3D seismic data is a technical problem that urgently needs to be solved. Summary of the Invention

[0006] To address the shortcomings of existing technologies, the purpose of this application is to provide a three-dimensional surface multiples and wavelet estimation method based on accelerated linear Bregman and GPU, which can effectively improve the accuracy and efficiency of surface multiples elimination in three-dimensional seismic data.

[0007] To achieve the above objectives, in a first aspect, this application provides a method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and a graphics processing unit (GPU), comprising the following steps: S10: Acquire the original three-dimensional seismic input data containing surface-correlated multiples; reorganize the original three-dimensional seismic input data in the dimensions of time, common receiver points, and common shot points to obtain a time-domain seismic wavefield; perform a one-dimensional fast Fourier transform (FFT) on the time-domain seismic wavefield along the time axis to obtain a three-dimensional data matrix in the frequency domain; and vectorize the original three-dimensional seismic input data to obtain the corresponding time-domain vectorized data; transfer the three-dimensional data matrix and the time-domain vectorized data to the GPU memory. S20, on the GPU, using the three-dimensional data matrix and the time-domain vectorized data as input, and under the condition of fixed wavelets, based on the feedback model expression, the sparse inversion problem is solved using the accelerated linear Bregman algorithm to obtain a sparse estimate of the Green's function without surface multiple interference; the feedback model expression describes the relationship between the total input seismic data containing all orders of surface-correlated multiples, the Green's function to be optimized, and the wavelets. S30, under the condition that the sparse estimate of the Green's function is fixed, the corresponding wavelet is directly estimated using the analytical solution formula of the wavelet based on the sparse estimate of the Green's function. S40, repeat S20 and S30, alternately update the sparse estimate of the Green's function and the wavelet until the set number of iterations or error tolerance is reached and the iteration stops. Based on the sparse estimate of the Green's function and the wavelet obtained when the iteration stops, separate the surface-related multiples and primarys from the original three-dimensional seismic input data.

[0008] As a further preferred embodiment, in step S20, the feedback model expression is a bilinear equation. ,in This represents the three-dimensional data matrix. Denotes the Green's function, Indicates the sub-wavelength, Represents the free surface reflection operator. Indicates a primary wave signal. Represents the full-order surface-correlated multiple wave; the bilinear equation in and When one variable is fixed, the free surface reflection operator degenerates into a linear form with respect to the other variable; Under the simplified assumption of the air-water interface reflectance coefficient, the negative value of the identity matrix is ​​taken.

[0009] As a further preferred embodiment, step S20, which involves using the accelerated linear Bregman algorithm to solve the sparse inversion problem, specifically includes: The feedback model expression is reconstructed into a system of linear equations in the time domain. ,in For the time-domain vectorized data, Let the vector of the Green's function to be determined in the time domain be . The forward operator includes time-frequency transformation, zero-filling, and matrix multiplication operations. Based on the sparse distribution of the Green's function in the time domain, the Green's function is inverted by solving a minimization problem based on L1 norm constraints. The objective function is: ,in Denotes the L1 norm of g. express The L2 norm error between d and d For sparse threshold parameters, Indicates the error tolerance.

[0010] As a further preferred embodiment, the specific steps of the accelerated linear Bregman algorithm are as follows: Obtaining forward operators Observation data vector Initial value vector threshold Maximum number of iterations ; Initialization: Iteration variables , Will Assign to ; Repeat the following steps Second-rate: calculate ; calculate ; Will Updated to +1; Derive sparse solutions = ;in, same , which is the initial solution; lIndex for iteration count; , These are temporary intermediate variables during the iteration process; This is the step size parameter.

[0011] As a further preferred embodiment, in step S30, the wavelet analytical solution formula is: superscript H Represents the complex conjugate transpose of a matrix. This represents the spectrum of a certain frequency component of the sparse estimate of the Green's function obtained through step S20. This represents the spectrum of a certain frequency component of the current estimated wavelet. This represents the spectrum of a certain frequency component of the input data. This represents the spectrum of a surface multiple signal at a certain frequency, predicted using a sparse estimate of the current best Green's function and the input data; using this formula, at each frequency... The wavelet in the current iteration loop is then solved.

[0012] As a further preferred embodiment, in step S10, the time axis of the original three-dimensional seismic input data Coordinates of the main survey line of the receiver point Coordinates of the detection point connection line Main survey line coordinates of the shot point and the coordinates of the shot point connection line The original 3D seismic input data is reshaped using a seismic line index, and its dimensionality is reduced and reorganized into a 3D input data matrix. ,in This represents the number of time sampling points; Total number of detector points ( ); Total gun points ( Under ideal sampling conditions, the shot points and receiver points are arranged to completely overlap, meaning the total number of shot points equals the total number of receiver points. Reorganized time-domain seismic wave field After performing a one-dimensional fast Fourier transform along the time axis, it is transformed into a three-dimensional data volume in the frequency domain. ,in This represents the total number of discrete frequency components.

[0013] As a further preferred embodiment, in step S10, after acquiring the original three-dimensional seismic input data, if the original three-dimensional seismic input data contains direct waves and ghost waves, then before performing reconstruction and fast Fourier transform, the original three-dimensional seismic input data is first preprocessed to remove direct waves and ghost waves.

[0014] As a further preferred embodiment, in step S20, within the GPU hardware, all matrix multiplications and Fourier forward and inverse transformations are mapped to multi-threaded high-concurrency processing, and intermediate variables are cleaned up in real time during the calculation process to release video memory.

[0015] In a second aspect, this application provides an electronic device including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method as described in any one of the above descriptions.

[0016] Thirdly, this application provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the steps of the method as described in any one of the above statements.

[0017] This application achieves the following effects: By directly inverting the Green's function without surface multiple interference based on a feedback model within a 3D framework, the amplitude and travel time prediction biases caused by 3D out-of-plane effects in 2D methods when processing 3D data can be overcome, avoiding damage to the effective primary wave and thus achieving complete elimination of surface-related multiples rather than limited suppression; By using the accelerated linear Bregman algorithm to replace the highly complex SPGL1 solver for sparse constraint inversion, the difficulty of solving large-scale 3D inversion problems and the debugging cost can be reduced. At the same time, by directly estimating the wavelet using the wavelet analytical solution formula, the high computational burden brought by the traditional least squares iterative algorithm can be avoided. The synergistic effect of the two helps to improve the convergence speed and algorithm simplicity of the overall inversion process; Mapping the matrix multiplication and Fourier forward and inverse transformations after data recombination to GPU multi-threaded concurrent execution can cope with the extremely large computational requirements of 3D seismic data volumes and effectively alleviate the computational efficiency bottleneck of conventional CPU architectures; The combined effect of the above methods can significantly reduce the processing time of industrial-grade massive data while ensuring the accuracy of multiple prediction, thereby effectively improving the accuracy and computational efficiency of surface multiple elimination in 3D seismic data. Attached Figure Description

[0018] Figure 1 This is a flowchart of a three-dimensional seismic surface multiple wavelet and wavelet estimation method based on the accelerated linear Bregman algorithm and GPU acceleration provided in the embodiments of this application; Figure 2 This is a three-dimensional visualization of the three-layer horizontal velocity model used for basic verification provided in the embodiments of this application; Figure 3 This is a 3D visualization of the Overthrust 3D velocity model used for constructing adaptive verification, provided in an embodiment of this application. Figure 4This is a three-dimensional display of the simulated total data shot set (including the primary wave and all multiple waves) obtained by excitation at the center position of the three-layer horizontal velocity model provided in the embodiments of this application; Figure 5 The provided three-dimensional coherent source gathers for the three-layer horizontal model used in the algorithm processing of this application are: (a) is the extraction result of the first three-dimensional coherent shot gather, (b) is the extraction result of the middle three-dimensional coherent shot gather, and (c) is the extraction result of the last three-dimensional coherent shot gather. Figure 6 This is a wavefield feature map of a time slice (horizontal slice) of the total data corresponding to the three-layer horizontal model provided in this application embodiment at time 0.5 seconds; Figure 7 This is a comparison of the effects of the central resonance source gather of the three-layer horizontal model provided in this application embodiment after processing by the algorithm of this application; wherein, (a) is the original input shot gather containing the primary wave and all multiple waves, (b) is the three-dimensional surface-related multiple waves accurately predicted and stripped by the method of this application, and (c) is the clean shot gather record output after eliminating surface multiple waves. Figure 8 The wavelet information diagram is estimated by the algorithm of this application after processing the three-layer horizontal model data provided in the embodiments of this application; wherein, (a) is the time domain waveform and (b) is the amplitude spectrum; Figure 9 These are the results of testing the three-layer horizontal model data provided in this application, showing the superiority of the method in this application over the three-dimensional SRME method; where (a) is the original common shot record, (b) is the free surface-related multiples estimated and eliminated by the 3D SRME method, (c) is the primary wave result obtained by 3D SRME processing, (d) is the difference between Figure (c) and Figure (f), (e) is the free surface-related multiples estimated and eliminated by the 3D EMS method, and (f) is the primary wave result obtained by 3D EMS processing; Figure 10 This is a three-dimensional display of the simulated total data shot set (including the primary wave and all multiple waves) excited at the center of the Overthrust model, provided in the embodiments of this application. Figure 11 The provided three-dimensional coherent source gathers for the Overthrust model used in the algorithm processing of this application are: (a) is the extraction result of the first three-dimensional coherent shot gather, (b) is the extraction result of the middle three-dimensional coherent shot gather, and (c) is the extraction result of the last three-dimensional coherent shot gather. Figure 12 This is a time-slice wavefield feature map of the total data of the Overthrust model provided in this application embodiment at time 0.5 seconds; Figure 13The above is a comparison of the effects of the central resonance source gather of the Overthrust model provided in this application processed by the algorithm of this application; where (a) is the original input shot gather containing the primary wave and all multiple waves, (b) is the three-dimensional surface-related multiple waves accurately predicted and stripped by the method of this application, and (c) is the clean shot gather record output after eliminating surface multiple waves. Figure 14 The wavelet information diagram is estimated by processing the Overthrust model data provided in this application embodiment through the algorithm of this application; wherein, (a) is the time domain waveform and (b) is the amplitude spectrum; Figure 15 These are the center survey shot gather processing results provided in this application embodiment on a complex Overthrust 3D model, respectively using the two-dimensional surface multiples elimination method and the three-dimensional surface-related multiples elimination method of this application; wherein, (a) to (c) are the processing results of the three-dimensional surface-related multiples elimination method of this application, (d) are the same original records, (e) are the multiples stripped by the two-dimensional surface multiples elimination method, and (f) are the primary waves after processing by the two-dimensional surface multiples elimination method; Figure 16 These are the results of testing the Overthrust model data provided in this application, showing the superiority of the proposed method over the 3D SRME method; where (a) is the original common shot record, (b) is the free surface-related multiples estimated and eliminated using the 3D SRME method, (c) is the primary wave result obtained by 3D SRME processing, (d) is the difference between Figure (c) and Figure (f), (e) is the free surface-related multiples eliminated by the 3D EMS method, and (f) is the primary wave result obtained by 3D EMS processing; Figure 17 These are diagrams showing the superior performance of the method described in this application compared to the 3D SRME method in the zero-shot-receiver offset data volume profile of the three-dimensional Overthrust model provided in this embodiment. Among them, (a) is the original noisy data volume profile, (b) is the free surface correlated multiples estimated and eliminated by the 3D SRME method, (c) is the primary wave zero offset data after 3D SRME processing, (d) is the difference between Figure (c) and Figure (f), (e) is the three-dimensional surface multiple noise data volume estimated by the 3D EMS method of this application, and (f) is the pure primary wave data volume after final processing at zero shot-receiver offset. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0020] To overcome the shortcomings of existing surface multiples elimination techniques in processing 3D pre-stack seismic data, such as large travel time and amplitude prediction deviations, loss of effective signal, and high computational and parameter tuning costs, and to improve the efficiency and accuracy of surface multiples elimination in 3D seismic data, this application provides a 3D seismic surface multiples and wavelet estimation method based on an accelerated linear Bregman algorithm and GPU acceleration, including the following steps: S1: Pre-stack data volumes acquired through actual 3D seismic exploration typically contain 5 dimensions, namely: time axis Coordinates of the main survey line of the receiver point Coordinates of the detection point connection line Main survey line coordinates of the shot point and the coordinates of the shot point connection line This application first reshapes the original 5D pre-stack seismic data using a seismic line index, reducing its dimensionality and reorganizing it into a three-dimensional input data matrix. .in, This represents the number of time sampling points; Total number of detector points ( ); Total gun points ( Under ideal sampling conditions, the shot points and receiver points are arranged to completely overlap, meaning the total number of shot points equals the total number of receiver points. Reorganized time-domain seismic wavefield Perform a one-dimensional FFT along the time axis to transform it into a three-dimensional data volume in the frequency domain. ,in This represents the total number of discrete frequency components. It should be clearly stated that the term "total data" used subsequently in this application refers only to seismic data containing primaries, free-surface correlated multiples, and inter-layer multiples, excluding direct arrivals and ghost waves. Therefore, for actual seismic data, if direct arrivals and ghost waves are present in the input data, direct arrival removal and ghost wave removal preprocessing should be performed before multiple wave suppression processing.

[0021] S2: This application transmits the seismic data obtained in step one and resides it in the GPU memory space. Internally on the GPU hardware, all matrix multiplications and Fourier transforms are mapped to multi-threaded, high-concurrency processing. Simultaneously, this application cleans up intermediate variables and releases memory in real-time during the computation process, which is beneficial for handling large-scale 3D seismic data.

[0022] S3: This application is based on a feedback model. It includes total input seismic data in the space-frequency domain, encompassing all orders of surface-correlated multiples. With the Green function to be optimized and wavelets The following bilinear equation applies between them: (1) In the formula, The representative free surface reflection operator (taking the negative value of the identity matrix under the simplified assumption of the air-water interface reflection coefficient, i.e.) The first term on the right side of equation (1) Defined as a primary wave signal, the second term This represents the full-order surface correlation multiple wave. In the implementation of this application, this application introduces in the forward modeling operator... Factors, to meet technical requirements, namely This should represent the normal derivative of Green's function with respect to the free surface. This treatment is consistent with the obliquity factor used in EPSI and REPSI. Angular frequency is defined as follows: , For frequency.

[0023] The main task of this application is to solve equation (1) from observation data. Estimating two unknowns and It can be seen that equation (1) in and The above is bilinear. In other words, when fixed... and At one time, the equation degenerates into a linear form with respect to another variable. This formulation is similar to the traditional blind deconvolution problem, which is typically affected by amplitude scale uncertainties and time / phase translation uncertainties. However, the existence of free-surface correlated multiples helps mitigate these uncertainties.

[0024] Equation (1) contains two unknowns to be solved. and To solve this equation, the sparsity of the Green's function in the time domain needs to be taken as a prerequisite. Equation (1) can be reconstructed as: (2) in, This is the vectorized form of the total input seismic data in the time domain. Let be the vector of the Green's function to be determined in the time domain. This is a forward modeling operator, which includes composite operations such as time-frequency transformation, zero-value filling, and matrix multiplication. Given the wavelet estimated in the current step... Then, equation (2) becomes a solution for the impulse response that can be obtained using an inversion algorithm. A system of linear equations.

[0025] Due to Green's function The objective function exhibits a sparse distribution in the time domain. This application inverts the Green's function by solving a minimization problem based on L1 norm constraints. To avoid using the complex SPGL1 solver of traditional methods, this application innovatively introduces the ALB algorithm. The objective function is: (3) In equation (3), For sparse threshold parameters, This indicates the error tolerance. The ALB algorithm specifically includes the following steps: a) Obtain the following parameters: forward operator Observation data vector Initial value vector threshold Maximum number of iterations ; b) Initialization: Iteration variables , Will Assign to ; c) Repeat steps d through f below. Second-rate; d) Calculation ; e) Calculation ; f) will Updated to +1; g) Obtain the sparse solution = This sparse solution This refers to the update performed using the linear Bregman algorithm. The result was as follows. same , which is the initial solution; l Index for iteration count; , These are temporary intermediate variables during the iteration process; This is the step size parameter.

[0026] S4: Obtain the estimated value of the three-dimensional Green's function matrix through step three. Subsequently, this application will Treating it as a known quantity, for the wavelet The solution is then performed. To avoid the high computational burden caused by multiple iterations of the least squares algorithm for wavelet solving in the traditional EPSI method, this application utilizes multiple wave terms in the feedback model to overcome the amplitude and phase ambiguity problem in blind deconvolution, and derives the wavelet solution formula: (4) In equation (4), the superscript Represents the complex conjugate transpose of a matrix; This is to predict surface multiple wave signals using the current best sparse Green's function and input data. Using this formula, at each frequency... The wavelet in the current loop is then solved.

[0027] S5: Repeat steps three and four, alternating between updating the Green's function. With wavelet After reaching the set number of iterations or error tolerance, the three-dimensional surface-related multiple waves and primary waves are finally separated.

[0028] The following are specific implementation examples provided in this application.

[0029] Figure 1 The flowchart illustrates the three-dimensional seismic surface multiple wavelet and wavelet estimation method based on the accelerated linear Bregman algorithm and GPU acceleration, as shown in this embodiment. In fact, the main steps of 3D EMS are essentially the same as those of 2D EMS. Figure 1 A flowchart of the proposed 3D EMS method is shown. Similar to the two-dimensional case, the algorithm estimates the Green's function and wavelet without a free surface in an alternating manner: first, in a fixed... Optimization under the condition Subsequently, in a fixed Optimization under the condition Both are initialized to zero. After obtaining a sparse estimate of the Green's function without free surfaces, a wavelet estimation step is performed, because problems with sparse constraints usually yield sparse solutions upon convergence. Consistent with the 2D EMS method, this application also employs a frequency downsampling strategy, significantly reducing computational costs without compromising overall accuracy. Examples include the ALB algorithm, Equation 4, and... Figure 1 As shown, the overall process of the 3D EMS method remains relatively simple and does not involve complex algorithms, which is one of the important advantages of this method.

[0030] To verify the method of this application, two representative three-dimensional geological models were used in the embodiments. For example... Figure 2 As shown, the first model is a simple three-layer horizontal velocity model, with spatial sampling intervals of 6.25 m in the x, y, and z directions. The first layer has a longitudinal wave velocity of 1800 m / s and a density of 1000 kg / m³; the second layer has a velocity of 2300 m / s and a density of 3000 kg / m³; and the third layer has a velocity of 2000 m / s and a density of 1100 kg / m³. Figure 3As shown, in order to test complex three-dimensional geological conditions, a three-dimensional Overthrust velocity model was further used, with a model size of 201×201×201 grid points.

[0031] It should be noted that in the simulations of both datasets, the source and detector were located on a free surface (depth of 0 meters). Therefore, the simulated data did not contain ghost waves, and no ghost wave removal processing was required. As shown in the subsequent results, the direct wave was also completely removed before the multiple wave removal processing was performed. Specifically, the direct wave was simulated using a uniform background model and then directly subtracted from the simulated data generated based on the real model. In other words, this application intentionally does not consider preprocessing steps such as ghost wave removal and direct wave removal in order to focus the research on the goal of removing multiple waves. Cases involving approximate preprocessing errors (e.g., ghost waves or direct waves not being accurately eliminated) are not discussed in this application.

[0032] Figure 4 for Figure 2 The model shown simulates the total shot gather, with the source located at the model center at depth z = 0. The wavelet used in the forward modeling of the wave equation is a 15 Hz Ricker wavelet, with a time sampling interval of 1 ms. The detector spacing is 25 m. The figure includes not only the effective primary reflection wave but also strong surface-correlated multiples. This figure visually illustrates the state of the original 3D data and the complex interference of the multiples.

[0033] for Figure 2 The three-dimensional horizontal layered model shown in this application does not employ conventional methods to simulate all the shot sets required to construct the total dataset D. To reduce the computational cost of forward modeling of the wave equation, only a shot set with a large spatial coverage was generated at the center of the model. Figure 4 Subsequently, based on the source-detector reciprocity principle, a complete total dataset D is constructed by translating the source position. Compared with shot-by-shot forward modeling, this strategy can significantly improve computational efficiency and reduce computational costs.

[0034] After acquiring the raw 3D data, the data volume, which typically has 5 dimensions (time, shot point x and y, receiver point x and y), needs to be reshaped into a 3D input matrix. (Arranged according to time axis, common receiver axis, and common shot axis).

[0035] Figure 5 The results of the extraction of the three-dimensional common shot point gathers for the first, middle, and last layers of the recombined three-layer horizontal model are presented.

[0036] Figure 6 Seismic data for a three-level horizontal model processed using the method of this application The time slice (at 0.5 seconds) used a total of 1681 sources and 1681 detectors (i.e., 41 main survey lines, each with 41 intersecting survey points). Theoretically, this matrix is ​​a square matrix with symmetry, its diagonals corresponding to zero-offset data. This structure reflects the three-dimensional effects produced by three-dimensional wavefield propagation; in contrast, the two-dimensional case corresponds to only a symmetrical square matrix.

[0037] The preprocessed seismic data is transmitted and stored in the GPU memory. Within the constructed GPU computing architecture, intensive multidimensional matrix multiplication and fast Fourier transform tasks are all executed concurrently by multi-threaded GPUs, improving computational efficiency.

[0038] Figure 7 for Figure 2 The results of processing the center resonance source gather of the three-layer horizontal model shown in this application are obtained using the method described in this application. Figure 10 A series of figures illustrate the comparison of the processing effects of the proposed method on the central resonance source gather. (a) shows the original input shot gather containing the primary wave and all multiples; (b) shows the three-dimensional surface-related multiples accurately predicted and separated by the proposed method; and (c) shows the clean shot gather record output after removing surface multiples. The comparison demonstrates that the proposed method can separate multiples with high fidelity in three-dimensional space and effectively recover primary wave information.

[0039] Figure 8 This shows the wavelet information estimated from the simulation data of the three-level model using the method of this application. The estimated wavelet of the three-level model dataset. (a) Time-domain waveform; (b) Amplitude spectrum. The gray curve represents the original estimated wavelet, and the black curve represents the estimated wavelet after correcting for the obliquity factor.

[0040] Figure 9 for Figure 7 The central coherent source gather of the three-layer horizontal model shown is the central shot gather. To provide clearer local details, Figure 9 A series of images were extracted Figure 7 The centerline of the coherent source gather is displayed. Figure 7Central common shot record of data and corresponding processing results. (a) Original common shot record; (b) Free surface correlated multiples estimated and eliminated by 3D SRME method; (c) Primary wave result obtained by 3D SRME processing; (d) Difference between Figure (c) and Figure (f); (e) Free surface correlated multiples estimated and eliminated by 3D EMS method; (f) Primary wave result obtained by 3D EMS processing. This figure further confirms that even in areas with strong multiple interference, the method of this application can achieve complete multiple elimination without damaging the amplitude and phase axis morphology of the effective primary wave at depth. The method of this application is superior to the "prediction-subtraction" type 3D SRME method. The surface correlated multiple elimination of the method of this application is cleaner. To verify the advantages of 3D EMS, this application independently implemented 3DSRME based on the methods proposed by Verschuur et al. (1992) and Berkhout and Verschuur (1997). The results of 3D SRME are compared with those of 3D SRME. Figure 9 Compared to (c), the results of 3D EMS ( Figure 9 The conclusion that (f) can more effectively suppress surface-correlated multiples can also be drawn from the difference plot (f) Figure 9 This is further confirmed in (d). Therefore, these results demonstrate that 3D EMS outperforms 3D SRME in terms of multiple wave cancellation accuracy, and can obtain more accurate multiple wave cancellation results.

[0041] Although both SRME and EMS are based on the same feedback model (Equation 1), they differ fundamentally in their implementation mechanisms. SRME is a typical prediction-subtraction method, and its core idea is to use input data... Approximate replacement of Green's function without free surface And construct R based on this. As a surface-correlated multiples prediction model, this approximation leads to prediction errors, and the subsequent adaptive subtraction process based on the least squares criterion often fails to completely eliminate residual multiples. In contrast, EMS directly inverts the wavelet by jointly solving the feedback model through ALB sparse inversion and one-step scaling estimation. Green's function without free surface Therefore, the constructed multiple wave prediction term R The correlation multiples of a more realistic free surface can theoretically provide more accurate multiple prediction and elimination effects than SRME.

[0042] Figure 10 for Figure 3The figure shows a 3D view of the total shot gather simulated by the model, with the source located at the center of the model at depth z = 0. The wavelet used in the forward modeling of the wave equation is a 13 Hz Ricker wavelet, with a time sampling interval of 0.5 ms. The detector spacing is 25 m. The figure contains complex primary reflection waves and strong surface-correlated multiple waves, visually demonstrating the severe interference caused by the interlacing, overlapping, and out-of-plane scattering of multiple wave fields under complex 3D undulating structures.

[0043] Figure 11 The results of the extraction of the first, middle and last three-dimensional common shot point gathers of the recombined Overthrust model are presented.

[0044] Figure 12 The time slice of the seismic data from the Overthrust model (at 0.5 seconds) uses a total of 1681 sources and 1681 geophones (i.e., 41 main survey lines, each with 41 intersecting survey points). This figure visually reflects the non-uniform distribution of the wavefield caused by complex subsurface structures.

[0045] Figure 13 for Figure 3 The results of processing the central resonance source gather of the Overthrust model using the method described in this application. Figure 13 A series of figures illustrate the comparison of the processing performance of the proposed method under complex structures. (a) shows the original input shot gather containing the primary wave and all multiples; (b) shows the three-dimensional surface-related multiples accurately predicted and separated by the proposed method; and (c) shows the clean shot gather record output after eliminating surface multiples. The comparison demonstrates that even under conditions of strong subsurface structural undulations and complex wavefield intersections, the proposed method can still separate multiples with high fidelity within a three-dimensional framework.

[0046] Figure 14 This shows the wavelet information estimated from simulation data of the 3D Overthrust model using the method of this application. (a) Time-domain waveform; (b) Amplitude spectrum. The gray curve represents the original estimated wavelet, and the black curve represents the estimated wavelet after correcting for the obliquity factor.

[0047] Figure 15The results of center line shot gather processing on a complex Overthrust 3D model were compared using a 2D surface multiples elimination method and the 3D surface-correlated multiples elimination method of this application. (a) to (c) show the processing results of the 3D surface-correlated multiples elimination method of this application, (d) shows the same original record, (e) shows the multiples removed by the 2D surface multiples elimination method, and (f) shows the primary wave after processing by the 2D surface multiples elimination method. This figure clearly reveals that due to the existence of 3D out-of-plane effects, conventional 2D algorithms have prediction errors in both dynamic amplitude and kinematic travel time, resulting in multiples remaining and damage to the deep primary wave signal, representing incomplete "suppression." This application, however, achieves the elimination of surface-correlated multiples in 3D seismic data.

[0048] For 3D seismic data acquired under complex subsurface structural conditions, the prediction of free-surface correlated multiples should be performed within a complete 3D framework to accurately describe its 3D wavefield propagation characteristics. This is also the fundamental premise of 3D EMS. For comparison, this application also applies 2D EMS to the same 3D dataset, and the central shot processing results are as follows: Figure 15 As shown in (b)–15 (c). Comparison with 3D EMS results ( Figure 15 Compared to (f), the results of 2D EMS ( Figure 15 In (c) there are more residual multiples, especially in the deep region, and this difference becomes more pronounced with the enhancement of three-dimensional effects. Because the two-dimensional approximation cannot accurately characterize the true three-dimensional wavefield propagation process, 2D EMS produces deviations in the prediction of amplitude and travel time of free-surface correlated multiples, thus achieving only a limited degree of multiple suppression and failing to achieve true multiple elimination. In contrast, 3D EMS can more accurately predict and reconstruct free-surface correlated multiples in the original data, thereby achieving more thorough multiple elimination and obtaining clearer imaging results. These results demonstrate that, under complex three-dimensional geological conditions, fully considering three-dimensional effects is crucial for achieving high-precision free-surface correlated multiple elimination, and further prove the necessity of developing the 3D EMS method.

[0049] Figure 16The center-line shot gather processing results on a complex Overthrust 3D model were compared using the 3D SRME method and the 3D surface-related multiples elimination method of this application. (a) to (c) show the processing results of the 3D SRME surface-related multiples elimination method; (b) shows the free surface-related multiples eliminated using the 3D SRME method; (c) shows the first-order wave result obtained by 3D SRME processing; (d) shows the difference between Figure (c) and Figure (f); (e) shows the free surface-related multiples eliminated using the 3D EMS method; and (f) shows the first-order wave result obtained by 3D EMS processing. This application compares the processing effects of 3D EMS and 3D SRME (see [link to application]). Figure 16 and Figure 17 Compared with the results of 3D SRME (Figure 16(c)), the results of 3D EMS ( Figure 16 (f) can more effectively suppress residual free surface correlated multiples. The multiple energy predicted by 3D SRME is weaker ( Figure 17 (b) reveals the inherent limitations of the “prediction-subtraction” method. Furthermore, in practical applications, adaptive least squares subtraction often fails to completely remove correlated multiples on free surfaces, which further validates the advantages of 3D EMS over 3D SRME.

[0050] Figure 17 This paper presents a complete zero-shot-receiver distance data volume profile extracted from the total data volume of the 3D Overthrust model. Figure 17 (a) shows the original noisy data. The arrows indicate locations with numerous strong 3D surface-dependent multiples caused by complex seafloor undulations and deep structures. (a) Zero-offset data volume extracted from the total data volume D of the 3D Overthrust model; (b) Free surface-dependent multiples estimated and eliminated using the 3D SRME method; (c) Zero-offset data of the first wave after 3D SRME processing; (d) Differences between Figure (c) and Figure (f). All subplots use the same amplitude display range. Arrows indicate areas where free surface-dependent multiples are effectively eliminated. Figure 17 (e) represents the three-dimensional surface multiple wave noise data volume estimated by 3D EMS using the method of this application. Figure 17 Figure (f) shows the final processed zero-shot-receiver distance pure primary wave data. By comparing the results with Figure (f), it can be seen that the strong surface-dependent multiple wave interference in the area indicated by the arrow has been cleared, and the deep primary wave signal has been recovered.

[0051] This application further compares the processing effects of 3D EMS and 3D SRME (see...). Figure 16 and Figure 17 ). Result with 3D SRME ( Figure 16Compared to (c), the results of 3D EMS ( Figure 16 (f) can more effectively suppress residual free surface-related multiples. Meanwhile, the energy of the multiples model predicted by 3D SRME is significantly weaker (f) Figure 17 (b) This reflects the inherent limitations of the prediction-subtraction method, namely that the multiple model obtained by approximate prediction is difficult to accurately represent the real free surface correlated multiples. Furthermore, even with adaptive least squares subtraction, it is often difficult to completely eliminate residual multiples in actual data. Therefore, 3D EMS outperforms 3D SRME in both multiple prediction accuracy and final elimination effect, further demonstrating its advantage over the traditional prediction-subtraction method.

[0052] Figure 16 and Figure 17 This further demonstrates that the 3D EMS method of this application outperforms the "prediction-subtraction" type 3D SRME method. The 3D EMS method of this application achieves cleaner surface correlation multiple wave elimination.

[0053] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for estimating multiple waves and wavelets of a three-dimensional surface based on accelerated linear Bregman and GPU, characterized in that, Includes the following steps: S10: Acquire the original three-dimensional seismic input data containing surface-correlated multiples; reorganize the original three-dimensional seismic input data in the dimensions of time, common receiver points, and common shot points to obtain a time-domain seismic wavefield; perform a one-dimensional fast Fourier transform on the time-domain seismic wavefield along the time axis to obtain a three-dimensional data matrix in the frequency domain; and vectorize the original three-dimensional seismic input data to obtain the corresponding time-domain vectorized data; transfer the three-dimensional data matrix and the time-domain vectorized data to the GPU memory. S20, on the GPU, using the three-dimensional data matrix and the time-domain vectorized data as input, and under the condition of fixed wavelets, based on the feedback model expression, the sparse inversion problem is solved using the accelerated linear Bregman algorithm to obtain a sparse estimate of the Green's function without surface multiple interference; the feedback model expression describes the relationship between the total input seismic data containing all orders of surface-correlated multiples, the Green's function to be optimized, and the wavelets. S30, under the condition that the sparse estimate of the Green's function is fixed, the corresponding wavelet is directly estimated using the analytical solution formula of the wavelet based on the sparse estimate of the Green's function. S40, repeat S20 and S30, alternately update the sparse estimate of the Green's function and the wavelet until the set number of iterations or error tolerance is reached and the iteration stops. Based on the sparse estimate of the Green's function and the wavelet obtained when the iteration stops, separate the surface-related multiples and primarys from the original three-dimensional seismic input data.

2. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 1, characterized in that, In step S20, the feedback model expression is a bilinear equation. ,in This represents the three-dimensional data matrix. Denotes the Green's function, Indicates the sub-wavelength, Represents the free surface reflection operator. Indicates a primary wave signal. Represents the full-order surface-correlated multiple wave; the bilinear equation in and When one variable is fixed, the free surface reflection operator degenerates into a linear form with respect to the other variable; Under the simplified assumption of the air-water interface reflectance coefficient, the negative value of the identity matrix is ​​taken.

3. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 1, characterized in that, In step S20, the step of using the accelerated linear Bregman algorithm to solve the sparse inversion problem specifically includes: The feedback model expression is reconstructed into a system of linear equations in the time domain. ,in For the time-domain vectorized data, Let the vector of the Green's function to be determined in the time domain be . The forward operator includes time-frequency transformation, zero-filling, and matrix multiplication operations. Based on the sparse distribution of the Green's function in the time domain, the Green's function is inverted by solving a minimization problem based on L1 norm constraints. The objective function is: ,in Denotes the L1 norm of g. express The L2 norm error between d and d For sparse threshold parameters, Indicates the error tolerance.

4. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 3, characterized in that, The specific steps of the accelerated linear Bregman algorithm are as follows: Obtaining forward operators Observation data vector Initial value vector threshold Maximum number of iterations ; Initialization: Iteration variables , Will Assign to ; Repeat the following steps Second-rate: calculate ; calculate ; Will Updated to +1; Derive sparse solutions = ;in, same , which is the initial solution; l Index for iteration count; , These are temporary intermediate variables during the iteration process; This is the step size parameter.

5. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 1, characterized in that, In step S30, the analytical solution formula for the wavelet is as follows: superscript H Represents the complex conjugate transpose of a matrix. This represents the spectrum of a certain frequency component of the sparse estimate of the Green's function obtained through step S20. This represents the spectrum of a certain frequency component of the current estimated wavelet. This represents the spectrum of a certain frequency component of the input data. This represents the spectrum of a surface multiple signal at a certain frequency, predicted using a sparse estimate of the current best Green's function and the input data; using this formula, at each frequency... The wavelet in the current iteration loop is then solved.

6. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 1, characterized in that, In step S10, the time axis of the original three-dimensional seismic input data Coordinates of the main survey line of the receiver point Coordinates of the detection point connection line Main survey line coordinates of the shot point and the coordinates of the shot point connection line The original 3D seismic input data is reshaped using a seismic line index, and its dimensionality is reduced and reorganized into a 3D input data matrix. ,in This represents the number of time sampling points; Total number of detector points ( ); Total gun points ( Under ideal sampling conditions, the shot points and receiver points are arranged to completely overlap, meaning the total number of shot points equals the total number of receiver points. Reorganized time-domain seismic wave field After performing a one-dimensional fast Fourier transform along the time axis, it is transformed into a three-dimensional data volume in the frequency domain. ,in This represents the total number of discrete frequency components.

7. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 1, characterized in that, In step S10, after acquiring the original three-dimensional seismic input data, if the original three-dimensional seismic input data contains direct waves and ghost waves, the original three-dimensional seismic input data is preprocessed by removing direct waves and ghost waves before performing reconstruction and fast Fourier transform.

8. The method for estimating three-dimensional surface multiple waves and wavelets based on accelerated linear Bregman and GPU as described in claim 1, characterized in that, In step S20, within the GPU hardware, all matrix multiplications and Fourier forward and inverse transformations are mapped to multi-threaded high-concurrency processing, and intermediate variables are cleaned up in real time during the calculation process to release video memory.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method as described in any one of claims 1 to 8.

10. A computer-readable storage medium storing computer instructions thereon, characterized in that, When executed by a processor, the computer instructions implement the steps of the method as described in any one of claims 1 to 8.

Images