Method for suppressing marchenko multiple based on wasserstein distance correction of wave field mismatch

Abstract

The present application belongs to the technical field of seismic data processing in geophysical exploration, and relates to a Marchenko multiple wave suppression method based on Wasserstein distance correction of wave field mismatch, comprising: obtaining a pulse reflection sequence through preprocessing and deconvolution; inputting the pulse reflection sequence to perform Marchenko iteration initialization; and remodeling the Marchenko iteration based on the Wasserstein distance correction of the wave field mismatch. The present application first systematically introduces the optimal transport theory into the field of Marchenko multiple wave suppression, while maintaining the advantage of the Marchenko method that it does not require a velocity model, significantly improves the adaptability of the method to actual data noise, defects and amplitude distortion, thereby realizing higher-precision multiple wave suppression and more reliable primary wave recovery, and is used for high-precision seismic imaging in a complex onshore exploration environment. The present application has significant advantages in constraint mechanism, error control and global feature protection.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of seismic data processing in geophysical exploration, and particularly relates to a Marchenko multiple wave suppression method based on Wasserstein distance constraint of optimal transport theory, and particularly relates to a Marchenko multiple wave suppression method based on Wasserstein distance correction of wave field mismatch. BACKGROUND

[0002] In seismic exploration, the underground structure is reconstructed mainly by analyzing the response of artificially excited seismic wave field, and the effectiveness and accuracy of the reconstruction process directly depend on the correct identification and understanding of different components in the wave field. Among them, the primary wave as the effective signal in the seismic wave field refers to the wave that is emitted from the source, reflected once by the underground reflection interface and then returned to the ground surface to be received. Because of its unique and direct propagation path, the travel time, amplitude and phase information it carries can truly and accurately reflect the depth, shape and wave impedance characteristics of the underground reflection interface, and the primary wave is also the core basic data for imaging and reservoir inversion. In contrast, the multiple wave is inherent interference noise in the seismic wave field, which is formed by multiple reflections between the wave and multiple strong wave impedance interfaces in the underground, such as the sea bottom, the basement or the main lithology interface or between the ground surface. This complex propagation mechanism results in a non-unique path, and the travel time overlaps with the primary wave, eventually forming false events on the seismic record that do not reflect the true stratigraphic structure. These events not only produce structural illusions in the seismic section and mislead geological interpretation, but more seriously, when the multiple wave and the primary wave are mixed in the time-space domain, the dynamics characteristics such as amplitude, frequency and phase of the effective wave will be distorted, resulting in inaccurate results of reservoir fluid identification and quantitative interpretation based on these characteristics.

[0003] Studies have shown that the essential difference between the primary wave and the multiple wave lies in the uniqueness and authenticity of their propagation paths. How to accurately separate the pure primary wave field in seismic data processing is the core prerequisite for improving the reliability of seismic exploration and achieving accurate oil and gas reservoir identification. In the prior art, two types of multiple wave suppression methods are mainly used. One is the filtering method based on the kinematic difference of wave field, such as Radon transform, F-K filtering, etc., which has a simple principle but has limitations: for interlayer multiple waves that are highly overlapped with the primary wave velocity spectrum, filtering will seriously damage the effective signal while removing the interlayer multiple waves. The second is the method based on wave equation theory, such as CFP method, ISS method, etc. Although the above methods can better describe the physical process of wave propagation, they face the problem of suppressing interlayer multiple waves while suppressing precision or requiring accurate prior model velocity information.

[0004] To this end, researchers have developed a Marchenko multiple elimination method based on Marchenko imaging theory. This method introduces an integral operator independent of the model to extrapolate the wave field to the observation surface, and theoretically completely eliminates the dependence on the velocity model. However, Marchenko multiple elimination method faces a series of practical challenges when moving from theory to practical application: Marchenko multiple elimination method needs to follow several assumptions. First, the method is based on ideal input data without noise interference; second, the theoretical derivation requires an infinite dense sampling point distribution of the observation system; and third, the algorithm design assumes that the source wavelet is an ideal pulse and the reflection coefficient sequence is an ideal pulse sequence. At the same time, a very important assumption is that the medium of the assumed data in Marchenko is an acoustic medium, and there is no attenuation in the medium. However, there is a significant gap between the above theoretical settings and the actual acquisition of land seismic data. First, the noise sources in the actual field observation environment are complex and have strong energy, especially the strong amplitude interference of surface waves is widely distributed. Even after systematic noise suppression processing, it is still difficult to completely eliminate all noise components. Second, due to the restriction of complex surface conditions and acquisition cost, the spatial sampling density of the actual observation system is often difficult to reach the theoretical ideal value required by the method, resulting in incomplete recording of wave field information. In addition, the actual source wavelet has obvious time domain continuation characteristics, which is essentially different from the ideal pulse assumption, which also makes the reflection response unable to converge to an ideal pulse sequence form. Moreover, a very important practical obstacle is that the actual acquisition medium is attenuated, which does not completely meet the Marchenko assumption of no attenuation of acoustic waves.

[0005] At present, the existing land-based measured data Marchenko multiple wave suppression method based on compressed sensing can make the noise level of the input actual data lower than that of the data after conventional denoising through sparse transformation from the perspective of high-quality preprocessing, the collection density is higher, the wavelet is approximately pulse sequence, and finally the input data is as close as possible to the Marchenko input assumption to realize the preliminary practicality. However, even if it is close, it is still difficult to achieve perfect assumption. Especially for attenuating media, even after a certain amplitude compensation, it is still difficult to restore the non-attenuating amplitude due to the complexity of the actual attenuation. These imperfections will cause errors in the Marchenko element-by-element corresponding iterative paradigm, and the errors will be treated as effective wave field values, participating in the subsequent wave field construction, and being gradually amplified in the iteration, and finally, leading to poor inversion effect or even false appearance. Therefore, it is not complete to improve the Marchenko method only from the perspective of preprocessing to realize the practicality. It is necessary to build a global constraint iterative mechanism with error suppression capability. The mechanism needs to meet the dual goals: on the one hand, it needs to effectively control the accumulation and transmission of errors in the iteration process; on the other hand, it must maintain the key physical properties of the wave field, including the macroscopic geometric structure of the same phase axis, such as continuity, curvature and the relative relationship of the amplitude.

[0006] The existing regularization method based on point-by-point comparison, such as Tikhonov regularization or sparse constraint method, has obvious limitations in this problem. Although the L2 norm constraint can effectively suppress noise interference, the uniformization process will make the wave front profile fuzzy and weaken the kinematic characteristics of the wave field; while the constraint based on L1 or total variation norm can better maintain the signal edge, but it will change the dynamics of the wave. The common limitation of the above methods is that the mathematical essence of point-by-point comparison cannot accurately distinguish the local amplitude error caused by noise from the global structural change caused by the kinematic characteristics of the wave field, such as time shift, stretching, etc. Therefore, additional wave distortion is often introduced in the noise suppression process.

[0007] In summary, there is an urgent need to develop a Marchenko multiple wave removal method that can simultaneously maintain the kinematic characteristics of the wave field and maintain the dynamics of the global constraint mechanism. This method should not only effectively suppress the error propagation in the iteration process, but also completely protect the key physical information contained in the wave field from the global perspective, so as to realize stable and high-fidelity multiple wave suppression. SUMMARY

[0008] The purpose of the present application is to provide a robust interlayer multiple wave suppression method for land measured seismic data by introducing a Wasserstein distance constraint mechanism based on optimal transmission theory, which introduces Wasserstein distance as a matching criterion and can suppress error propagation, so as to solve the problem of error propagation sensitivity of the existing Marchenko multiple wave suppression method in the application of actual land seismic data.

[0009] The present application is realized by the following technical solutions:

[0010] A Marchenko multiple wave suppression method based on Wasserstein distance correction of wave field mismatch, comprising the following steps:

[0011] a: obtaining the impulse reflection sequence by preprocessing and deconvolution:

[0012] First, the original collected seismic data is reorganized into a three-dimensional data body, and a regular grid is established; then, the data is converted to a sparse domain by using a Shearlet transform with multi-scale and multi-direction analysis capability; then, hard threshold denoising processing and wave field reconstruction are sequentially performed in the Shearlet domain to suppress noise and repair missing channel data; finally, high-resolution impulse reflection response is obtained by sparse constraint deconvolution technology ;

[0013] b: input the impulse reflection sequence, and perform Marchenko iteration initialization:

[0014] The impulse reflection sequence obtained in step a is used as input to define the initial wave field, set the initial primary wave field , set the initial multiple wave prediction term , set the initial iteration number , set the maximum iteration number , and the iteration process is cyclically performed under the condition of , and is divided into odd and even iterations; wherein, the odd iteration updates the primary wave field: the original wave field minus the predicted multiple wave field; the even iteration updates the predicted multiple wave field: the time window term acts on the last even predicted multiple wave field;

[0015] c: remodeled Marchenko iteration based on Wasserstein distance correction of wave field mismatch:

[0016] c1: remodeled Marchenko odd iteration based on Wsserstein distance correction of predicted wave field amplitude and waveform:

[0017] ​c11: Standardization is performed on the input data of the original wave field and the predicted multiple wave field. First, non-negative standardization is performed on the original seismic trace to eliminate negative amplitudes; then amplitude normalization is performed;

[0018] c12: A corresponding target function is set to find an optimal transport mapping to minimize the total transportation cost:

[0019] c13: The space-time-amplitude index is extracted, and a space-time-amplitude composite cost function is designed:

[0020] c14: Sinkhorn iteration is used to solve the optimal transport mapping;

[0021] c15: The subtracted multiple wave field is updated;

[0022] c2: Reshaped Marchenko even iteration based on Wsserstein distance correction of the predicted wave field phase:

[0023] c21: Standardization is performed on the input data of the time window item and the last predicted multiple wave field, and non-negative standardization is performed on each trace;

[0024] c22: Exponential coding enhances the phase feature;

[0025] c23: The CDF inverse function of the standardized time window item and the multiple wave field is calculated to obtain the candidate time shift;

[0026] c24: The optimal time shift mapping that minimizes the Wasserstein cost function is solved, and the predicted multiple wave field is updated;

[0027] c3: Marchenko iteration based on Wsserstein distance correction of wave field mismatch is completed through odd-even cycle iteration, and the primary wave field is finally output.

[0028] Further, step a, high-resolution pulse reflection response is obtained by the following formula:

[0029] ;

[0030] wherein, is the original observation data, is a sampling operator representing the data missing position, is a source wavelet convolution matrix, represents the convolution process, is an unknown reflection response to be solved, is a reflection coefficient sparse term, is a Shearlet transform, represents a seismic record sparse term, and denote the L1 and L2 norm, respectively, with is a regularization parameter;

[0031] data fidelity term to ensure that the processed result can fit the original observation data ; is a reflection coefficient sparsity term, which is constrained by the L1 norm to make the solution itself sparse in the specific Shearlet domain; denotes the seismic record sparsity term, which constrains the synthesized seismic record to be also sparse in the Shearlet domain, and simultaneously realizes effective wavefield reconstruction and noise suppression.

[0032] Further, step b, the initial primary wavefield and multiple wavefield are defined: ; ;

[0033] wherein, is a specific time window, which is used to cut off the part of the data outside the time window is the initial primary wavefield, is the initial multiple wavefield; the odd-numbered iteration mainly performs multiple wave subtraction, that is, subtracts the multiple wavefield predicted in the last even-numbered iteration from the primary wavefield in the last odd-numbered iteration, and the traditional Marchenko multiple wave removal method is described as:

[0034] ;

[0035] wherein, is the updated primary wavefield in the current odd-numbered iteration, is the updated primary wavefield in the last odd-numbered iteration, is the multiple wavefield predicted in the last even-numbered iteration;

[0036] the even-numbered iteration, is even, and mainly performs multiple wave prediction:

[0037] ;

[0038] wherein, is the predicted multiple wavefield in the current even-numbered iteration, is the predicted multiple wavefield in the last even-numbered iteration, and denote the cross-correlation and convolution forms of respectively, and for the convenience of representation, is​​ ;

[0039] The conventional Marchenko iteration is to perform the above-mentioned even-odd iteration loop, when the iteration stops, and all time slices of updated in the last iteration are collected as the final one-way wavefield.

[0040] Further, the step c11 specifically comprises: normalizing the two-dimensional data participating in the operation, denoted as , , wherein respectively represent the number of time sampling points and the number of spatial channels, and the normalization process is written in the form of discrete points as follows:

[0041] ;

[0042] ;

[0043] , , denotes the minimum value in , denotes the maximum value in , denotes the minimum value in , denotes the maximum value in , is a very small positive number to prevent zero.

[0044] Further, the step c12, the total transmission cost function is expressed as follows:

[0045] ;

[0046] , wherein and total transmission cost function, and transmission cost function of a single point .

[0047] Further, the step c13 specifically comprises the following steps:

[0048] , , and , :​​​​​​

[0049] ;

[0050] wherein, and are respectively and time sampling point indexes, and are respectively and spatial channel number indexes, and are respectively and amplitude sizes corresponding to time and space indexes, is the number of time and space indexes, the higher the number means the greater the gradient of the subsequent space-time cost function, so that the amplitude correction will be limited to a very small range;

[0051] controlling the time and space range in the composite cost function term is determined as:

[0052] ;

[0053] controlling the amplitude in the composite cost function term is determined as:

[0054] ;

[0055] the composite cost function is determined as:

[0056] ;

[0057] and are two weight constants for controlling the space-time cost and the amplitude cost, which are artificially set.

[0058] Further, the step c14 specifically includes the following steps:

[0059] an entropy penalty term is introduced into the objective function, and the problem is converted into the following form and then solved efficiently by Sinkhorn iteration:

[0060] ;

[0061] wherein, is the entropy of the transmission plan , and is a regularization coefficient. controlling the transmission mapping matrix The sum of each row element of must equal the corresponding component of the source distribution un, and the sum of each column element of must equal the corresponding component of the target distribution mn, where is a column vector with all elements equal to 1;

[0062] The Sinkhorn iteration is used to solve the problem, as follows:

[0063] ;

[0064] ;

[0065] where, is the intermediate constructed matrix for computing the optimal transport mapping , which is completely determined by the cost function . The larger the cost function is, the smaller the value of K corresponding to is. In the equation, and are optional scaling vectors, which are used to ensure that satisfies the marginal constraint of mass conservation, denotes a diagonal matrix form, denotes ;

[0066] The iteration process is as follows:

[0067] First, initialize ; then, expand the alternating iteration for :

[0068] ;

[0069] ;

[0070] where, denotes element-wise division, and denote the scaling vectors updated in the current iteration;

[0071] Iterate until the changes of and are smaller than a certain threshold, and the final solution of the optimal transport mapping is:

[0072] ;

[0073] For this optimal transport mapping , the optimization of the predicted multiple wavefield is:

[0074] ;

[0075] where, is an artificially set constant that controls the degree of update, Indicates the direction of transport, i.e., whether the amplitude decreases or increases; the subscript here... and express and The element index in the middle, For the optimally mapped The elements in the corrected standardized multiple wave field, and the overall updated multiple wave field, need to be inversely standardized to be restored to a state where multiple wave subtraction can be performed;

[0076] Step c15 specifically includes the following steps:

[0077] For each of the collected of Multiple wave fields after overall elemental correction Perform inverse standardization:

[0078] ;

[0079] This refers to the multiple wave field subtracted from the final corrected multiple wave field.

[0080] Further, step c21 specifically includes the following steps:

[0081] In even-numbered iterations, the phase-matching correction process is performed track-by-track, and the selected object is... and The continuous function forms of the single-channel records in the corresponding time slices are as follows: and Similarly, because the signal amplitude in seismic data fluctuates (positive and negative) and the total energy is uncertain, direct application of Wasserstein distance calculation requires prior non-negative data standardization.

[0082] ;

[0083] ;

[0084] in, To prevent extremely small positive numbers from being divided by zero, for The minimum value in, for The minimum value in, for The maximum value in the above standardization preprocessing and The two corresponding standardized single tracks and All are non-negative and have the same measure. After non-negative normalization, the original Dao unifies the sums to... of the standard.

[0085] Further, the step c22 specifically comprises the following steps:

[0086] After standardization, the target function of the even times of iteration amplitude correction is defined as the Wasserstein distance:

[0087] ;

[0088] A rigid time shift mapping is set , , :

[0089] ;

[0090] The candidate phase difference is required for the non-standardized and , and the following exponential coding is performed:

[0091] ;

[0092] wherein is a predetermined normal number to control the upper limit of the power of the numerical value, which is artificially set, and the larger the value, the more obvious the phase feature; the processed and are single-channel data that highlight the phase feature.

[0093] Further, the step c23 specifically comprises the following steps:

[0094] The corresponding cumulative distribution function (CDF) and are calculated and , and the difference between the inverse functions of the CDF and is obtained to obtain the candidate time shift, and the candidate set is :

[0095] ;

[0096] ;

[0097] Then, among all the candidate , an optimal time shift is found, which makes the above target function value minimum, and the corresponding transmission plan :

[0098] ;

[0099] Step c24 specifically comprises the following steps:

[0100] Constraint condition Wherein, The maximum movable phase, in order to ensure the accuracy of the algorithm, here Take half of the source wavelet Length , record The optimized normalized data Then the following inverse normalization process can be used instead of Participate in the convolution correlation operation:

[0101] .

[0102] Compared with the prior art, the beneficial effects of the present application are:

[0103] 1、The present application first systematically introduces the optimal transmission theory into the field of Marchenko multiple wave suppression, and fundamentally solves the inherent limitation of the traditional point-by-point comparison method which cannot distinguish between local amplitude error and global waveform difference through the global constraint characteristics of Wasserstein distance;

[0104] 2、Through the division and optimization of odd and even iterations, the method is corrected for phase mismatch and amplitude distortion respectively, so that it has stronger robustness for noise, irregular sampling and amplitude attenuation in actual land data;

[0105] 3、Actual test shows that the multiple wave suppression precision of the present application is improved compared with the traditional method, and it has a unique advantage in overcoming the lossless medium assumption condition, and compared with the leading method, the improved Marchenko interlayer multiple wave suppression method based on compressed sensing for land actual data, the results of the present application can obviously eliminate the artifacts caused by amplitude mismatch caused by actual attenuation ring, and provide reliable technical support for high-precision seismic imaging under complex geological conditions;

[0106] 4、The present application successfully improves the effect of Marchenko method in land exploration, and can realize high-precision multiple wave suppression without velocity model, and significantly improves the reliability of seismic data processing and the accuracy of oil and gas exploration. BRIEF DESCRIPTION OF DRAWINGS

[0107] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed to be used in the embodiments will be briefly introduced as follows. It should be understood that the following drawings only show some of the embodiments of the present application, and therefore should not be regarded as a limitation on the scope. Other related drawings can also be obtained by those of ordinary skill in the art without creative labor on the basis of these drawings.

[0108] Figure 1 The step flow chart of the Marchenko multiple suppression method based on the wave field mismatch correction of the Wasserstein distance of the present application;

[0109] Figure 2 The velocity model corresponding to the model data;

[0110] Figure 3 The single shot processing effect comparison chart of the model data, wherein (a) is the original input single shot record, (b) is the single shot record after processing based on the multiple suppression method of compressive sense, (c) is the processing result of the present application, (d) is the local amplification of the original input single shot record, (e) is the local amplification of the single shot record after processing based on the multiple suppression method of compressive sense, and (f) is the local amplification of the processing result of the present application;

[0111] Figure 4 The zero offset gather before and after processing of the model data, wherein (a) is the zero offset gather of the original input data, (b) is the zero offset gather of the data after processing based on the multiple suppression method of compressive sense, and (c) is the zero offset gather of the data processed by the present application;

[0112] Figure 5 The imaging result before and after processing of the model data, wherein (a) is the original velocity model, (b) is the imaging result of the original input data, (c) is the imaging result of the data after processing based on the multiple suppression method of compressive sense, and (d) is the imaging result of the data processed by the present application;

[0113] Figure 6 The single shot processing effect comparison chart of the actual land data, wherein (a) is the original input single shot record, (b) is the single shot record after processing based on the multiple suppression method of compressive sense, and (c) is the processing result of the present application;

[0114] Figure 7 The zero offset gather before and after processing of the actual land data, wherein (a) is the zero offset gather of the original input data, (b) is the zero offset gather of the data after processing based on the multiple suppression method of compressive sense, and (c) is the zero offset gather of the data processed by the present application;

[0115] Figure 8The imaging results of the actual land data before and after processing, wherein (a) is the original input data imaging result, (b) is the data imaging result after processing based on the compression sense multiple wave suppression method, and (c) is the data imaging result after processing by the method. DETAILED DESCRIPTION

[0116] The application will be further described below in conjunction with the embodiments.

[0117] The application will be further described below in conjunction with the embodiments.

[0118] It should be noted that: similar labels 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 the subsequent drawings. Meanwhile, in the description of the application, the terms "first", "second" and the like are only used for distinguishing description, and cannot be understood as indicating or implying relative importance.

[0119] The application provides a multiple wave suppression method for correcting wave field mismatch based on Wasserstein distance of optimal transmission theory, which innovatively introduces Wasserstein distance into Marchenko iteration framework, and effectively solves the error propagation problem existing in the application of traditional methods in actual land seismic data.

[0120] The method comprises three links of data preprocessing, Marchenko iteration initialization and iteration of wave field mismatch correction of Wasserstein distance. In the data preprocessing link, the application adopts a complete data processing flow based on the compression sensing theory, including observation system regularization, Shearlet domain denoising, wave field reconstruction and sparse constraint deconvolution and the like, so as to provide high-quality input data for subsequent iteration and the possibility of comparison with the result of the front-line method, that is, the improved Marchenko interlayer multiple wave suppression method for land actual data based on compression sensing. In the core iteration link of wave field mismatch correction of Wasserstein distance, the application designs an optimization strategy of alternating even and odd iterations according to different characteristics of the error generation mechanism in the Marchenko iteration process: in the even iteration, the global phase correction of the predicted multiple wave field is realized by using the sensitivity of the Wasserstein distance to the phase change; in the odd iteration, the amplitude-waveform accurate matching is realized by constructing a composite Wasserstein distance objective function to simultaneously constrain the space-time structure and amplitude characteristics of the wave field.

[0121] The Marchenko multiple wave suppression method based on the wave field mismatch correction of Wasserstein distance comprises the following steps:

[0122] a: obtaining the impulse reflection sequence by preprocessing and deconvolution:

[0123] Since the Marchenko iteration needs to input the impulse reflection sequence processed by preprocessing and deconvolution, the original seismic data collected cannot be directly used as the input of the Marchenko multiple wave removal method. Therefore, the application adopts the existing method based on the compression sensing theory to systematically process the original land seismic data, so as to improve the data quality.

[0124] Specifically, the original collected seismic data is first reorganized into a three-dimensional data body, and a regular grid is established; then the Shearlet transform with multi-scale and multi-direction analysis capability is used to convert the data to a sparse domain; then the hard threshold denoising processing and wave field reconstruction are sequentially performed in the Shearlet domain, so as to suppress noise and repair missing channel data; finally, the high-resolution impulse reflection response is obtained through the sparse constraint deconvolution technology , so as to provide high-quality input data close to the basic assumption of the Marchenko method for the subsequent multiple wave suppression link. This process can be summarized as the following formula:

[0125] ;

[0126] Wherein, is the original observation data, is a sampling operator, indicating the data missing position, is the source wavelet convolution matrix, denotes the convolution process, is the unknown reflectivity to be solved, is the reflectivity sparsity, is the Shearlet transform, denotes the seismic record sparsity, and denote the L1 norm and L2 norm, respectively, and is the regularization parameter.

[0127] The core of the pre-processing procedure can be summarized as solving a unified optimization problem, which aims to obtain a high-quality impulse reflectivity . The optimization problem consists of three terms: 1. The data fidelity term ensures that the processed result fits the original observed data . 2. is the reflectivity sparsity, which is constrained by the L1 norm to force the solution to be sparse in the specific Shearlet domain, aiming to strip the influence of the source wavelet and thus improve the resolution of the reflectivity. 3. denotes the seismic record sparsity, which constrains the seismic record synthesized by the solution to be also sparse in the Shearlet domain, simultaneously achieving effective wavefield reconstruction and noise suppression. Through the adjustment of the regularization parameters and , the three terms work together to ensure that the final output has the characteristics of high fidelity, high resolution, and low noise, providing ideal input data for subsequent Marchenko multiple suppression.

[0128] b: input impulse reflectivity sequence, initialize Marchenko iteration:

[0129] Take the obtained in step a as the input to initialize the wavefield definition, set the initial primary wavefield , set the initial multiple prediction term , set the initial iteration number and the maximum iteration number , and the iteration process is looped under the condition of , and is divided into odd and even iterations. Among them, the odd iteration updates the primary wavefield: the original wavefield minus the predicted multiple wavefield; the even iteration updates the predicted multiple wavefield: the time window term acts on the last odd predicted multiple wavefield.

[0130] Specifically, ; ;

[0131] wherein, is a certain time window for cutting off the part of the data outside . . is the initial primary wave field, is the initial multiple wave field;

[0132] In the present application, the odd-numbered iteration mainly performs multiple wave subtraction, that is, subtracting the multiple wave field predicted in the last even-numbered iteration from the primary wave field in the last odd-numbered iteration, and the traditional Marchenko multiple wave removal method is described as follows:

[0133] ;

[0134] wherein, is the updated primary wave field in the current odd-numbered iteration, is the updated primary wave field in the last odd-numbered iteration, is the multiple wave field predicted in the last even-numbered iteration.

[0135] The even-numbered iteration (n is even) mainly performs multiple wave prediction:

[0136] ;

[0137] wherein, is the predicted multiple wave field in the current even-numbered iteration, is the predicted multiple wave field in the last even-numbered iteration, and respectively represent the cross-correlation and convolution forms of , and for the convenience of representation, is denoted as .

[0138] The conventional Marchenko iteration is to perform the above-mentioned odd-even iteration cycle, and when , the iteration stops, and all time slices of updated in the last iteration are collected as the final primary wave field.

[0139] c: Marchenko iteration with remolding correction of wave field mismatch based on Wasserstein distance:

[0140] c1: Marchenko odd-numbered iteration with remolding correction of predicted wave field amplitude and waveform based on Wasserstein distance:

[0141] ​Observing the regular odd-numbered iterations of Marchenko, imperfections in the actual data can lead to mismatches in the two terms being subtracted, resulting in multiple predicted wave fields appearing during the subtraction process. It will still be related to the original wave field The multiple wave components in the wavefield may exhibit mismatches in waveform and amplitude. Such subtraction can damage the primary wave or cause incomplete attenuation of the multiple waves. Therefore, this invention should match the amplitude and waveform characteristics of the predicted wavefield and the original wavefield in the odd-numbered reshaping iteration system. To maintain inter-channel continuity, global optimization using two-dimensional spatiotemporal data is required. Since the predicted multiple wavefield contains both multiple and primary wave components, direct full-wavefield matching would damage the effective signal. Therefore, this invention designs a composite transmission mapping cost: amplitude-waveform mapping is only performed on regions in the predicted wavefield that are highly similar in shape to the observed data, while significantly different regions are identified as primary waves and retained. This process is achieved by constructing a Wasserstein distance objective function weighted by spatiotemporal cost and amplitude cost, where the spatiotemporal term constrains the wavefront geometry and the amplitude term controls energy distribution, thereby ensuring the integrity of the primary wave and the spatial continuity of the wavefield while suppressing multiple waves. Finally, an optimal mapping is obtained to transform the originally erroneous multiple wavefield into an optimal wavefield.

[0142] c11: Standardize the input data of the original wavefield and the predicted multiple wavefield. First, perform non-negative standardization on the original seismic trace to eliminate negative amplitude values; then, perform amplitude normalization to make the data meet the basic mathematical requirements of distribution non-negativity and total mass conservation in the optimal transmission theory.

[0143] Wasserstein distance calculations require non-negative data and a uniform total measure. However, in seismic data, due to positive and negative fluctuations in signal amplitude and uncertain total energy, data standardization is necessary before directly applying Wasserstein distance calculations. Preprocessing ensures the validity of the Wasserstein distance metric and the rationality of its physical meaning. The two-dimensional data used in the calculations are standardized, and... for , for ,in, These represent the number of time sampling points and the number of spatial channels, respectively. The standardization process, expressed in discrete point form, is as follows:

[0144] ;

[0145] ;

[0146] For standardization , For standardization , denotes the minimum value in , denotes the maximum value in , denotes the minimum value in , denotes the maximum value in , is a very small positive number to prevent division by zero.

[0147] c12: Set up the corresponding objective function to find an optimal transmission mapping that minimizes the total transportation cost:

[0148] ;

[0149] where, is and the total transportation cost function, is and the transportation cost function of a single point .

[0150] c13: Extract the space-time-amplitude index, and design the space-time-amplitude composite cost function:

[0151] Let and be the time index, trace (space) index, and amplitude value of and :

[0152] ;

[0153] where, and are the time sampling point indices of and , and are the space trace number indices of and , and are the amplitude values corresponding to the time and space indices of and , is the degree of space-time index, and the higher the degree means the greater the gradient of the subsequent space-time cost function, so the amplitude correction will be limited to a very small range.

[0154] Determine the space-time range control term in the composite cost function as:

[0155] ;

[0156] The amplitude control term in the composite cost function is determined as:

[0157] ;

[0158] The composite cost function is determined as:

[0159] ;

[0160] and are two weight constants for the spatio-temporal cost and the amplitude cost, which are artificially set.

[0161] c14: Sinkhorn iteration to solve the optimal transport map:

[0162] Since there are multiple point distributions in two-dimensional data, for point distribution, the number of variables of the linear programming is , which means that a linear programming problem containing millions of variables needs to be solved, which is completely infeasible in calculation. And the solution of linear programming is usually sparse, which is not conducive to guiding the smooth transition of data. In order to overcome the above bottleneck, the present application introduces an entropy penalty term in the objective function, and the problem is transformed into the following form and then solved efficiently by Sinkhorn iteration:

[0163] ;

[0164] wherein, is the entropy of the transport plan , and is a regularization coefficient. The sum of each row element of the transport map matrix must be equal to the corresponding component of the source distribution un, and the sum of each column element must be equal to the corresponding component of the target distribution mn, wherein represents a column vector with all elements being 1, which ensures mass conservation.

[0165] The Sinkhorn iteration is used to solve the solution of the problem, as follows:

[0166] ;

[0167] ;

[0168] wherein, is an intermediate construction matrix for calculating the optimal transport map , which is completely determined by the cost function The decision is that the larger the cost function, the smaller the value of K. In the formula, and This is an optional scaling vector, its purpose is to ensure... To satisfy the marginal constraint of mass conservation, Representing a diagonal matrix, express .

[0169] The iterative process is as follows:

[0170] First, initialize Next, the alternating iterations are performed for... :

[0171] ;

[0172] ;

[0173] in, This represents element-wise division. and This represents the scaling vector updated in the current iteration.

[0174] Iterate until and The change is less than a certain threshold. The final solution for the optimal transport mapping is:

[0175] ;

[0176] For this optimal transport mapping In this regard, the optimization for predicting multiple wave fields is as follows:

[0177] ;

[0178] in, A constant is manually set to control the update rate. Indicates the direction of transport, i.e., whether the amplitude decreases or increases; the subscript here... and express and The element index in the array. For the optimally mapped The elements in the corrected standardized multiple wave field, and the overall updated multiple wave field, need to be inversely standardized to restore it to a state where multiple wave subtraction can be performed.

[0179] c15: Updates were subtracted from Tron multiple times:

[0180] For each of the collected of Multiple wave fields after overall elemental correction Perform inverse standardization:

[0181] ;

[0182] This refers to the multiple wave field subtracted from the final corrected multiple wave field.

[0183] c2: The Marchenko even-numbered iterations of the reconstructed wavefield phase based on Wsserstein distance correction:

[0184] Observing the regular Marchenko even-numbered iterations, due to the complex convolution and correlation processes, it is easy to cause... (Right now )and The phase mismatch between them, in The process involves multiple time windows, specifically a strictly causal time window. The core mathematical principle of this algorithm is to construct multiple wave paths based on the temporal (phase) relationship of events. Inaccurate phase relationships can lead to errors in the window, resulting in the removal of more or less data. Accurate phase relationships are crucial for generating physically reliable data. This is a prerequisite. This phase mismatch will produce errors in phase or amplitude in the calculated values. These amplitude errors can be adjusted in subsequent odd-numbered iterations. In this step, correcting the phase mismatch is considered the core issue. Therefore, this invention proceeds in this step. Phase correction.

[0185] c21: Standardize the input data for the time window term and the previous prediction of multiple wave fields, performing non-negative standardization on a channel-by-channel basis:

[0186] In even-numbered iterations, the phase-matching correction process is performed track-by-track, and the selected object is... and The continuous function forms of the single-channel records in the corresponding time slices are as follows: and Similarly, because the signal amplitude in seismic data fluctuates (positive and negative) and the total energy is uncertain, direct application of Wasserstein distance calculation requires prior non-negative data standardization.

[0187] ;

[0188] ;

[0189] in, To prevent extremely small positive numbers from being divided by zero, for The minimum value in, for The minimum value in, is the maximum value in , the and are the corresponding two normalized traces and are both non-negative and have the same measure, the original traces are normalized to be after non-negative normalization.

[0190] c22: Exponential coding enhanced phase feature:

[0191] After normalization, the Wassertein distance objective function of even iteration amplitude correction is defined as:

[0192] ;

[0193] An optimal transport mapping is expected to satisfy the above objective function, and a rigid time shift mapping is set, which is used to shift on the time axis by :

[0194] ;

[0195] Obviously, to determine the optimal mapping, first, the candidate phase difference is required. To enhance the phase feature of the seismic trace, the present application adopts an exponential coding method. For the non-normalized and , the following exponential coding is performed:

[0196] ;

[0197] wherein is a specified normal number to control the upper bound of the power of the numerical precision, which is artificially set, and the larger the value, the more obvious the phase feature. Since the fluctuation of the exponential function on the negative half real axis is weak, the above procedure distinguishes the positive and negative parts of the seismic data to ensure that the fluctuation degree of the positive and negative parts of the function is increased, and the processed and are single-channel data that highlight the phase feature.

[0198] c23: Calculate the CDF inverse function of the normalized time window item and the multiple wave field to obtain the candidate time shift:

[0199] Calculate the cumulative distribution functions (CDFs) of and corresponding to and , and obtain the CDF inverse functions and Difference is obtained to obtain candidate time shift Set is :

[0200] ;

[0201] ;

[0202] Then, find an optimal time shift in all candidate The objective function value is the smallest The corresponding transmission plan :

[0203] .

[0204] c24: optimal time shift mapping that minimizes the Wasserstein cost function, update the predicted multiple wave field.

[0205] There are constraints Where, Is the maximum movable phase, in order to ensure the accuracy of the algorithm here Take half the length of the source wavelet , , record The normalized data after optimization , then the following inverse normalization process can replace Participate in convolution correlation operation:

[0206] ;

[0207] c3: odd and even cycle iterations can complete Marchenko iteration based on Wasserstein distance correction wave field mismatch, finally output the primary wave field.

[0208] The innovation of the application is that the wave field is updated by remodeling Marchenko odd iteration based on Wsserstein distance measurement and remodeling Marchenko even iteration based on Wsserstein distance measurement, so as to achieve more stable and less error interlayer multiple wave elimination results. The Wsserstein distance optimization strategy is to find a mapping strategy that minimizes the mapping cost and minimizes the Wasserstein distance for two different distributions, so that one distribution is mapped to a zero-one distribution through the mapping.

[0209] Compared with existing methods, Wasserstein distance provides a fundamental theoretical breakthrough. Although the traditional L2 norm constraint can suppress noise, its global average effect will cause wavefront blurring and loss of key kinematic features; L1 or TV norm can maintain edges, but it will cause systematic distortion of amplitude information. These point-by-point comparison-based methods cannot effectively distinguish between "local amplitude error" and "global waveform difference". In contrast, Wasserstein distance measures the overall shape similarity between probability distributions to achieve distribution matching with the minimum transmission cost, which not only sensitively captures kinematic features such as phase shift and time delay, but also completely maintains the relative relationship of the amplitude in the matching process. This feature makes it particularly suitable for the matching problem of seismic wave field with strong physical constraints.

[0210] Compared with the current most advanced land-based measured data Marchenko multiple wave suppression method based on compressed sensing, the method realizes the paradigm shift from "data preprocessing" to "inversion mechanism". CS-MME mainly improves the quality of input data through sparse constraints such as Shearlet transform, but its iterative kernel still follows the traditional point-by-point comparison method, which cannot fundamentally solve the error propagation problem. The purpose of the present invention is to innovate directly from the level of inversion mechanism, and a unique two-way constraint framework is constructed. In even-numbered iterations, the present invention introduces Wasserstein distance as the phase correction criterion, and realizes the global calibration of the phase mismatch between the predicted multiple wave field and the observed data by optimizing the transmission plan, ensuring the causality and physical reasonableness of wave field propagation; in odd-numbered iterations, a composite Wasserstein objective function is used to simultaneously constrain the amplitude characteristics and spatio-temporal structure of the wave field, and by introducing a transmission mechanism weighted by the spatio-temporal cost and amplitude cost, the dynamic information of the primary wave is effectively protected while suppressing multiple waves.

[0211] This constraint mechanism based on overall distribution matching makes the method inherently robust to data loss and noise. The present method breaks through the limitations of traditional constraint methods, while maintaining the advantage of Marchenko method without velocity model, significantly improving its adaptability to noise, loss and amplitude distortion of actual data, thereby realizing higher precision multiple wave suppression and more reliable primary wave recovery, and providing effective technical support for high-precision seismic imaging in complex land exploration environment. The present method shows significant advantages in constraint mechanism, error control and global feature protection.

[0212] Embodiment 1:

[0213] This embodiment has been applied to model data and a land-based measured data and the effect is analyzed.

[0214] The model data is from Figure 2The velocity model shown is obtained by forward modeling. The model is based on a partial version of the Pluto model published by SMAART. The Pluto model is a standardized typical benchmark for testing the effectiveness of multiple wave suppression algorithms in the industry. Its partial version contains salt dome structures with significant velocity differences and complex internal structures, with a velocity range of 1500 meters / second to 4500 meters / second. The existence of high-speed salt bodies and underlying high-speed reflection layers forms a strong wave impedance contrast with the surrounding strata. This large velocity jump and multi-layer characteristics are conducive to generating complex and significant interlayer multiple waves between the two high-speed layers, providing a good basis for subsequent testing. Based on the model, a finite difference method is used to forward simulate a data set rich in interlayer multiple waves. The data set contains 201 shots and 201 channels per shot. The shot points and receivers are evenly distributed with a spacing of 25 meters. The receiver records have a time sampling interval of 4 milliseconds, containing 1751 time samples.

[0215] The actual land data acquisition parameters are a shot point spacing of 200 meters and a receiver spacing of 25 meters. After regularization processing, an interpolated data volume is formed with shot points and receivers staggered at intervals of 12.5 meters. The final data set contains 331 shots and 331 channels per shot. The total recording time is 6 seconds, with an 8 millisecond sampling interval containing 750 time samples. The effective signal is mainly concentrated in the first 4 seconds, while the interlayer multiple waves are mainly developed after 2.5 seconds.

[0216] As shown in Figure 1 , a Marchenko multiple wave suppression method based on Wasserstein distance correction of wave field mismatch includes the following steps:

[0217] a: Marchenko iteration needs to input preprocessed and deconvolved impulse reflection sequences. Raw seismic data obtained by original acquisition cannot be directly used as input for Marchenko multiple wave removal methods. We use existing front-line methods based on the technical process of compressive sensing theory to systematically process raw land seismic data to improve data quality while ensuring the comparability of the present invention with existing advanced methods, proving that the present method has further improved effect compared to the front-line method. Preprocessing first reorganizes the original data into a three-dimensional data volume and establishes a regularized grid through observation system regularization. Then, the Shearlet transform with multi-scale and multi-directional analysis capability is used to convert the data to the sparse domain. In the Shearlet domain, hard threshold denoising and wave field reconstruction are performed in sequence to effectively suppress noise and repair missing channel data. Finally, sparse constraint deconvolution technology is used to obtain high-resolution impulse responses, providing high-quality input data close to the basic assumptions of the Marchenko method for subsequent multiple wave suppression steps. This process can be summarized as follows:

[0218] ;

[0219] The core of this preprocessing procedure can be summarized as solving a unified optimization problem, whose goal is to obtain high-quality pulse reflection responses . This optimization problem consists of three key components: a data fidelity term ensuring that the processed result faithfully fits the original observed data , where is the sampling operator, indicating the position of data missing, is the source wavelet convolution matrix, denoting the convolution process, is the unknown reflection response to be solved. is a reflection coefficient sparsity term, and denote the L1 norm and L2 norm, respectively, is the Shearlet transform, which is constrained by the L1 norm to make the solution itself sparse in the specific Shearlet domain, aiming to strip the influence of the source wavelet and thus improve the resolution of the reflection coefficient; and denotes a seismic record sparsity term, which constrains the seismic record synthesized by the solution to also be sparse in the Shearlet domain, which can simultaneously achieve effective wavefield reconstruction and noise suppression. Through the adjustment of the regularization parameters and , these three terms work together to ensure that the final output has the characteristics of high fidelity, high resolution, and low noise, providing ideal input data for subsequent Marchenko multiple suppression. In the model data of the embodiment of the present application, and are 0.3 and 0.5, and in actual land data, and are 0.4 and 0.4.

[0220] b: Perform Marchenko iteration initialization after preprocessing:

[0221] First, the basic formula for initializing the conventional Marchenko iteration is given: Marchenko iteration first initializes the wavefield definition according to the only input obtained in the first step , and sets the initial primary wavefield:

[0222] ;

[0223] Here, is the initial primary wavefield.

[0224] Set initial multiple prediction term:

[0225]

[0226] where, is the initial set multiple prediction term, is a specific time window, and data outside the time window is cut off. Set initial iteration number

[0227] and maximum iteration number , and the iteration process is looped under the condition that and is divided into odd and even iterations: Odd iteration mainly performs multiple subtraction, that is, subtracts the multiple field predicted by the last even iteration from the primary field of the last odd iteration. The traditional Marchenko multiple removal method describes odd iteration as (this will be reshaped as innovative Marchenko iteration based on Wasserstein distance in S3):

[0228]

[0229]

[0230] where, is the updated primary field in the current odd iteration, is the updated primary field in the last odd iteration, is the multiple field predicted by the last even iteration.

[0231] Even iteration (n is even) mainly performs multiple prediction:

[0232]

[0233] where, is the predicted multiple field in the current even iteration, is the predicted multiple field in the last even iteration, and respectively represent the cross-correlation and convolution forms of , and for convenience of representation,

[0234] The conventional Marchenko iteration is to perform the above odd-even iteration loop, and when the iteration stops, and the of all time slices updated in the last iteration is collected as the final primary field.

[0235] ​​​​​​​c: Marchenko iteration loop with remodeled Wsserstein distance correction of wavefield mismatch:

[0236] The innovation of the present application is that the interbed multiple elimination result with more robust and less error is achieved by the S3 remodeled Marchenko odd iteration based on Wsserstein distance measurement and S4 remodeled Marchenko even iteration based on Wsserstein distance measurement, and the basic idea of Wsserstein distance optimization strategy is to find a mapping strategy with the lowest mapping cost and the minimum Wasserstein distance for two different distributions, so that one distribution is mapped to a zero-one distribution through this mapping, and specifically:

[0237] c1: Marchenko odd iteration with remodeled Wsserstein distance correction of predicted wavefield amplitude and waveform: observing the conventional Marchenko odd iteration, the actual data imperfection factors described in the background will cause mismatch factors in the two items to be subtracted, resulting in predicted multiple wavefield in the subtraction process Still exist waveform and amplitude mismatch with multiple wave components in the original wavefield Such subtraction will damage the primary wave in the wavefield or make the multiple wave attenuation incomplete. Therefore, the amplitude and waveform characteristics of the predicted wavefield and the original wavefield should be matched in the odd iteration system of the remodeled iteration of the present application. In this link, in order to maintain inter-channel continuity, two-dimensional space-time data is used for global optimization. Since the predicted multiple wavefield contains multiple wave and primary wave components, direct full wavefield matching will damage the effective signal. Therefore, the present application designs a composite transmission mapping cost: only the amplitude-waveform mapping is performed on the area in the predicted wavefield that is highly similar to the observed data form, and the significantly different area is determined as a primary wave to be reserved. This process is realized by constructing a Wsserstein distance objective function weighted by space-time cost and amplitude cost, wherein the space-time item constrains the wavefront geometry, and the amplitude item controls the energy distribution, so as to suppress multiple waves while ensuring the integrity of primary waves and the spatial continuity of wavefield. Finally, an optimal mapping is obtained to transform the originally error-prone multiple wavefield into an optimal wavefield.

[0238] c11: Wasserstein distance calculation requires data to be non-negative and have the same total measure. In seismic data, due to the positive and negative fluctuations of signal amplitudes and the uncertainty of total energy sum, direct application of Wasserstein distance calculation requires data standardization. The specific process includes: first, non-negative processing of the original seismic trace to eliminate negative amplitudes; then, amplitude normalization is implemented to make the data meet the basic mathematical requirements of the optimal transport theory about the non-negativity of distribution and the conservation of total mass. The preprocessing ensures the effectiveness of the Wasserstein distance measurement and the rationality of the physical meaning. The two-dimensional data participating in the operation are standardized, denoted as respectively, the number of time sampling points and the number of spatial traces. The standardization process is written in the form of discrete points as follows:

[0239]

[0240]

[0241]

[0242] c12: After standardizing the input data, the corresponding objective function needs to be set. Similar to the classical optimal transport problem, the objective function is still to find an optimal transport mapping to minimize the total cost of transportation:

[0243]

[0244]

[0245] ​​​​​​​​​​​​​​​​​​​​​​​​​​​​​c13: in order to design a suitable , we will unfolding the time index, trace (space) index and amplitude value of into a three-dimensional point set and :

[0246] ;

[0247] wherein, and are the time sampling point indexes of and respectively, and are the space (trace number) indexes of and respectively, and are the amplitude values corresponding to the time and space indexes of and respectively, is the degree of spatiotemporal index, the higher the degree means the greater the gradient of the subsequent spatiotemporal cost function, so that the amplitude correction will be limited in a very small range, in the model data of the present application is 3, and in the actual land data is 2.

[0248] The spatiotemporal range control term in the composite cost function is determined as:

[0249] ;

[0250] The amplitude control term in the composite cost function is determined as:

[0251] ;

[0252] The composite cost function is determined as:

[0253] ;

[0254] and are two weight constants for controlling the spatiotemporal cost and the amplitude cost, which are artificially set, and in the model data of the present application and are 3 and 5 respectively.

[0255] c14: solving the transmission mapping , the purpose is to obtain the optimal transmission mapping Since there are multiple points in the two-dimensional data, for the distribution of points, the number of variables of the linear programming is , which means that a linear programming problem containing millions of variables needs to be solved, which is completely infeasible in calculation. And the solution of linear programming is usually sparse, which is not conducive to guiding the smooth transition of data. In order to overcome the above bottleneck, we introduce an entropy penalty term in the objective function, and convert the problem into the following form and solve it efficiently using Sinkhorn iteration:

[0256] ;

[0257] Where, is the entropy of the transmission plan , and is the regularization coefficient. Control the sum of each row element of the transmission mapping matrix T_S, which must be equal to the corresponding component of the source distribution un, and the sum of each column element, which must be equal to the corresponding component of the target distribution mn, among which represents a column vector with all elements being 1, which ensures mass conservation.

[0258] The solution of this problem solved by Sinkhorn iteration can be written in the following form:

[0259] ;

[0260] ;

[0261] Where, is the intermediate construction matrix for calculating the optimal transmission mapping , which is completely determined by the cost function . Here is a preliminary explanation of a rule, that is, when the cost function is larger, the value corresponding to K is smaller. and are optional scaling vectors, which ensure that satisfy the marginal constraints of mass conservation, represents a diagonal matrix form, represents .

[0262] The iteration process is:

[0263] First, initialize . Then expand the alternating iteration for :

[0264] ;

[0265] ;

[0266] where, denotes element-wise division, and denotes the scaling vector updated in the current iteration.

[0267] iterating until and change less than a certain threshold. The final solution of the optimal transport mapping is:

[0268] ;

[0269] For this optimal transport mapping , the optimization of the predicted multiple wavefield is:

[0270] ;

[0271] where, is an artificial constant to control the updating degree, and the model data and the actual land data in the present application , indicates the carrying direction, i.e. the amplitude decreases or increases, where the subscript and denote and . is the element in the normalized multiple wavefield after the correction of the optimal mapping , and the overall updated multiple wavefield needs to be inversely normalized to restore to a state where multiple wave subtraction can be performed.

[0272] c15: inversely normalizing the overall corrected multiple wavefield of each element collected :

[0273] ;

[0274] is the multiple wavefield after the final correction which participates in multiple wave subtraction.

[0275] c2: reshaping the Marchenko even-time iteration based on the Wsserstein distance corrected predicted wavefield phase: observing the conventional Marchenko even-time iteration, due to the complex convolution and correlation process, it is easy to cause (i.e. ) and phase mismatch between The item contains many time window participations, the window is a strict causal time window, the mathematical kernel of this algorithm is based on the time (phase) relationship of events to construct multiple wave paths, the error caused by inaccurate phase will cause the window to cut more or less data, accurate phase relationship is the premise of generating physically reliable . This phase mismatch will cause errors in the calculated phase or amplitude, which can be adjusted in subsequent odd iterations, in this step we consider the correction of phase mismatch as the core problem, so in this step the invention carries out Phase correction; Specifically:

[0276] c21: In even iterations, the phase matching correction process is carried out channel by channel, and the selected object is And Its continuous function form is respectively And Similarly, due to the positive and negative fluctuations of signal amplitude and the uncertainty of total energy in seismic data, direct application of Wasserstein distance calculation needs to be standardized first:

[0277] ;

[0278] ;

[0279] Where, is a small positive number to prevent division by zero, is the minimum value in , is the minimum value in , is the maximum value in , and the standardized and corresponding to the two standardized channels and are non-negative and have the same measure. After non-negative normalization, the original channel will be unified to Standard.

[0280] c22: After such standardization, the Wasserstein distance objective function of even iteration amplitude correction is defined as:

[0281] ;

[0282] We hope to find an optimal transport mapping that satisfies the above objective function, and set a rigid time shift mapping , which is used to shift on the time axis :

[0283] ;

[0284] It is obvious that to determine the optimal mapping, first requires the candidate phase difference In order to enhance the phase characteristics of seismic traces, we use an exponential encoding method, for the non-standardized And , the following exponential encoding is performed:

[0285] ;

[0286] Wherein, Is a predetermined normal number to control the upper bound of the power of the numerical precision, artificially set, the model data of the present application and the actual land data The larger the value, the more obvious the phase characteristics. Since the exponential function fluctuates weakly on the negative half real axis, the above program treats the positive and negative parts of the seismic data differently to ensure that the fluctuation degree of the positive and negative parts of the function has been increased, and the processed And Is a single-channel data that highlights the phase characteristics.

[0287] c23: Calculate And The corresponding cumulative distribution function (CDF) And , the candidate time shift is obtained by the difference between the inverse function of CDF And Get the candidate Set :

[0288] ;

[0289] ;

[0290] Then we find an optimal time shift in all candidate The above objective function value is minimum Corresponding to the transmission plan :

[0291] ;

[0292] c24: There is a constraint condition , wherein Is the maximum movable phase, in order to ensure the accuracy of the algorithm, here Take the length of half of the source wavelet , , record for the optimized standardized data then the following inverse standardization process can be used to replace the data Participate in the convolution correlation operation:

[0293] ;

[0294] Marchenko iteration based on Wasserstein distance correction of wave field mismatch can be completed by iterating between step c1 and step c2.

[0295] The application effect of the present application can be divided into two groups, the first group is model data, and the second group is actual land data:

[0296] For model data:

[0297] As shown in Figure 3 , the original data and the results of the front method and the method proposed in the present application are compared. By observing Figure 3 (a), Figure 3 (b), Figure 3 (c) complete shot gather data, it is found that interlayer multiple waves mainly develop after 3 seconds. In order to facilitate display, the area after 3 seconds is enlarged to Figure 3 (d), Figure 3 (e), Figure 3 (f), and the results of processing the original data by the two systems are compared and analyzed. Figure 3 (d) can be seen, the front method result (e) can achieve a certain degree of multiple wave suppression, but compared with the area indicated by the arrow Figure 3 (f) shows that the method proposed in the present application achieves better suppression effect, and it can be observed that the front method result produces false images due to obvious over subtraction, while the same position in Figure 4 (f) is very clean.

[0298] To verify the effectiveness of the method proposed in the present application on full shot gather data, the zero offset gather of the original data is calculated, as shown in Figure 4 (a), the zero offset gather of the data processed by the front method is shown in Figure 4 (b), and the zero offset gather of the data processed by the method proposed in the present application is shown in Figure 4 (c), for clearer display effect, the data is locally enlarged, and the position indicated by the arrow in Figure 4 (a), Figure 4 (b) and Figure 3 (c) can be seen, the proposed method has better full shot gather suppression effect, and the position of multiple wave removal coincides with the event of eliminating multiple wave events in the single shot record shown in Figure 5 .

[0299] The final purpose of verifying the method proposed in the present application is to eliminate the false events caused by multiple waves. The original data and the results of the method proposed in the present application are imaged by RTM. Figure 5 (a) shows the real model event information, in Figure 5 (b) shows the imaging results of the original data, in which many false events inconsistent with the real events can be seen. Comparing Figure 5 (c) and Figure 6 (d), it can be seen that in the results of the front method, there are false events or false events in the positions indicated by the black solid arrows in the enlarged area, while in the imaging results of the method proposed in the present application, there are no false events in the same position. The event information here is consistent with the model, and the comparison is very obvious. At the same time, it is found that in the position of the primary wave event indicated by the black hollow arrow, the result of the front method is discontinuous, while in the result of the method proposed in the present application, the continuity of the event is obviously improved, which proves that the proposed method has the ability to restore the primary wave.

[0300] For actual land data:

[0301] As Figure 6 shown, the single-shot results of the field original data, the front method results and the method proposed in the present application are compared. In order to compare more clearly, we enlarged the area where multiple waves develop; it can be seen that Figure 6 (a) has two obvious multiple wave development times indicated by arrows, in the processing results of the front method Figure 6 (b), the position indicated by the arrow is suppressed to a certain extent, but the residual of this multiple wave time is also obviously visible, while in the processing results of the method proposed in the present application Figure 7 (c), the two multiple wave events are suppressed very thoroughly.

[0302] The method proposed in the present application shows excellent multiple wave suppression effect in full shot data processing. By comparing the field original data, the front method results and the method proposed in the present application, as Figure 6 shown, it can be seen that the removal position of the partial enlarged view is highly consistent with Figure 8 the multiple wave elimination events in the single-shot record, effectively verifying the processing consistency of the method under different data scales.

[0303] To verify the ultimate goal of the proposed system - to eliminate the false events caused by multiple waves, the present application will stack the pre-stack time migration imaging results of the field raw data, the front-line method results and the data results after processing by the method proposed by the present application. For the convenience of display, we will magnify the data locally, and it can be observed that compared with the processing results of the front-line method, the method proposed by the present application exhibits more excellent ability to suppress false events in the position indicated by the black solid arrow. At the same time, more importantly, it is found that compared with Figure 8 (a) The field raw data shown are different from Figure 8 (b) The place indicated by the green arrow in the front-line method processing results not only fails to achieve the effect of event elimination, but also appears some false events with stronger amplitude. The reason is that the Marchenko method follows the assumption of acoustic wave without attenuation, while the actual data has attenuation effect, and the multiple wave field amplitude calculated by the front-line method is inaccurate non-attenuation wave field, so there is obvious over subtraction in the position of the green arrow, while the imaging results of the method proposed by the present application ​ (c) There is no false event in the same position, and all the multiple waves are correctly calculated and subtracted by the method proposed by the present application. It is proved that the method proposed by the present application can well relax the restriction of the acoustic wave without attenuation assumption of the Marchenko method.

[0304] Note that the design of the above numerical model, the parameter setting of each material, and the selection of the measurement plane are only the preferred embodiments of the present application and the technical principles applied. Those skilled in the art will understand that the present application is not limited to the specific embodiments described herein, and various obvious changes, readjustments and substitutions can be made by those skilled in the art without departing from the scope of the present application. Therefore, although the present application has been described in more detail through the above embodiments, the present application is not limited to the above embodiments, and can include more other equivalent embodiments without departing from the concept of the present application, and the scope of the present application is determined by the scope of the appended claims.

Claims

1. A Marchenko multiple suppression method based on Wasserstein distance correction of wave field mismatch, characterized in that, The method comprises the following steps: a: obtaining a pulse reflection sequence by pre-processing and deconvolution: Firstly, the original seismic data is reorganized into a three-dimensional data volume and a regular grid is established; secondly, the data is converted into a sparse domain by using a Shearlet transform with multi-scale and multi-direction analysis capability; thirdly, hard threshold denoising and wave field reconstruction are sequentially performed in the Shearlet domain to suppress noise and repair missing trace data; finally, high-resolution impulse reflection response is obtained by using sparse constraint deconvolution technology ​ b: inputting the pulse reflection sequence, and performing Marchenko iteration initialization: The result obtained in step a Initialize the wave field definition for the input, setting the initial primary wave field. Then set the initial multiple wave prediction term. Set the initial number of iterations. and maximum number of iterations The iterative process is in The process is repeated under certain conditions, and is divided into odd-numbered and even-numbered iterations. In the odd-numbered iterations, the wave field is updated once: the original wave field is reduced by the predicted multiple wave field. In the even-numbered iterations, the predicted multiple wave field is updated: the time window term is applied to the previous even-numbered predicted multiple wave field. c: Marchenko iteration of remolding and correcting wave field mismatch based on Wasserstein distance: c1: Marchenko odd-numbered iteration of remolding and correcting prediction wave field amplitude and waveform based on Wasserstein distance: c11: standardizing input data of the original wave field and the predicted multiple wave field, performing non-negative standardization on the original seismic trace to eliminate negative amplitude, and then performing amplitude normalization; c12: planning to set a corresponding objective function to find an optimal transport mapping to minimize the total transport cost: c13: extracting a space-time-amplitude index, and designing a space-time-amplitude composite cost function: c14: solving the optimal transport mapping by Sinkhorn iteration: c15: updating the subtracted multiple wave field; c2: Marchenko even-numbered iteration of remolding and correcting prediction wave field phase based on Wasserstein distance: c21: standardizing input data of the time window item and the last predicted multiple wave field, and performing non-negative standardization on each trace; c22: performing exponential coding to enhance phase characteristics; c23: calculating a CDF inverse function of the standardized time window item and the multiple wave field to obtain a candidate time shift; c24: solving an optimal time shift mapping that minimizes the Wasserstein cost function, and updating the predicted multiple wave field; c3: completing the Marchenko iteration of correcting wave field mismatch based on Wasserstein distance by odd-even cycle iteration, and finally outputting a primary wave field.

2. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step a, high resolution pulse reflectance response Obtained from the following formula: ； where, is the original observed data, is the sampling operator, which represents the data missing position, is the source wavelet convolution matrix, represents the convolution process, is the unknown reflectivity to be solved, is the reflectivity sparsity, is the Shearlet transform, represents the seismic record sparsity, and represent the L1 norm and L2 norm, respectively, and is the regularization parameter; data fidelity term , to ensure that the processed result fits the original observed data ; , is a reflectivity sparsity term, which is enforced by an L1 norm constraint to promote a solution that is inherently sparse in the specific Shearlet domain; , is a seismic record sparsity term, which constrains the synthesized seismic record to be also sparse in the Shearlet domain , and simultaneously enables efficient wavefield reconstruction and noise suppression.

3. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that: Step b, initial primary and multiple wavefield definitions are performed: ; ; wherein, is a certain time window for cutting off data outside of ; is an initial primary wavefield, is an initial multiple wavefield; The odd-numbered iteration mainly performs multiple wave subtraction, that is, subtracting the multiple wave field predicted in the last even-numbered iteration from the primary wave field in the last odd-numbered iteration, and the traditional Marchenko multiple wave removal method is described as follows: ； wherein, is the one-way wavefield updated in the current odd iteration, is the one-way wavefield updated in the previous odd iteration, is the multiple wavefield predicted in the previous even iteration; even-numbered iterations, For even-numbered iterations, mainly the multiples are predicted: ； wherein, is the predicted multiple wavefield for the current even iteration, is the predicted multiple wavefield for the previous even iteration, and respectively represent the cross-correlation and convolution form of is noted as ; The conventional Marchenko iteration is just to perform the above-mentioned odd-even iteration cycle, when the iteration stops, and all time slices of the last iteration update are collected as the final one-way wavefield.

4. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c11, specifically: standardization of the two-dimensional data involved in the operation, denoted as For , For Wherein, The number of time sampling points and the number of spatial channels, respectively. The standardization process is written in the form of discrete points as follows: ； ； is the standardized , is the standardized , denotes the minimum value in , denotes the maximum value in , denotes the minimum value in , denotes the maximum value in , is a very small positive number to prevent division by zero.

5. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c12, the total transport cost function is expressed as follows: ； wherein is and the total transmission cost function, is and the transmission cost function for a single point. the transmission cost function for a single point.

6. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c13, specifically comprising the following steps: :​​​​ ； wherein, is the time sampling point index, is the spatial channel index, is the amplitude size corresponding to the time and spatial index pair, is the time sampling point index, is the spatial channel index, is the amplitude size corresponding to the time and spatial index pair, is the time sampling point index, is the spatial channel index, is the amplitude size corresponding to the time and spatial index pair, is the time and spatial index pair, the higher the degree means the larger the gradient of the subsequent time and spatial cost function, then the amplitude correction will be limited in a very small range.​​​ Controlling spatiotemporal extent in composite cost functions is determined as: ； controlling the amplitude term in the composite cost function it is determined that: ； The composite cost function is defined as: f = f1 + f2 ； and are two control weight constants for the spatio-temporal cost and the amplitude cost, respectively, and are set artificially.

7. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c14, specifically comprising the following steps: An entropy penalty term is introduced into the objective function, and the problem is converted into the following form and then solved by Sinkhorn iteration: ； wherein is the entropy of the transmission plan , is a regularization coefficient; controls the sum of the elements of each row of the transmission mapping matrix , which must equal the corresponding component of the source distribution un, and the sum of the elements of each column, which must equal the corresponding component of the target distribution mn, where denotes a column vector of ones. The solution of the problem is solved by Sinkhorn iteration, as follows: ； ； wherein, To compute the optimal transmission map of the intermediate construction matrix, is completely determined by the cost function , the greater the cost function the smaller the value of K corresponding to, in which, and are optional scaling vectors whose role is to ensure that the marginal constraint of mass conservation is satisfied, denotes a diagonal matrix form, denotes ; The iteration process is as follows: First, initialize ; then, unroll the alternating iteration over : ； ； wherein, denotes element-wise division, and denotes the scaling vector of the current iteration update; iterating until and the change is less than some threshold, the final solution optimal transport map is: ； For this optimal transmission mapping The optimization for predicting the multiple wavefield is: ； in, A constant is manually set to control the update rate. Indicates the direction of transport, i.e., whether the amplitude decreases or increases; the subscript here... and express and The element index in the middle, For the optimally mapped The elements in the corrected standardized multiple wave field, and the overall updated multiple wave field, need to be inversely standardized to be restored to a state where multiple wave subtraction can be performed; Step c15, specifically comprising the following steps: The collected data for each The element after overall correction of multiple wave field Perform inverse normalization: ； Subtract the multiple wave field corrected finally from the participation multiple wave.

8. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c21, specifically comprising the following steps: In even iteration, the correction phase matching process is performed channel by channel, and the selected object is With The continuous function form of a single channel in the corresponding time slice is and Similarly, due to the positive and negative fluctuations of signal amplitude and the uncertainty of total energy in seismic data, direct application of Wasserstein distance calculation requires data non-negative normalization: ； ； wherein, To prevent the division by zero of a very small positive number, is the minimum value in is the minimum value in is the maximum value in and corresponding to two normalized single channels and are both non-negative and have the same measure, the original channels after non-negative normalization will all be unified to the standard.

9. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c22, specifically comprising the following steps: After standardization, the Wasserstein distance objective function of even-numbered iteration amplitude correction is defined as follows: ； Setting a rigid time-shift mapping , the : shifting in time axis : ； requesting the candidate phase difference , for which no standardization has been performed and , the following exponential coding is performed: ； wherein, is a predetermined normal number to control the upper bound of the power of the numerical value, which is artificially set, and the larger the value, the more obvious the phase feature is; the processed and is the single-channel data highlighting the phase feature.

10. The Marchenko multiple attenuation method based on Wasserstein distance correction of wavefield mismatch according to claim 1, characterized in that, Step c23, specifically comprising the following steps: Computations and The corresponding cumulative distribution function (CDF) and The candidate time shifts are obtained by taking the difference between the inverse CDFs and The set of candidate time shifts is {t} = {t | F(t) - F(t - 1) > 0} : ； ； Next, find an optimal time shift among all candidates which makes the above objective function value minimum corresponding transmission schedule : ； Step c24, specifically comprising the following steps: Existence constraints where, is the maximum movable phase, to guarantee the accuracy of the algorithm here Take half of the source wavelet length , record as the optimized normalized data , then the following inverse normalization process of this data can replace Participate in the convolution correlation operation: 。

Images