CN115808713A - Seismic prestack data optimization method and device based on improved BEMD algorithm - Google Patents
Seismic prestack data optimization method and device based on improved BEMD algorithm Download PDFInfo
- Publication number
- CN115808713A CN115808713A CN202111071173.9A CN202111071173A CN115808713A CN 115808713 A CN115808713 A CN 115808713A CN 202111071173 A CN202111071173 A CN 202111071173A CN 115808713 A CN115808713 A CN 115808713A
- Authority
- CN
- China
- Prior art keywords
- algorithm
- bimf
- component
- bemd
- prestack
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Image Processing (AREA)
Abstract
The invention belongs to the field of seismic data processing and data optimization, and particularly relates to a seismic prestack data optimization method and device based on an improved BEMD algorithm. The method decomposes the prestack gather into characteristic signals with different scales by using an improved BEMD algorithm and a self-adaptive denoising algorithm. And then, carrying out threshold-based orthogonal wavelet transform denoising on each component, and removing most of noise. And then reconstructing the data by solving the correlation coefficient between each component and the original data and based on the correlation coefficient. The method has the advantages of reserving effective signals to the maximum extent, removing the interference of noise signals, improving the signal-to-noise ratio of the prestack gather and providing a good data basis for the earthquake follow-up prediction algorithm.
Description
Technical Field
The invention relates to the field of seismic data processing and data optimization, in particular to a seismic prestack data optimization method and device based on an improved BEMD algorithm.
Background
Seismic data has been a vital position for people as a medium for subsurface geology. The quality of the seismic data directly affects the accuracy of seismic imaging and seismic inversion. The prestack data is the basis of all methods, the post-stack seismic data can be stacked from the prestack data, and the prestack data is data directly obtained after the acquired data are processed.
However, in the prior art, due to the accuracy of the acquisition equipment, defects of processing algorithms and other problems, prestack gather often has the problems of low signal-to-noise ratio, unobvious AVO effect, discontinuous phase axis and the like. In addition, the conventional optimization processing flow has more and tedious steps, and the pre-stack data has large scale, so that the problems of low processing speed, more accumulated errors, weak adaptability and the like often exist. The patent with the application number of CN201410415741.6 provides a prestack gather optimization method based on wavelet packet decomposition, wherein wavelet packets are used for conducting adaptive decomposition on a seismic gather, the seismic gather is reconstructed after homophase axis flattening is conducted on different frequency bands, and the method does not optimize the whole prestack gather. There is a need for a fast, stable, and efficient method of prestack gather optimization.
Disclosure of Invention
The invention aims to overcome the problems of low processing speed, more accumulated errors, weak adaptability and the like in the prior art, and provides a seismic prestack data optimization method and equipment based on an improved BEMD (two-dimensional empirical mode decomposition) algorithm, so that the internal trend of prestack gather signals is effectively enhanced, the signal noise is reduced, a better data base can be provided for prestack inversion methods and the like, and the accuracy of reservoir prediction is improved.
In order to achieve the above purpose, the invention provides the following technical scheme:
a seismic prestack data optimization method based on an improved BEMD algorithm comprises the following steps:
s1: inputting a prestack gather;
s2: performing empirical Mode decomposition on the pre-stack gather by improving a BEMD algorithm to obtain a plurality of BIMF (binary Intrinsic Mode functions) and residual components; the improved BEMD algorithm is characterized in that a corrosion and expansion operation calculation extreme point of a mathematical morphology algorithm and a cubic spline interpolation algorithm are added into the BEMD algorithm for envelope fitting;
s3: carrying out threshold-based orthogonal wavelet transform denoising processing on the BIMF component to obtain the denoised BIMF component;
s4: and acquiring a correlation coefficient of the denoised BIMF component and the residual component, performing weighted stacking processing according to the correlation coefficient, and outputting data after the weighted stacking processing as an optimized pre-stack gather. The invention improves the two-dimensional empirical mode algorithm, improves the adaptability of the two-dimensional empirical mode algorithm through a mathematical morphology algorithm and a cubic spline interpolation algorithm, and greatly improves the data processing speed by adopting the cubic spline interpolation algorithm for calculation; and then decomposing the pre-stack trace set by the algorithm to obtain the intrinsic mode components of different scale characteristics, and more adapting to the change characteristics of the pre-stack signals and the change characteristics of the waveform along with the offset. The invention carries out denoising processing on the first component by an orthogonal wavelet transform method based on a threshold value, effectively removes the noise of the signal and furthest reserves the effective signal. And finally, performing weighted stacking reconstruction processing based on the factorization, reserving effective signals to the maximum extent by utilizing the correlation coefficient, removing the interference of noise signals, and providing a good data basis for a seismic follow-up prediction algorithm. The method has higher calculation precision and high speed, and is suitable for processing massive seismic data.
As a preferred embodiment of the present invention, the step S1 further includes performing pre-analysis on the input prestack gather; the pre-analysis comprises the following steps:
and judging the data quality of the prestack gather, and when the prestack gather is unclear in the position of the same-phase axis, the same-phase axis is in an unhevel discontinuous state and/or the signal-to-noise ratio is lower than a preset threshold, carrying out denoising, prediction anti-convolution and layer flattening processing on the prestack gather. The quality evaluation standard of the prestack gather is as follows: the in-phase axes at the equal time points should be kept on a horizontal plane, be in a continuous change state and be clearly visible, so that the quality of the in-phase axes does not meet the expected standard when the continuity is weak (the in-phase axes are in a state of non-horizontal discontinuity) and/or the in-phase axes are not clear (the positions of the in-phase axes are not clearly seen) and/or the signal-to-noise ratio is lower than the preset threshold value. According to the method, the quality analysis is carried out on the input prestack gather data, so that the problem of large result error caused by directly using low-quality data for calculation is effectively avoided, the problem of overlarge workload caused by preprocessing all data is also avoided, the workload of the method is effectively reduced on the premise of ensuring the data quality, and the working efficiency is improved.
As a preferred embodiment of the present invention, the step S2 includes the following steps:
s21: let i =1,r i-1 (x,y)＝f(x,y)，h i-1 (x,y)＝r i-1 (x, y) wherein r i (x, y) is the matrix to be processed after the i-layer BIMF component is decomposed, r 0 (x, y) = f (x, y), f (x, y) is the prestack gather, h i (x, y) is an original matrix after the BIMF component is iterated for i times;
s22: performing corrosion and expansion operation on the original matrix by adopting a mathematical morphology algorithm to obtain a local maximum value point and a local minimum value point of the original matrix of the BIMF component iteration for i-1 times;
s23: envelope fitting is carried out on the local maximum value point and the local minimum value point by adopting a surface interpolation algorithm based on a cubic basis function, and the BIMF component Bimf of iteration i times is calculated i (x, y) and decomposing the i-layer BIMF component to-be-processed matrix;
s24: judging whether a preset termination condition is met, and if the preset termination condition is met, entering the step S25; if the termination condition is not satisfied, let i = i +1, go to step S22;
s25: calculating a residual component r N (x, y) and outputting the decomposed BIMF component and the residual component; wherein the residual component r N And (x, y) is the matrix to be processed after the i-layer BIMF component is decomposed. The invention adopts the improved BEMD algorithm to quickly realize the separation of signals, and the key points of the method are the acquisition of local extreme points and the acquisition of an extreme value profile, so the method adopts the mathematical morphology algorithm to extract the extreme points, carries out envelope fitting through a surface interpolation algorithm based on a cubic basis function, and simultaneously improves the stability of the algorithm and the adaptability to seismic data.
As a preferable embodiment of the present invention, the step S23 includes:
s231: enveloping and fitting the maximum value points and the minimum value points to construct a two-dimensional maximum enveloping surface u (i-1)max (x, y) and two-dimensional minimum envelope surface u (i-1)min (x,y)；
S232: calculating the two-dimensional maximum value enveloping surface u (i-1)max (x, y) and the two-dimensional minimum envelope surface u (i-1)min (x, y) mean;
s233: computing iterations i times according toBIMF component BIMF i (x, y) and decomposing the i-layer BIMF component to be processed matrix:
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
r i (x,y)＝r i-1 (x,y)-Bimf i (x,y)。
as a preferred embodiment of the present invention, the decomposed BIMF component and residual component are expressed as:
n is the number of the BIMF components.
As a preferable embodiment of the present invention, the termination conditions are:
wherein r is a preset value.
As a preferred embodiment of the present invention, the calculation formula of the optimized prestack gather in step S4 is:
wherein, I (x, y) is the optimized prestack gather, and cor (A.B) represents the correlation coefficient of A and B. The invention carries out reconstruction optimization on the pre-stack trace set by a reconstruction superposition method based on the factorization, has better adaptability, improves the self-adaptive optimization characteristic of the signal, more effectively improves the signal-to-noise ratio of the pre-stack signal, suppresses noise and enhances the AVO characteristic.
As a preferable embodiment of the present invention, the step S3 includes:
s31: performing orthogonal dyadic wavelet transform on the BIMF component to obtain wavelet coefficients and scale coefficients of the BIMF component on each scale;
s32: selecting a wavelet space threshold, and performing threshold processing on the wavelet coefficient through a nonlinear threshold function to obtain a wavelet coefficient after threshold processing;
s33: and reconstructing the wavelet coefficient and the scale coefficient after the threshold processing to obtain the denoised BIMF component.
As a preferred embodiment of the present invention, the BIMF component in step S3 is a first BIMF component. Most noise information is contained in the small-scale component, but if the data reconstruction is carried out after the small-scale component information is removed, some effective detail information on the small-scale component information is inevitably lost, and the first BIMF component is the minimum-scale component, so the denoising method effectively carries out denoising after denoising the first BIMF component. While wavelet transform can effectively concentrate the energy of the signal on a few wavelet bases while still preserving white noise in the entire wavelet threshold. Therefore, aiming at the signals after wavelet transformation, a proper threshold value is selected to reserve the coefficient of the effective signals to the maximum extent, the noise coefficient is removed, and then the original signals are reconstructed through the wavelet coefficient, so that the purpose of denoising is achieved.
An electronic device comprising at least one processor, and a memory communicatively coupled to the at least one processor; the memory stores instructions executable by the at least one processor to enable the at least one processor to perform any of the methods described above.
Compared with the prior art, the invention has the beneficial effects that:
1. the invention improves the two-dimensional empirical mode algorithm, improves the adaptability of the two-dimensional empirical mode algorithm through a mathematical morphology algorithm and a cubic spline interpolation algorithm, and greatly improves the data processing speed by adopting the cubic spline interpolation algorithm for calculation; and then decomposing the pre-stack trace set by the algorithm to obtain the intrinsic mode components of different scale characteristics, and more adapting to the change characteristics of the pre-stack signals and the change characteristics of the waveform along with the offset. The invention further carries out denoising processing on the first component by an orthogonal wavelet transform method based on a threshold value, effectively removes the noise of the signal and furthest reserves the effective signal. And finally, performing weighted stacking reconstruction processing based on the factorization, reserving effective signals to the maximum extent by utilizing the correlation coefficient, removing the interference of noise signals, and providing a good data basis for a seismic follow-up prediction algorithm. The method has higher calculation precision and high speed, and is suitable for processing massive seismic data.
2. According to the method, the quality analysis is carried out on the input prestack gather data, so that the problem of large result error caused by directly using low-quality data for calculation is effectively avoided, the problem of overlarge workload caused by preprocessing all data is also avoided, the workload of the method is effectively reduced on the premise of ensuring the data quality, and the working efficiency is improved.
3. The invention adopts the improved BEMD algorithm to quickly realize the separation of signals, and the key points of the method are the acquisition of local extreme points and the acquisition of an extreme value profile, so the method adopts the mathematical morphology algorithm to extract the extreme points, carries out envelope fitting through a surface interpolation algorithm based on a cubic basis function, and simultaneously improves the stability of the algorithm and the adaptability to seismic data.
4. The invention carries out reconstruction optimization on the pre-stack trace set by a reconstruction superposition method based on the factorization, has better adaptability, improves the self-adaptive optimization characteristic of the signal, more effectively improves the signal-to-noise ratio of the pre-stack signal, suppresses noise and enhances the AVO characteristic.
5. The invention effectively carries out denoising by reconstructing the first BIMF component after denoising. While wavelet transform can effectively concentrate the energy of the signal on a few wavelet bases while still preserving white noise in the entire wavelet threshold. Therefore, aiming at the signals after wavelet transformation, a proper threshold value is selected to reserve the coefficient of the effective signals to the maximum extent, the noise coefficient is removed, then the original signals are reconstructed through the wavelet coefficient, most noise information is removed on the basis that the effective signals are reserved through the first BIMF component after denoising processing, and white noise and random noise are effectively suppressed.
Drawings
Fig. 1 is a schematic flow chart of a seismic prestack data optimization method based on an improved BEMD algorithm according to embodiment 1 of the present invention;
fig. 2 is a diagram (target interval) of a well side-track prestack gather of a certain group of a place in the sikawa basin in the seismic prestack data optimization method based on the improved BEMD algorithm according to embodiment 2 of the present invention;
fig. 3 is a schematic diagram of a prestack gather and a BEMD result in the seismic prestack data optimization method based on the improved BEMD algorithm according to embodiment 2 of the present invention;
fig. 4 is a schematic diagram of a denoising result and removed noise of a first component by using a wavelet threshold denoising algorithm based on an orthogonal threshold in the seismic prestack data optimization method based on the improved BEMD algorithm according to embodiment 2 of the present invention;
fig. 5 is a correlation coefficient table and a factorization weight value of each component and an original component in the seismic prestack data optimization method based on the improved BEMD algorithm according to embodiment 2 of the present invention;
fig. 6 is a schematic diagram of a prestack gather optimized based on a factorization weight in a seismic prestack data optimization method based on an improved BEMD algorithm according to embodiment 2 of the present invention;
FIG. 7 is a schematic diagram of a comparison graph before and after optimization of gather processing and an AVO characteristic analysis in a seismic prestack data optimization method based on an improved BEMD algorithm according to embodiment 2 of the present invention;
fig. 8 is an electronic device according to embodiment 3 of the present invention, which utilizes the seismic prestack data optimization method based on the improved BEMD algorithm described in embodiment 1.
Detailed Description
The present invention will be described in further detail with reference to test examples and specific embodiments. It should be understood that the scope of the above-described subject matter is not limited to the following examples, and any techniques implemented based on the disclosure of the present invention are within the scope of the present invention.
The EMD algorithm, namely the empirical mode decomposition algorithm, is widely developed in the aspect of time-frequency analysis and has more applications; compared with Fourier change and wavelet change which depend on prior function base, the method adaptively decomposes the image and noise into a finite number of intrinsic mode function sub-images with different characteristic scales step by step according to the intrinsic characteristics of data, and the decomposed IMF component not only reveals the intrinsic real physical information of the noisy image, but also reduces the interference and coupling of the characteristic information between the images. The present application is based on the BEMD algorithm.
Extracting extreme points based on a mathematical morphology algorithm:
mathematical morphology (Mathematical morphology) is an image analysis subject based on the theory of lattice and topology, and is the basic theory of Mathematical morphology image processing. The matrix is subjected to mathematical morphology operation by defining structural elements, so that the extraction of extreme points can be quickly and accurately realized, and the method mainly comprises the following steps: corrosion and swelling operations.
Corrosion (oxidation): the definition of the corrosion of structure a by structural element B is:
it will be understood that moving structure B, if the intersection of structure B and structure a falls completely within the area of structure a, holds the location point, and all points that satisfy the condition constitute the result of structure a being eroded by structure B. The value of the output pixel is the minimum of all the input pixel values, and in a binary image, if one pixel value in the field is 0, the output pixel value is 0.
Expansion (deformation): the definition of structure a expanded by structure B is:
it can be understood that, when the convolution operation is performed on the structure B on the structure a, if there is an overlapping area with the structure a during the movement of the structure B, the position is recorded, and the set of all the positions where the movement of the structure B intersects with the structure a is the expansion result of the structure a under the action of the structure B. The value of the output pixel is the maximum of all the input pixel values. In the binary image, if one pixel value is 1 in the field, the output pixel value is 1.
The method is characterized in that the method respectively corresponds to the extraction of the maximum value point and the extraction of the minimum value point based on the corrosion and expansion treatment of morphology, and can quickly and efficiently enhance the characteristics of the maximum value point and ensure the accuracy of the extraction result through the treatment of structural elements. Through a series of experiments and verifications, the structural elements with the cross shape of 4 multiplied by 4 are selected to carry out morphological operation on the structural elements, so that the transverse resolution and the longitudinal resolution are protected and improved to a certain extent, and the accuracy of extreme point extraction is improved.
The surface interpolation algorithm based on the cubic basis function is as follows:
the algorithm utilizes the values of 16 points around a to-be-sampled point to perform cubic interpolation, not only considers the influence of 4 directly adjacent points, but also considers the influence of the change rate between the adjacent points. The cubic basis function similar to the seismic wavelet is selected, and the characteristics of the seismic section can be better restored and embodied. In the cubic polynomial interpolation method, the interpolation function and the first derivative thereof are continuous, so that the interpolation result is smooth, and the operation speed is slightly high.
The mathematical expression of the cubic basis function is as follows:
the interpolation formula is as follows:
f(i+u,j+v)＝ABC
wherein A, B, C are matrices of the following form:
A＝[S(1+u) S(u) S(1-u) S(2-u)]
C＝[S(1+v) S(v) S(1-v) S(2-v)] T
where f (i, j) represents the value at the original data (i, j). The found point is the interpolation point.
The cubic function-based surface interpolation can consider various factors, and the basis function matched with the seismic wavelet can be used for more strongly matching the seismic section characteristics, so that the integrity of a maximum (small) envelope surface is ensured. The algorithm is high in calculation speed, the construction of the envelope surface with the maximum (small) value can be rapidly realized, and a good foundation is laid for the processing of the pre-stack seismic mass data.
Denoising based on threshold orthogonal wavelet transform:
the quality of the wavelet threshold denoising effect mainly depends on three aspects: (1) selecting the number of signal decomposition layers; (2) determining a threshold value; and (3) selecting a threshold processing function. The method mainly refers to a practical self-adaptive determination method based on white noise inspection and 3q criterion of wavelet space threshold values of all layers at present, and combines the advantages and disadvantages of hard threshold value denoising and soft threshold value denoising.
1) Selecting the number of signal decomposition layers;
the wavelet transform theory shows that white noise is still white noise after orthogonal wavelet transform, and the energy of the white noise is mainly distributed in most wavelet spaces, so that the white noise plays a dominant role in the wavelet spaces of the orders, and the wavelet coefficient has obvious white noise characteristics; after the useful signal is wavelet-transformed, its energy is compressed into a few wavelet coefficients with larger value in large-scale wavelet space, and the wavelet coefficients of the useful signal are dominant, so that these wavelet coefficient sequences represent non-white noise sequences. Therefore, by judging whether the wavelet space coefficient sequences of each layer have white noise characteristics or not, reasonable decomposition levels can be determined in a self-adaptive mode, and the purposes of removing noise and reserving useful signals as much as possible are achieved.
From the knowledge of the random process, white noise is a purely random process, which is composed of a set of random variable sequences that are not related to each other. The autocorrelation sequence of the discrete white noise is:
the white noise test can be carried out on the wavelet coefficient sequence decomposed by each layer by utilizing the characteristic of the discrete white noise autocorrelation sequence. The test method comprises the following steps:
let the wavelet coefficient of the j-th layer be WY j,k (k＝1,2,...,N j )，N j The number of wavelet coefficients of the j-th layer is the autocorrelation sequence of rho i (i =1,2,.. Times, M), M is typically 5-10, if ρ i Satisfies the following formula:
then consider WY j,k For a white noise sequence, the wavelet decomposition is continued until the wavelet coefficient sequence shows a non-white noise sequence, and the decomposition is terminated.
2) Determining a threshold value;
the wavelet coefficient sequence obtained by decomposing Gaussian white noise on each scale is subjected to Gaussian distribution, namely W s N(x)～N(μ,σ 2 ). As can be seen from the 3 σ criterion in mathematical statistics:
p{-3σ≤W s N(x)-μ≤3σ}＝0.9974
w in this study s N(x)～N(0,σ 2 ). Therefore, the wavelet coefficient sequence after whitening test can estimate the sigma of the white noise sequence by repeatedly calculating the root mean square value sigma of the coefficient sequence, and the threshold value is 3 sigma.
High frequency coefficient WY decomposed for each layer j,k The σ value is determined using the following steps (i.e., the 3 σ criterion):
(1) calculating an initial root mean square value:
(2) obtaining | WY j,k Maximum value | WY of | (k =1,2,.., N) j,k | max If | WY j,k | max If the signal is more than 3 sigma, the wavelet coefficient is regarded as the wavelet coefficient of the signal, and the wavelet coefficient is removed.
(3) The root mean square value is recalculated.
(4) Repeating (2) and (3) until the absolute value WY j,k | max < 3 σ, it is considered that
The value of (c) is an estimated value of the white noise sequence σ. And taking 3 sigma as a threshold for processing each layer of high-frequency coefficient, substituting the threshold into a threshold function to obtain the optimal estimation of the wavelet coefficient of the original signal, and performing wavelet inverse transformation to obtain the signal after denoising.
3) Selection of a thresholding function.
The non-linear threshold function f (WY, t) can be classified into a hard threshold method and a soft threshold method according to their selection. In the hard threshold method, in the case of the hard threshold method,
discontinuity at t may give some oscillation to the reconstructed signal, so that the reconstructed signal does not have the desired smoothness as the original signal; and estimated by soft threshold method
Although the overall continuity is good, when
When the temperature of the water is higher than the set temperature,
there is always a constant deviation from WY, reducing the accuracy of the reconstructed signal. Combining the advantages and disadvantages of the hard threshold and soft threshold methods, an improved wavelet coefficient threshold estimation model is developed:
wherein:
when | WY j,k |/t j When the ratio is less than 1, the reaction solution is,
the method is a weighted average of wavelet coefficients estimated by a soft threshold and a hard threshold, wherein the soft threshold estimation coefficients adopt a square-difference-evolution processing method, so that the deviation degree of the wavelet coefficients and the threshold is increased, and the signal-noise separation is promoted.
When | WY j,k |/t j To reduce the loss of detail information of useful signal as much as possible, and to keep the characteristics of signal at singular point, on the scale of j ≧ j', the wavelet coefficient of signal is considered to be absolutely dominant, even if the coefficient is less than threshold, some information of signal is considered to be contained in it, then the propagation characteristic of wavelet coefficient of signal increasing with scale is considered, the method uses
As estimated coefficients, on other scales
And (5) setting to zero.
Example 1
As shown in fig. 1, a seismic prestack data optimization method based on the improved BEMD algorithm includes the following steps:
s1: inputting a prestack gather, and performing pre-analysis on the input prestack gather; the pre-analysis comprises:
and judging the data quality of the prestack gather, and when the same phase axis of the prestack gather is unclear, the continuity is extremely weak and/or the signal-to-noise ratio is lower than a preset threshold value, carrying out denoising, prediction anti-folding and layer flattening processing on the prestack gather.
S2: and performing empirical mode decomposition on the prestack gather by improving a BEMD method to obtain a plurality of BIMF components and residual components. The core of the BEMD algorithm lies in the maximum and minimum value solving and the envelope surface constructing of the two-dimensional pre-stack data. Therefore, the invention adds mathematical morphology algorithm to carry out erosion and dissolution in the extreme point calculation of BEMD algorithm, and then carries out the construction of maximum and minimum enveloping surface based on the cubic spline interpolation algorithm of waveform mode. Considering the prestack gather as an M × N matrix, i.e., f (x, y), x =1, …, M; y =1, …, N, comprising the steps of:
s21: let i =1,r i-1 (x,y)＝f(x,y)，h i-1 (x,y)＝r i-1 (x, y) wherein r i (x, y) is the matrix to be processed after the i-layer BIMF component is decomposed, r 0 (x, y) = f (x, y), f (x, y) is the prestack gather, h i (x, y) is an original matrix after the BIMF component is iterated for i times;
s22: carrying out corrosion and expansion operation on the original matrix by adopting a mathematical morphology algorithm to obtain a local maximum value point and a local minimum value point of the original matrix of the BIMF component iteration for i-1 times;
s23: envelope fitting is carried out on the local maximum value point and the local minimum value point by adopting a surface interpolation algorithm based on a cubic basis function, and the BIMF component Bimf of iteration i times is calculated i (x, y) and decomposing the i-layer BIMF component to-be-processed matrix;
s231: enveloping and fitting the maximum value points and the minimum value points to construct a two-dimensional maximum enveloping surface u (i-1)max (x, y) and two-dimensional minimum envelope surface u (i-1)min (x,y)；
S232: calculating the two-dimensional maximum value enveloping surface u (i-1)max (x, y) and the two-dimensional minimum envelope surface u (i-1)min (x, y) mean;
s233: calculating BIMF component Bimf of iteration i times according to the following formula i (x, y) and the matrix to be processed after decomposing the i-layer BIMF component:
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
r i (x,y)＝r i-1 (x,y)-Bimf i (x,y)。
s24: judging whether a preset termination condition is met, and if the preset termination condition is met, entering the step S25; if the termination condition is not satisfied, let i = i +1, go to step S22;
wherein the termination condition is:
S25: calculating a residual component r N (x, y) and outputting the decomposed BIMF component and the residual component; wherein the residual component r N (x, y) is a matrix to be processed after the i-layer BIMF component is decomposed, and the expressions of the decomposed BIMF component and the residual component are as follows:
n is the number of the BIMF components.
S3: and carrying out threshold-based orthogonal wavelet transform denoising processing on the first BIMF component. The first BIMF component is the smallest scale component, and often the small scale component is the component that contains mainly noise, since noise is a relatively weak interfering signal. And the orthogonal wavelet transform processing based on the threshold value can well separate noise signals and effectively remove the noise signals.
S31: for the BiMF component Bimf i (x, y) orthogonal dyadic wavelet transform to obtain wavelet coefficient w of the BIMF component in each scale k (i, j) and a scale factor h k (i,j)；
S32: selecting a wavelet space threshold value, and applying a nonlinear threshold function to the wavelet coefficient w k (i, j) carrying out threshold processing to obtain wavelet coefficient after threshold processing
S33: wavelet coefficient after threshold processing
And the scale factor h k And (i, j) reconstructing to obtain the denoised BIMF component.
S4: and (5) performing coefficient weighted superposition processing. And acquiring a correlation coefficient of the denoised BIMF component and the residual component, performing self-adaptive fitting through the magnitude relation of the correlation coefficient, and outputting data after weighted stacking processing as an optimized pre-stack gather. A large correlation coefficient has a large specific gravity, whereas a small correlation coefficient has the opposite. Therefore, the optimization processing of the gather is realized, the signal-to-noise ratio is improved, and the internal relation of the wave group is enhanced.
Wherein, the calculation formula of the optimized prestack gather is as follows:
where I (x, y) is the optimized prestack gather, cor (A.B) denotes the correlation coefficient of A and B, example 2
This embodiment is a specific application example of the method described in embodiment 1.
As shown in fig. 2, it is a diagram of a set of well side-track prestack (target interval) in a certain group of the site of the sichuan basin, wherein the left diagram is the diagram of the well side-track prestack, and the right diagram is the enlarged result. It can be seen that there are cases where the gather is uneven and the signal-to-noise ratio is low in the prestack gather.
As shown in fig. 3, a BEMD map of a well side-track prestack gather. It can be seen that the components obtained by the two-dimensional EMD decomposition have different characteristics at different scales. The IMF1 component under the small scale contains more detailed information, the noise and other information also contain more, the IMF2 and IMF3 components under the large scale have good continuity of the same phase axis, the internal information of the original seismic signal is effectively mined and reflected, and the AVO characteristic is more obvious compared with the original gather.
As shown in fig. 4, the denoising result and the removed noise of the first component are obtained by using the wavelet threshold denoising algorithm based on the orthogonal threshold. Therefore, most noise information is removed on the basis that the effective signal is reserved for the IMF1 component subjected to denoising processing, and white noise and random noise are effectively suppressed. It can be seen from the removed noise profile that most of the removed noise is useless information, the continuous signal is less, and although the signal also contains a part of continuous information, because the signal is effective information on a small-scale component, the proportion of the signal to the original signal is smaller, the loss of part of the effective information does not bring great influence to the result, and the removal of most of the noise information effectively improves the signal-to-noise ratio.
As shown in fig. 5, the correlation coefficient table and the scaling weight value of each component and the original component are shown. The related parameters also represent the characteristics of different components, and the weighted values obtained after the coefficient processing can highlight effective signals and suppress noise.
As shown in fig. 6, is the prestack gather optimized based on the factorized weights. Compared with the original gather, the optimized gather has the advantages that the integral signal-to-noise ratio is effectively improved, noise is suppressed, and the continuity of the same phase axis is stronger. From the partial enlarged image comparison, the optimization result is better, and some disordered same-phase axes are also rounded and compensated.
As shown in FIG. 7, a pre-and post-optimization comparison graph and an AVO feature analysis graph are processed for a gather. The AVO analysis technology is a seismic technology which utilizes the relation between seismic reflection Amplitude and Offset change (AVO for short), namely, the lithology is identified and the gas content is detected by analyzing the seismic reflection of different Offset in a CDP channel set. The result comparison shows that the AVO characteristics before and after optimization keep better consistency, the optimized AVO effect is effectively enhanced, the amplitude changes along with the offset distance are more obvious, and the characteristic fitting degree is higher.
Example 3
As shown in fig. 8, an electronic device includes at least one processor, and a memory communicatively coupled to the at least one processor; the memory stores instructions executable by the at least one processor to enable the at least one processor to perform a method of seismic prestack data optimization based on the modified BEMD algorithm as described in the previous embodiments. The input and output interface can comprise a display, a keyboard, a mouse and a USB interface and is used for inputting and outputting data; the power supply is used for supplying electric energy to the electronic equipment.
Those skilled in the art will understand that: all or part of the steps for realizing the method embodiments can be completed by hardware related to program instructions, the program can be stored in a computer readable storage medium, and the program executes the steps comprising the method embodiments when executed; and the aforementioned storage medium includes: various media that can store program codes, such as a removable Memory device, a Read Only Memory (ROM), a magnetic disk, or an optical disk.
When the integrated unit of the present invention is implemented in the form of a software functional unit and sold or used as a separate product, it may also be stored in a computer-readable storage medium. Based on such understanding, the technical solutions of the embodiments of the present invention may be essentially implemented or a part contributing to the prior art may be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device) to execute all or part of the methods described in the embodiments of the present invention. And the aforementioned storage medium includes: a removable storage device, a ROM, a magnetic or optical disk, or other various media that can store program code.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.
Claims (10)
1. A seismic prestack data optimization method based on an improved BEMD algorithm is characterized by comprising the following steps:
s1: inputting a prestack gather;
s2: performing empirical mode decomposition on the pre-stack gather by improving a BEMD algorithm to obtain a plurality of BIMF components and residual components; the improved BEMD algorithm is characterized in that a corrosion and expansion operation calculation extreme point of a mathematical morphology algorithm and a cubic spline interpolation algorithm are added into the BEMD algorithm for envelope fitting;
s3: carrying out threshold-based orthogonal wavelet transform denoising processing on the BIMF component to obtain the denoised BIMF component;
s4: and acquiring a correlation coefficient of the denoised BIMF component and the residual component, performing weighted stacking processing according to the correlation coefficient, and outputting data after the weighted stacking processing as an optimized pre-stack gather.
2. The seismic prestack data optimization method based on the improved BEMD algorithm as claimed in claim 1, wherein the step S1 further comprises pre-analyzing the input prestack gather; the pre-analysis comprises the following steps:
and judging the data quality of the prestack gather, and when the prestack gather is unclear in the position of the same-phase axis, the same-phase axis is in an unhevel discontinuous state and/or the signal-to-noise ratio is lower than a preset threshold, carrying out denoising, prediction anti-convolution and layer flattening processing on the prestack gather.
3. The seismic prestack data optimization method based on the improved BEMD algorithm as claimed in claim 1, wherein the step S2 comprises the following procedures:
s21: let i =1,r i-1 (x,y)＝f(x,y)，h i-1 (x,y)＝r i-1 (x, y) wherein r i (x, y) is the matrix to be processed after the i-layer BIMF component is decomposed, r 0 (x, y) = f (x, y), f (x, y) is the prestack gather, h i (x, y) is an original matrix after the BIMF component is iterated for i times;
s22: carrying out corrosion and expansion operation on the original matrix by adopting a mathematical morphology algorithm to obtain a local maximum value point and a local minimum value point of the original matrix of the BIMF component iteration for i-1 times;
s23: envelope fitting is carried out on the local maximum value point and the local minimum value point by adopting a surface interpolation algorithm based on a cubic basis function, and the BIMF component Bimf of iteration i times is calculated i (x, y) and decomposing the i-layer BIMF component to-be-processed matrix;
s24: judging whether a preset termination condition is met, and if the preset termination condition is met, entering the step S25; if the termination condition is not satisfied, let i = i +1, go to step S22;
s25: calculating a residual component r N (x, y) and outputting the decomposed BIMF component and the residual component; wherein the residual component r N And (x, y) is the matrix to be processed after the i-layer BIMF component is decomposed.
4. The method for optimizing seismic prestack data based on the improved BEMD algorithm as claimed in claim 3, wherein the step S23 comprises:
s231: enveloping and fitting the maximum value points and the minimum value points to construct a two-dimensional maximum enveloping surface u (i-1)max (x, y) and two-dimensional minimum envelope surface u (i-1)min (x,y)；
S232: calculating the two-dimensional maximum value enveloping surface u (i-1)max (x, y) and the two-dimensional minimum envelope surface u (i-1)min (x, y) mean;
s233: calculating the BIMF component Bimf iterated for i times according to the following formula i (x, y) and the matrix to be processed after decomposing the i-layer BIMF component:
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
r i (x,y)＝r i-1 (x,y)-Bimf i (x,y)。
5. the seismic prestack data optimization method based on the improved BEMD algorithm as claimed in claim 4, wherein the decomposed BIMF component and residual component expressions are:
n is the number of the BIMF components.
6. The seismic prestack data optimization method based on the improved BEMD algorithm as claimed in claim 3, wherein the termination condition is:
wherein r is a preset value.
7. The seismic prestack data optimization method based on the improved BEMD algorithm as claimed in claim 3, wherein the calculation formula of the optimized prestack gather in the step S4 is as follows:
wherein, I (x, y) is the optimized prestack gather, and cor (A.B) represents the correlation coefficient of A and B.
8. The method for optimizing seismic prestack data based on the improved BEMD algorithm as claimed in claim 1, wherein the step S3 comprises:
s31: performing orthogonal dyadic wavelet transform on the BIMF component to obtain wavelet coefficients and scale coefficients of the BIMF component on each scale;
s32: selecting a wavelet space threshold, and performing threshold processing on the wavelet coefficient through a nonlinear threshold function to obtain a wavelet coefficient after threshold processing;
s33: and reconstructing the wavelet coefficient after threshold processing and the scale coefficient to obtain the denoised BIMF component.
9. The seismic prestack data optimization method based on the improved BEMD algorithm as claimed in claim 1, wherein the BIMF component in step S3 is the first BIMF component.
10. An electronic device comprising at least one processor, and a memory communicatively coupled to the at least one processor; the memory stores instructions executable by the at least one processor to enable the at least one processor to perform the method of any one of claims 1 to 9.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111071173.9A CN115808713A (en) | 2021-09-13 | 2021-09-13 | Seismic prestack data optimization method and device based on improved BEMD algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111071173.9A CN115808713A (en) | 2021-09-13 | 2021-09-13 | Seismic prestack data optimization method and device based on improved BEMD algorithm |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115808713A true CN115808713A (en) | 2023-03-17 |
Family
ID=85481361
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111071173.9A Pending CN115808713A (en) | 2021-09-13 | 2021-09-13 | Seismic prestack data optimization method and device based on improved BEMD algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115808713A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116819627A (en) * | 2023-06-30 | 2023-09-29 | 中海石油(中国)有限公司深圳分公司 | Method, device, equipment and medium for enhancing weak earthquake signal |
| CN118145825A (en) * | 2024-03-07 | 2024-06-07 | 北京顺政水环境有限公司 | Sewage treatment method and system |
-
2021
- 2021-09-13 CN CN202111071173.9A patent/CN115808713A/en active Pending
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116819627A (en) * | 2023-06-30 | 2023-09-29 | 中海石油(中国)有限公司深圳分公司 | Method, device, equipment and medium for enhancing weak earthquake signal |
| CN118145825A (en) * | 2024-03-07 | 2024-06-07 | 北京顺政水环境有限公司 | Sewage treatment method and system |
| CN118145825B (en) * | 2024-03-07 | 2024-08-13 | 北京顺政水环境有限公司 | Sewage treatment method and system |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106204467B (en) | Image denoising method based on cascade residual error neural network | |
| CN112946749B (en) | Method for suppressing seismic multiples based on data augmentation training deep neural network | |
| CN115808713A (en) | Seismic prestack data optimization method and device based on improved BEMD algorithm | |
| CN110058307B (en) | Full waveform inversion method based on fast quasi-Newton method | |
| CN111337980B (en) | Carbon dioxide sequestration monitoring method and system based on time-shift full-waveform inversion | |
| CN113238189A (en) | Sound source identification method and system based on array measurement and sparse prior information | |
| AU2019406345B2 (en) | Methods and systems for calibrating depth in a well to seismic data in a subsurface volume of interest | |
| CN114418886B (en) | Robust denoising method based on depth convolution self-encoder | |
| Wu et al. | An effective approach for underwater sonar image denoising based on sparse representation | |
| CN106443787A (en) | Prestack seismic gather noise suppression method and device | |
| CN108230280A (en) | Image speckle noise minimizing technology based on tensor model and compressive sensing theory | |
| CN108428221A (en) | A kind of neighborhood bivariate shrinkage function denoising method based on shearlet transformation | |
| CN109212608B (en) | Borehole microseismic signal antinoise method based on 3D shearlet transformation | |
| Sun et al. | Full-waveform inversion using a learned regularization | |
| Huang et al. | Multiplicative noise removal based on unbiased box-cox transformation | |
| CN113204051A (en) | Low-rank tensor seismic data denoising method based on variational modal decomposition | |
| CN112630824B (en) | Discrete point spread function generation method and system in seismic imaging | |
| CN110261912B (en) | Interpolation and denoising method and system for seismic data | |
| CN107704831A (en) | A kind of gas density data noise reduction based on singular value decomposition median method | |
| CN113031072B (en) | Multiple wave pressing method, device and equipment between virtual phase axis layers | |
| CN102509268A (en) | Immune-clonal-selection-based nonsubsampled contourlet domain image denoising method | |
| CN113009564A (en) | Seismic data processing method and device | |
| CN118311669B (en) | Distributed optical fiber sensing coupling noise forward modeling and suppressing method | |
| CN112363217A (en) | Random noise suppression method and system for seismic data | |
| CN110687605A (en) | Improved K-SVD algorithm-based algorithm analysis application in seismic signal processing |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |