A seismic prestack data optimization method and device based on an improved BEMD algorithm

Abstract

The present application belongs to the field of seismic data processing and data optimization, in particular to a method and device for seismic prestack data optimization based on improved BEMD algorithm. The method of the present application uses improved BEMD algorithm and adaptive denoising algorithm to decompose prestack gathers into characteristic signals of different scales. Then, orthogonal wavelet transform denoising based on threshold is carried out on each component to remove most of the noise. Then, the correlation coefficient between each component and the original data is calculated, and the data is reconstructed based on the correlation coefficient. The effective signal is retained to the greatest extent, the interference of noise signal is removed, the signal-to-noise ratio of prestack gathers is improved, and a good data basis is provided for subsequent seismic prediction algorithms.

Description

Technical Field

[0001] The present invention relates to the field of seismic data processing and data optimization, and in particular to a seismic pre-stack data optimization method and device based on an improved BEMD algorithm. Background Art

[0002] Seismic data, as a medium for analyzing subsurface geological conditions, has always held a crucial position. The quality of seismic data directly impacts the accuracy of seismic imaging and inversion. Prestack data is the foundation of all methods. Post-stack seismic data can be derived by superimposing prestack data, which is the data directly obtained after processing the acquired data.

[0003] However, at present, due to the accuracy of the acquisition equipment and the defects of the processing algorithm, pre-stack gathers often have problems such as low signal-to-noise ratio, unclear AVO effect, and discontinuous event axes. In addition, the conventional optimization processing process has many steps and is cumbersome. Coupled with the large scale of pre-stack data, there are often problems such as slow processing speed, large cumulative errors, and weak adaptability. Patent application number CN201410415741.6 provides a pre-stack gather optimization method based on wavelet packet decomposition. The application uses wavelet packets to adaptively decompose seismic gathers, and reconstructs seismic gathers after flattening the event axes in different frequency bands. However, this method does not optimize the pre-stack gathers as a whole. Therefore, a fast, stable, and efficient pre-stack gather optimization method is needed. Summary of the Invention

[0004] The purpose of the present invention is to overcome the problems of slow processing speed, large cumulative errors, and weak adaptability existing in the prior art, and to provide a seismic pre-stack data optimization method and equipment based on an improved BEMD (Bidimensional empirical mode decomposition) algorithm to effectively enhance the internal trend of the pre-stack gather signal and reduce signal noise, thereby providing a better data basis for methods such as pre-stack inversion and improving the accuracy of reservoir prediction.

[0005] In order to achieve the above-mentioned object of the invention, the present invention provides the following technical solutions:

[0006] A seismic pre-stack data optimization method based on an improved BEMD algorithm comprises the following steps:

[0007] S1: input pre-stack gathers;

[0008] S2: performing empirical mode decomposition on the pre-stack gathers using an improved BEMD algorithm to obtain multiple BIMF components (Bidimensional Intrinsic Mode Function) and residual components; wherein the improved BEMD algorithm is a BEMD algorithm that adds erosion and dilation operations of a mathematical morphology algorithm to calculate extreme points and a cubic spline interpolation algorithm to perform envelope fitting;

[0009] S3: performing threshold-based orthogonal wavelet transform denoising on the BIMF component to obtain the denoised BIMF component;

[0010] S4: Obtain the correlation coefficient between the denoised BIMF component and the residual component, perform weighted stacking based on the correlation coefficient, and output the weighted stacked data as the optimized pre-stack gather. The present invention improves the two-dimensional empirical mode algorithm by using a mathematical morphology algorithm and a cubic spline interpolation algorithm to enhance its adaptability. The use of the cubic spline interpolation algorithm for calculation also significantly increases data processing speed. The algorithm then decomposes the pre-stack gather to obtain intrinsic mode components with different scale characteristics, which are more adaptable to the changing characteristics of the pre-stack signal and the waveform's variation with offset. The present invention then denoises the first component using a threshold-based orthogonal wavelet transform method, effectively removing signal noise and retaining the effective signal to the greatest extent possible. Finally, based on coefficientized weighted stacking reconstruction, the correlation coefficient is utilized to maximize the retention of the effective signal and remove the interference of noise signals, providing a good data foundation for subsequent earthquake prediction algorithms. The present invention offers higher computational accuracy and faster speed, making it suitable for processing massive amounts of seismic data.

[0011] As a preferred solution of the present invention, the step S1 further includes performing a pre-analysis on the input pre-stack gathers; the pre-analysis includes the following steps:

[0012] The data quality of the pre-stack data set is judged. When the position of the event axis of the pre-stack data set is unclear, the event axis is in a non-horizontal and discontinuous state, and / or the signal-to-noise ratio is lower than a preset threshold, the pre-stack data set is subjected to denoising, predictive deconvolution, and layer flattening. The quality judgment criteria of the pre-stack data set are: the event axis of equal time points should be kept on a horizontal plane, be in a continuous changing state, and be clearly visible. Therefore, when the continuity is weak (the event axis is in a non-horizontal and discontinuous state) and / or the event axis is unclear (the position of the event axis cannot be seen clearly), and / or the signal-to-noise ratio is lower than a preset threshold, its quality does not meet the expected standard. The present invention effectively avoids the problem of large error in the result caused by directly using low-quality data for calculation by performing quality analysis on the input pre-stack data, and also avoids the problem of excessive workload caused by pre-processing all data. It effectively reduces the workload of the method and improves work efficiency while ensuring data quality.

[0013] As a preferred embodiment of the present invention, step S2 includes the following process:

[0014] S21: Let i=1, r i-1 (x,y)=f(x,y),h i-1 (x,y)=r i-1 (x,y), where r i (x, y) is the matrix to be processed after decomposing the BIMF component of layer i, r0(x, y) = f(x, y), f(x, y) is the pre-stack gather, h i (x,y) is the original matrix after the BIMF component is iterated i times;

[0015] S22: performing erosion and dilation operations on the original matrix using a mathematical morphology algorithm to obtain local maximum points and local minimum points of the original matrix after the BIMF component is iterated i-1 times;

[0016] S23: Using a surface interpolation algorithm based on cubic basis function to perform envelope fitting on the local maximum point and the local minimum point, and calculate the BIMF component Bimf of iteration i i (x, y) and the matrix to be processed after decomposing the i-layer BIMF components;

[0017] S24: Determine whether the preset termination condition is met. If the termination condition is met, proceed to step S25; if the termination condition is not met, set i=i+1 and proceed to step S22;

[0018] S25: Calculate the residual component r N (x, y), and output the decomposed BIMF component and residual component; wherein the residual component r N(x, y) is the matrix to be processed after decomposing the i-layer BIMF components. This paper uses an improved BEMD algorithm to rapidly separate signals. The key lies in obtaining local extreme points and extreme profiles. Therefore, this paper uses a mathematical morphology algorithm to extract extreme points and performs envelope fitting using a surface interpolation algorithm based on cubic basis functions. This also improves the algorithm's stability and adaptability to seismic data.

[0019] As a preferred embodiment of the present invention, step S23 includes:

[0020] S231: Perform envelope fitting on the maximum and minimum points to construct a two-dimensional maximum envelope surface u (i-1)max (x,y) and the two-dimensional minimum envelope surface u (i-1)min (x,y);

[0021] S232: Calculate the two-dimensional maximum envelope surface u (i-1)max (x, y) and the two-dimensional minimum envelope surface u (i-1)min The mean of (x,y);

[0022] S233: Calculate the BIMF component Bimf of iteration i according to the following formula i (x,y) and the matrix to be processed after decomposing the i-layer BIMF component:

[0023] Bimf i (x,y)=h i (x,y)=h i-1 (x,y)-[u (i-1)max (x,y)+u (i-1)min (x,y)] / 2

[0024] r i (x,y)=r i-1 (x,y)-Bimf i (x,y).

[0025] As a preferred embodiment of the present invention, the decomposed BIMF components and residual components are expressed as follows:

[0026]

[0027] n is the number of BIMF components.

[0028] As a preferred embodiment of the present invention, the termination condition is:

[0029] Where r is a preset value.

[0030] As a preferred solution of the present invention, the calculation formula of the optimized pre-stack gather in step S4 is:

[0031]

[0032] Where I(x, y) is the optimized prestack gather, and cor(A·B) represents the correlation coefficient between A and B. The present invention optimizes the reconstruction of prestack gathers by using a coefficientized reconstruction stacking method, which has better adaptability and enhances the adaptive optimization characteristics of the signal, more effectively improving the signal-to-noise ratio of the prestack signal, suppressing noise, and enhancing the AVO characteristics.

[0033] As a preferred embodiment of the present invention, step S3 includes:

[0034] S31: performing orthogonal binary wavelet transform on the BIMF component to obtain the wavelet coefficients and scale coefficients of the BIMF component at each scale;

[0035] S32: selecting a wavelet space threshold, performing threshold processing on the wavelet coefficients using a nonlinear threshold function, and obtaining threshold-processed wavelet coefficients;

[0036] S33: reconstructing the wavelet coefficients and the scaling coefficients after the threshold processing to obtain the denoised BIMF components.

[0037] As a preferred embodiment of the present invention, the BIMF component described in step S3 is the first BIMF component. Since most noise information is contained in the small-scale component, if only the small-scale component information is removed before reconstructing the data, some effective detail information on the small-scale component information will inevitably be lost. Since the first BIMF component is the minimum-scale component, the present invention uses the first BIMF component to be denoised and then reconstructed, effectively performing denoising. The wavelet transform can effectively concentrate the signal energy on a few wavelet bases while still preserving the white noise within the entire wavelet threshold. Therefore, for the wavelet-transformed signal, selecting a more appropriate threshold can maximize the retention of the effective signal coefficient, eliminate the noise coefficient, and then reconstruct the original signal using the wavelet coefficients to achieve the purpose of denoising.

[0038] An electronic device comprises at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to execute any of the methods described above.

[0039] Compared with the prior art, the present invention has the following beneficial effects:

[0040] 1. The present invention improves the two-dimensional empirical mode algorithm, and improves the adaptability of the two-dimensional empirical mode algorithm through the mathematical morphology algorithm and the cubic spline interpolation algorithm. The use of the cubic spline interpolation algorithm for calculation also greatly improves the data processing speed. Then, the pre-stack gather is decomposed by the algorithm to obtain the intrinsic mode components of different scale characteristics, which are more adaptable to the changing characteristics of the pre-stack signal and the changing characteristics of the waveform with the offset distance. The present invention then performs denoising processing on the first component through the threshold-based orthogonal wavelet transform method, effectively removing the noise of the signal and retaining the effective signal to the maximum extent. Finally, based on the coefficientized weighted superposition reconstruction processing, the correlation coefficient is used to retain the effective signal to the maximum extent, remove the interference of the noise signal, and provide a good data foundation for the subsequent earthquake prediction algorithm. The present invention has higher calculation accuracy and fast speed, and is suitable for processing massive seismic data.

[0041] 2. The present invention effectively avoids the problem of large error in the result caused by directly using low-quality data for calculation by performing quality analysis on the input pre-stack gather data, and also avoids the problem of excessive workload caused by pre-processing all data. It effectively reduces the workload of the method while ensuring data quality and improves work efficiency.

[0042] 3. The present invention adopts an improved BEMD algorithm to quickly realize signal separation. The key point is the acquisition of local extreme points and extreme profiles. Therefore, the present invention adopts a mathematical morphology algorithm to extract extreme points and performs envelope fitting through a surface interpolation algorithm based on cubic basis functions. At the same time, it also improves the stability of the algorithm and its adaptability to seismic data.

[0043] 4. The present invention optimizes the reconstruction of pre-stack gathers by using a coefficient-based reconstruction and stacking method, which has better adaptability and improves the adaptive optimization characteristics of the signal, more effectively improving the signal-to-noise ratio of the pre-stack signal, suppressing noise, and enhancing the AVO characteristics.

[0044] 5. The present invention effectively denoises the first BIMF component by performing denoising and subsequent reconstruction. Wavelet transforms effectively concentrate the signal's energy on a few wavelet bases while preserving white noise within the entire wavelet threshold. Therefore, selecting a suitable threshold for the wavelet-transformed signal maximizes the retention of valid signal coefficients while eliminating noise coefficients. The original signal is then reconstructed using the wavelet coefficients. The denoised first BIMF component retains the valid signal while removing most of the noise, effectively suppressing both white noise and random noise. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1This is a flow chart of a seismic pre-stack data optimization method based on an improved BEMD algorithm according to Example 1 of the present invention;

[0046] Figure 2 This is a pre-stack gather map (target layer interval) of a certain group of wellside traces in a certain location in the Sichuan Basin in a seismic pre-stack data optimization method based on an improved BEMD algorithm described in Example 2 of the present invention;

[0047] Figure 3 This is a schematic diagram of pre-stack gathers and BEMD results in a seismic pre-stack data optimization method based on an improved BEMD algorithm according to Example 2 of the present invention;

[0048] Figure 4 Schematic diagram of the denoising result of the first component and the removed noise using a wavelet threshold denoising algorithm based on an orthogonal threshold in a seismic prestack data optimization method based on an improved BEMD algorithm described in Example 2 of the present invention;

[0049] Figure 5 A table of correlation coefficients and coefficientized weight values ​​between each component and the original component in a seismic prestack data optimization method based on an improved BEMD algorithm as described in Example 2 of the present invention;

[0050] Figure 6 Schematic diagram of pre-stack gathers after coefficientized weight optimization in a seismic pre-stack data optimization method based on an improved BEMD algorithm according to Example 2 of the present invention;

[0051] Figure 7 A comparison diagram of gather processing before and after optimization and a schematic diagram of AVO feature analysis in a seismic pre-stack data optimization method based on an improved BEMD algorithm described in Example 2 of the present invention;

[0052] Figure 8 This is an electronic device described in Example 3 of the present invention that utilizes the seismic pre-stack data optimization method based on the improved BEMD algorithm described in Example 1. DETAILED DESCRIPTION

[0053] The present invention will be further described in detail below in conjunction with test examples and specific embodiments. However, this should not be understood as limiting the scope of the present invention to the following embodiments, and all technologies implemented based on the present invention fall within the scope of the present invention.

[0054] The EMD algorithm, or empirical mode decomposition, has been widely used in time-frequency analysis. The BEMD (Bidimensional empirical mode decomposition) algorithm, an adaptive time-frequency localized multiscale analysis algorithm, adaptively decomposes images and noise into a finite number of intrinsic mode function subimages with different characteristic scales based on the inherent characteristics of the data, compared to Fourier transforms and wavelet transforms that rely on prior function bases. The decomposed IMF components not only reveal the true physical information inherent in the noisy image, but also reduce the interference and coupling of feature information between images. Therefore, this application is based on the BEMD algorithm.

[0055] Extract extreme points based on mathematical morphology algorithm:

[0056] Mathematical morphology is an image analysis discipline based on the principles of grunt and topology. It is the fundamental theory of mathematical morphology image processing. By defining structuring elements and performing mathematical morphological operations on matrices, extreme points can be quickly and accurately extracted. The main operations are erosion and dilation.

[0057] Erosion: The erosion of structure A by structural element B is defined as:

[0058]

[0059] It can be understood as moving structure B. If the intersection of structure B and structure A is completely within the area of ​​structure A, then the position point is saved. All points that meet the conditions constitute the result of structure A being eroded by structure B. The value of the output pixel is the minimum value of all input pixel values. In a binary image, if there is a pixel value of 0 in the area, the output pixel value is 0.

[0060] Dilation: The dilation of structure A by structure B is defined as:

[0061]

[0062] This can be understood as performing a convolution operation on structure B. If there is an overlap with structure A during the movement of structure B, that location is recorded. The set of all locations where the movement of structure B intersects with structure A is the result of the expansion of structure A under the influence of structure B. The output pixel value is the maximum value of all input pixel values. In a binary image, if there is a pixel value of 1 in the area, the output pixel value is 1.

[0063] Morphologically based erosion and dilation processes, corresponding to the extraction of maximum and minimum points, respectively, can quickly and efficiently enhance the characteristics of extreme points through structuring element processing, ensuring the accuracy of the extraction results. After a series of experiments and verification, the selection of a cross 4×4 structuring element for morphological operations has a certain effect on protecting and improving both horizontal and vertical resolution, and improving the accuracy of extreme point extraction.

[0064] Surface interpolation algorithm based on cubic basis function:

[0065] This algorithm uses cubic interpolation using the values ​​of 16 points surrounding the sampled point. This method not only considers the influence of the four directly adjacent points but also the influence of the rate of change between the values ​​of each neighboring point. Selecting cubic basis functions, similar to seismic wavelets, allows for a better reconstruction and representation of seismic profile characteristics. In cubic polynomial interpolation, both the interpolation function and its first-order derivative are continuous, resulting in smoother interpolation results and slightly faster computation speed.

[0066] The mathematical expression of the cubic basis function is as follows:

[0067]

[0068] The interpolation formula is as follows:

[0069] f(i+u,j+v)=ABC

[0070] Among them, A, B, and C are all matrices in the following form:

[0071] A=[S(1+u) S(u) S(1-u) S(2-u)]

[0072]

[0073] C=[S(1+v) S(v) S(1-v) S(2-v)] T

[0074] Where f(i, j) represents the value of the original data at (i, j). The obtained point is the interpolation point.

[0075] Surface interpolation based on cubic functions can take multiple factors into account. Using basis functions that match seismic wavelets can more effectively match seismic profile characteristics and ensure the integrity of the maximum (minimum) envelope surface. Furthermore, this algorithm is computationally fast and can quickly construct the maximum (minimum) envelope surface, laying a solid foundation for processing massive amounts of prestack seismic data.

[0076] Threshold-based orthogonal wavelet transform denoising:

[0077] The effectiveness of wavelet threshold denoising depends mainly on three aspects: (1) the selection of the number of signal decomposition layers; (2) the determination of the threshold; and (3) the selection of the threshold processing function. This application mainly refers to the currently more practical method based on white noise test and adaptive determination of the wavelet spatial threshold of each layer based on the 3q criterion, and combines the advantages and disadvantages of hard threshold and soft threshold denoising to comprehensively improve the wavelet coefficient threshold denoising method.

[0078] 1) Selection of the number of signal decomposition layers;

[0079] According to wavelet transform theory, white noise remains white noise after orthogonal wavelet transform, and its energy is mainly distributed in most wavelet spaces. Therefore, in these order wavelet spaces, white noise plays a dominant role, so the wavelet coefficients show obvious white noise characteristics; while the useful signal after wavelet transform, its energy is compressed into a few large-scale wavelet coefficients with large values. The wavelet coefficients of the useful signal dominate, making these wavelet coefficient sequences show non-white noise sequences. Therefore, by judging whether the wavelet space coefficient sequence of each layer has white noise characteristics, we can adaptively determine the reasonable decomposition level, thereby achieving the purpose of removing noise and retaining as much useful signal as possible.

[0080] From the knowledge of random processes, we know that white noise is a purely random process, which is composed of a set of unrelated random variable sequences. The autocorrelation sequence of discrete white noise is:

[0081]

[0082] By using this characteristic of discrete white noise autocorrelation sequence, we can perform white noise test on the wavelet coefficient sequence of each layer decomposition. The test method is as follows:

[0083] Assume that the wavelet coefficient of layer j is WY j,k (k=1,2,...,N j ), N j is the number of wavelet coefficients in the jth layer, and its autocorrelation sequence is ρ i (i=1,2,....,M), M is usually 5-10, if ρ i Satisfy the following formula:

[0084]

[0085] It is believed that WY j,k If it is a white noise sequence, it is necessary to continue the wavelet decomposition until the wavelet coefficient sequence shows a non-white noise sequence and the decomposition is terminated.

[0086] 2) Determination of threshold value;

[0087] Since the wavelet coefficient sequences decomposed from Gaussian white noise at each scale obey Gaussian distribution, that is, W s N(x)~N(μ,σ 2 ). According to the 3σ criterion in mathematical statistics:

[0088] p{-3σ≤W s N(x)-μ≤3σ=0.9974

[0089] In this study, W s N(x)~N(0,σ 2 ). Therefore, the wavelet coefficient sequence that has undergone whitening test can be used to estimate the σ of the white noise sequence by repeatedly calculating the root mean square value σ of the coefficient sequence, and the threshold value is 3σ.

[0090] The high-frequency coefficient WY of each layer decomposition j,k Use the following steps (i.e., the 3σ criterion) to determine the σ value:

[0091] ①Calculate the initial RMS value:

[0092] ②Ask for|WY j,k |(k=1,2,...,N) the maximum value|WY j,k | max , if |WY j,k | max If it is greater than 3σ, it is considered to be the wavelet coefficient of the signal and is eliminated.

[0093] ③Recalculate the RMS value.

[0094] ④Repeat ②, ③, until |WY j,k | max <3σ, then it is considered that the The value of is the estimated value of the white noise sequence σ. Taking 3σ as the threshold for processing the high-frequency coefficients of each layer, substituting it into the threshold function, we can get the optimal estimate of the wavelet coefficients of the original signal. Then, we can perform the inverse wavelet transform to get the denoised signal.

[0095] 3) Selection of threshold processing function.

[0096] According to the different choices of nonlinear threshold function f(WY,t), it can be divided into hard threshold method and soft threshold method. In the hard threshold method, It is discontinuous at t, which will bring some oscillations to the reconstructed signal, so that the reconstructed signal does not have the same smoothness as the original signal; the soft threshold method estimates Although the overall continuity is good, hour, There is always a constant deviation from WY, which reduces the accuracy of the reconstructed signal. Combining the advantages and disadvantages of the hard threshold and soft threshold methods, an improved model for wavelet coefficient threshold estimation is developed:

[0097]

[0098] in:

[0099] σ j =median(|W j,k |).

[0100] When|WY j,k | / t j <1 o'clock, It is the weighted average of the soft threshold and hard threshold estimated wavelet coefficients. The soft threshold estimated coefficients use the "square-difference-square root" processing method, which increases the degree of deviation between the wavelet coefficients and the threshold and promotes signal-noise separation.

[0101] When|WY j,k | / t j >1, in order to minimize the loss of useful signal details and maintain the characteristics of the signal at the singular point, at the scale of j≥j′, the wavelet coefficients of the signal are considered to be absolutely dominant. Even the coefficients smaller than the threshold are considered to contain some information of the signal. Considering the propagation characteristics of the signal wavelet coefficients as the scale increases, the use As the estimated coefficients, other scales Set to zero.

[0102] Example 1

[0103] like Figure 1 As shown, a seismic pre-stack data optimization method based on an improved BEMD algorithm includes the following steps:

[0104] S1: Input pre-stack gathers and perform pre-analysis on the input pre-stack gathers; the pre-analysis includes:

[0105] The data quality of the pre-stack gather is judged. When the events of the pre-stack gather are unclear, the continuity is extremely weak, and / or the signal-to-noise ratio is lower than a preset threshold, the pre-stack gather is subjected to denoising, predictive deconvolution, and layer flattening processing.

[0106] S2: Perform empirical mode decomposition on the pre-stack data set by improving the BEMD method to obtain multiple BIMF components and residual components. Since the core of the BEMD algorithm lies in the determination of the maximum and minimum values ​​of the two-dimensional pre-stack data and the construction of the envelope surface. Therefore, the present invention adds a mathematical morphology algorithm to the extreme point determination in the BEMD algorithm for "erosion" and "dissolution", and then constructs the maximum and minimum envelope surface based on the waveform-based cubic spline interpolation algorithm. The pre-stack data set is regarded as an M×N matrix, that is, f(x,y), x=1,…,M; y=1,…,N, specifically including the following steps:

[0107] S21: Let i=1, r i-1 (x,y)=f(x,y),h i-1 (x,y)=r i-1 (x,y), where r i (x, y) is the matrix to be processed after decomposing the BIMF component of layer i, r0(x, y) = f(x, y), f(x, y) is the pre-stack gather, h i (x,y) is the original matrix after the BIMF component is iterated i times;

[0108] S22: performing erosion and dilation operations on the original matrix using a mathematical morphology algorithm to obtain local maximum points and local minimum points of the original matrix after the BIMF component is iterated i-1 times;

[0109] S23: Using a surface interpolation algorithm based on cubic basis function to perform envelope fitting on the local maximum point and the local minimum point, and calculate the BIMF component Bimf of iteration i i (x, y) and the matrix to be processed after decomposing the i-layer BIMF components;

[0110] S231: Perform envelope fitting on the maximum and minimum points to construct a two-dimensional maximum envelope surface u (i-1)max (x,y) and the two-dimensional minimum envelope surface u (i-1)min (x,y);

[0111] S232: Calculate the two-dimensional maximum envelope surface u (i-1)max (x, y) and the two-dimensional minimum envelope surface u (i-1)min The mean of (x,y);

[0112] S233: Calculate the BIMF component Bimf of iteration i according to the following formula i (x,y) and the matrix to be processed after decomposing the i-layer BIMF component:

[0113] Bimf i (x,y)=h i(x,y)=h i-1 (x,y)-[u (i-1)max (x,y)+u (i-1)min (x,y)] / 2

[0114] r i (x,y)=r i-1 (x,y)-Bimf i (x,y).

[0115] S24: Determine whether the preset termination condition is met. If the termination condition is met, proceed to step S25; if the termination condition is not met, set i=i+1 and proceed to step S22;

[0116] The termination condition is:

[0117] r is a preset value (r is generally set to 0.2-0.3). This application takes 0.2 to ensure the number and quality of IMFs and to ensure that they can better reflect the details of the waveform.

[0118] S25: Calculate the residual component r N (x, y), and output the decomposed BIMF component and residual component; wherein the residual component r N (x,y) is the matrix to be processed after decomposing the BIMF components of layer i. The expressions of the decomposed BIMF components and residual components are:

[0119]

[0120] n is the number of BIMF components.

[0121] S3: Perform threshold-based orthogonal wavelet transform denoising on the first BIMF component. The first BIMF component is the smallest-scale component, which often primarily contains noise because noise is a relatively weak interference signal. The threshold-based orthogonal wavelet transform effectively separates and removes the noise signal.

[0122] S31: the BIMF component Bimf i (x, y) is transformed by orthogonal binary wavelet to obtain the wavelet coefficients w of the BIMF components at each scale. k (i,j) and scale factor h k (i,j);

[0123] S32: Select the wavelet space threshold value, and use the nonlinear threshold function to calculate the wavelet coefficient w k (i, j) is thresholded to obtain the wavelet coefficients after thresholding

[0124] S33: processing the wavelet coefficients after threshold processing and the scale factor h k (i, j) is reconstructed to obtain the BIMF component after denoising.

[0125] S4: Coefficient-based weighted stacking. The correlation coefficients of the denoised BIMF components and the residual components are obtained. Adaptive fitting is then performed based on the magnitude of the correlation coefficients. The weighted stacked data is then output as the optimized prestack gathers. Larger correlation coefficients give higher weights, while smaller ones have the opposite effect. This optimizes the gathers, improves the signal-to-noise ratio, and enhances the internal relationships within the wave group.

[0126] The calculation formula of the optimized pre-stack gather is:

[0127]

[0128] Where I(x,y) is the optimized pre-stack gather, cor(A·B) represents the correlation coefficient between A and B. Example 2

[0129] This embodiment is a specific application example of the method described in Example 1.

[0130] like Figure 2 The following figure shows a pre-stack gather of a wellside trace from a formation in a certain area of ​​the Sichuan Basin (target interval). The left figure shows the pre-stack gather of the wellside trace, and the right figure shows the magnified result. It can be seen that the pre-stack gather is uneven and the signal-to-noise ratio is low.

[0131] like Figure 3 The following figure shows the BEMD plot for the pre-stack gather of the wellside trace. It can be seen that the components derived from the 2D EMD decomposition have different characteristics at different scales. The small-scale IMF1 component contains more detailed information, as well as more noise and other information. The large-scale IMF2 and IMF3 components have good event continuity, effectively mining and reflecting the internal information of the original seismic signal. The AVO characteristics are also more pronounced compared to the original gather.

[0132] like Figure 4Figure 2 shows the denoising results and noise removal of the first component using a wavelet threshold denoising algorithm based on orthogonal thresholds. It can be seen that after denoising, the IMF1 component retains the valid signal while removing most of the noise information, effectively suppressing white noise and random noise. The removed noise profile shows that most of the removed information is useless, with less continuous signal. Although some continuous information is present, this signal is valid information on small-scale components and accounts for a smaller proportion of the original signal. Therefore, the loss of some valid information does not significantly affect the result, while the removal of most of the noise effectively improves the signal-to-noise ratio.

[0133] like Figure 5 The following table shows the correlation coefficients of each component and the original components, as well as the coefficientized weights. The correlation parameters also characterize the characteristics of different components. The weights obtained after coefficientization can better highlight the effective signal and suppress the noise.

[0134] like Figure 6 The figure shows the prestack gathers after coefficientized weight optimization. Compared to the original gathers, the optimized gathers significantly improve the overall signal-to-noise ratio, suppress noise, and enhance the continuity of events. A comparison of the zoomed-in images further demonstrates the excellent optimization results, with some chaotic events being consolidated and compensated.

[0135] like Figure 7 The following figure shows a comparison of gather processing before and after optimization and an AVO feature analysis chart. AVO analysis utilizes the relationship between seismic reflection amplitude and offset (Amplitude-Versus-Offset, AVO for short). Specifically, it analyzes seismic reflections at different offsets in CDP gathers to identify lithology and detect gas content. A comparison of the results shows that the AVO features maintain good consistency before and after optimization. The optimized AVO effect is significantly enhanced, with a more pronounced amplitude-versus-offset variation and a higher degree of feature fit.

[0136] Example 3

[0137] like Figure 8 As shown, an electronic device includes at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the seismic prestack data optimization method based on the improved BEMD algorithm described in the aforementioned embodiment. The input and output interfaces may include a display, a keyboard, a mouse, and a USB interface for inputting and outputting data; and the power supply is used to provide power to the electronic device.

[0138] Those skilled in the art will understand that all or part of the steps of implementing the above-mentioned method embodiment can be completed by hardware related to program instructions, and the aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it executes the steps of the above-mentioned method embodiment; and the aforementioned storage medium includes: mobile storage devices, read-only memories (ROM), magnetic disks or optical disks, and other media that can store program codes.

[0139] When the above-mentioned integrated unit of the present invention is implemented in the form of a software functional unit and sold or used as an independent product, it can also be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the embodiment of the present invention, or the part that contributes to the prior art, can be embodied in the form of a software product, which is stored in a storage medium and includes a number of instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the methods described in each embodiment of the present invention. The aforementioned storage medium includes: various media that can store program codes, such as mobile storage devices, ROMs, magnetic disks, or optical disks.

[0140] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A seismic pre-stack data optimization method based on an improved BEMD algorithm, characterized in that: The following steps are involved: S1: input pre-stack gathers; S2: performing empirical mode decomposition on the pre-stack gathers by using an improved BEMD algorithm to obtain multiple BIMF components and residual components; the improved BEMD algorithm is a BEMD algorithm that adds erosion and dilation operations of a mathematical morphology algorithm to calculate extreme points and a cubic spline interpolation algorithm to perform envelope fitting; S3: performing threshold-based orthogonal wavelet transform denoising on the BIMF component to obtain the denoised BIMF component; S4: obtaining the correlation coefficient between the BIMF component and the residual component after denoising, performing weighted stacking processing according to the correlation coefficient, and outputting the weighted stacked data as the optimized pre-stack gather; Wherein, the step S2 includes the following process: S21: Order , , ,in, To decompose The matrix to be processed after the BIMF component layer, , is the pre-stack gather, Iterate over BIMF components The original matrix after the second time; S22: Using a mathematical morphology algorithm to perform corrosion and expansion operations on the original matrix to obtain the BIMF component iteration The local maximum and local minimum points of the original matrix; S23: Using a surface interpolation algorithm based on cubic basis function to perform envelope fitting on the local maximum point and the local minimum point, and calculating the iteration The BIMF component and decomposition The matrix to be processed after the BIMF component layer; S24: Determine whether the preset termination condition is met. If the termination condition is met, proceed to step S25; if the termination condition is not met, , proceed to step S22; S25: Calculate residual components , and output the decomposed BIMF component and residual component; wherein the residual component To decompose The matrix to be processed after the BIMF component layer; The step S3 comprises: S31: performing orthogonal binary wavelet transform on the BIMF component to obtain the wavelet coefficients and scale coefficients of the BIMF component at each scale; S32: selecting a wavelet space threshold, performing threshold processing on the wavelet coefficients using a nonlinear threshold function, and obtaining threshold-processed wavelet coefficients; S33: reconstructing the wavelet coefficients and the scaling coefficients after the threshold processing to obtain the denoised BIMF components.

2. The seismic prestack data optimization method based on the improved BEMD algorithm according to claim 1, characterized in that: The step S1 further includes performing a pre-analysis on the input pre-stack gathers; the pre-analysis includes the following steps: The data quality of the pre-stack gather is judged. When the event position of the pre-stack gather is unclear, the event is in a non-horizontal and discontinuous state, and / or the signal-to-noise ratio is lower than a preset threshold, the pre-stack gather is subjected to denoising, predictive deconvolution, and layer flattening processing.

3. The seismic prestack data optimization method based on the improved BEMD algorithm according to claim 1, characterized in that: The step S23 includes: S231: Perform envelope fitting on the maximum and minimum points to construct a two-dimensional maximum envelope surface and the two-dimensional minimum envelope surface ; S232: Calculating the two-dimensional maximum envelope surface and the two-dimensional minimum envelope surface The mean of S233: Calculate the iteration according to the following formula The BIMF component and decomposition The matrix to be processed after the BIMF component layer: 。 4. The seismic prestack data optimization method based on the improved BEMD algorithm according to claim 3, characterized in that: The decomposed BIMF components and residual components are expressed as follows: ， is the number of BIMF components.

5. The seismic pre-stack data optimization method based on the improved BEMD algorithm according to claim 1, characterized in that: The termination conditions are: , in, is the default value.

6. The seismic prestack data optimization method based on the improved BEMD algorithm according to claim 1, characterized in that: The calculation formula of the optimized pre-stack gather in step S4 is: ， in, is the optimized pre-stack gather, Indicates seeking A and B The correlation coefficient of .

7. The seismic pre-stack data optimization method based on the improved BEMD algorithm according to claim 1, characterized in that: The BIMF component in step S3 is the first BIMF component.

8. An electronic device, characterized in that: The invention comprises at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to execute the method according to any one of claims 1 to 7.

Images