Fourier finite difference deconvolution imaging method based on inclination angle adaptive noise suppression

Abstract

The invention discloses a Fourier finite difference deconvolution imaging method based on inclination angle adaptive noise suppression, and relates to the technical field of geophysical exploration of petroleum. The method specifically comprises the following steps: designing a scattering point speed model to solve a partial point spread function of an underground space; solving a partial point spread function by a speed disturbance method; solving point spread functions of all positions of the underground space by carrying out Gaussian filtering technology and inverse distance weighted interpolation on part of the point spread functions; using a plane wave destruction method to solve the inclination angle attribute of the work area, and constructing an inclination angle self-adaptive diagonal diffusion Fourier finite difference migration method based on the inclination angle attribute to solve a high-quality initial image suitable for an inversion system; and performing multi-dimensional deconvolution in the wave number domain to obtain a final migration imaging result. According to the method, the point spread function is successfully introduced into Fourier finite difference imaging, and the resolution and deep amplitude of a conventional Fourier finite difference imaging method can be effectively improved.

Description

Technical Field

[0001] This invention relates to the field of petroleum geophysical exploration technology, and in particular to a Fourier finite difference deconvolution imaging method based on tilt-adaptive noise suppression. Background Technology

[0002] Conventional Fourier finite-difference imaging (FTEMI) achieves an excellent balance between accuracy and efficiency and is widely used in petroleum geophysical exploration. However, factors such as acquisition noise and uneven underground lighting can still lead to insufficient accuracy in the imaging results. While data-domain least squares methods based on the linear inversion concept can effectively compensate for these shortcomings of conventional Fourier finite-difference migration, their enormous computational cost limits their application in practical production.

[0003] Therefore, in order to better apply the Fourier finite difference imaging method in the field of oil and gas exploration, it is particularly important to control the computational cost while improving the effect of the conventional Fourier finite difference imaging method. It is urgent to study a better Fourier finite difference inversion imaging method. Summary of the Invention

[0004] The purpose of this invention is to improve the performance of conventional Fourier finite difference imaging methods and control the overall computational cost. A Fourier finite difference deconvolution imaging method based on tilt angle adaptive noise suppression is proposed.

[0005] The present invention specifically adopts the following technical solution:

[0006] A Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression includes:

[0007] (1) Input velocity field, frequency domain source wavelet, and observation system parameters;

[0008] (2) Construct the velocity field of equally spaced scattering points based on the velocity field, and use the Fourier finite difference continuation method to obtain the point spread function at these scattering points;

[0009] (3) Apply Gaussian filtering to the obtained partial point spread function to suppress artifacts far from the center point;

[0010] (4) The point spread function of the entire work area is obtained by applying the inverse distance weighted interpolation method;

[0011] (5) Use the plane wave failure method to determine the dip angle property of the work area;

[0012] (6) Design a tilt-adaptive diagonal diffusion Fourier finite difference migration method to obtain high-quality initial images;

[0013] (7) Perform multidimensional deconvolution on the point spread function of the initial image and the entire work area in the wavenumber domain to output the final offset imaging result.

[0014] Preferably, in step (2), the conventional migration imaging results Considered as the true underground reflectance After Hessian matrix The blurred image is represented as follows:

[0015] (1);

[0016] The point spread function is equivalent to the image calculated by mapping a column of the Hessian matrix to the imaging domain. Approximating the Hessian matrix using the local point spread function and spatial interpolation, followed by deblurring, is a more efficient method. The local point spread function can be expressed as:

[0017] (2);

[0018] In the formula For local point spread function, For spatial coordinates, The equation for efficiently calculating the local point spread function using the Fourier finite difference method, where the spatial perturbation distance is:

[0019] (3);

[0020] In the formula, Re represents taking the real part. For frequency band range, Angular frequency, For the source wavelet, The spatial coordinates of the receiver point, The spatial coordinates of the firing point, For the Fourier finite difference depth propagation operator, To represent the transpose, the depth extension operator of the Fourier finite difference method is specifically expressed as:

[0021] (4);

[0022] In the formula For two-dimensional wave field values, The horizontal position Vertical depth position The imaginary unit, For depth Background wave propagation speed, This represents the wave propagation speed at the corresponding spatial location. and represents the matching parameters corresponding to the order of the Fourier finite difference, where By setting scattering points at equal intervals in the initial velocity model, and taking the velocity of the point as 1.5 times that of the initial model, the scattering point velocity model is obtained. Then, the imaging value at the location of these scattering points is obtained by using equation (3), which is approximately the point spread function of that point.

[0023] Preferably, in step (3), the Gaussian filter expression applied to the local point spread function is:

[0024] (5);

[0025] In the formula, These are the two-dimensional local point spread function values ​​after Gaussian filtering. is the standard deviation of the Gaussian kernel function.

[0026] Preferably, in step (4), the formula for the inverse distance weighted interpolation method is:

[0027] (6);

[0028] The formula selects the four known points closest to the interpolation points. Perform interpolation. For these four known Spatial location, The point spread function value for the entire region is obtained after inverse distance weighted interpolation. This represents the spatial location of the initially calculated partial point spread function. The interpolation weighting coefficients are expressed as follows:

[0029] (7);

[0030] In the formula The value can be any positive real number; the larger the value, the more emphasis is placed on the role of the nearest control point.

[0031] Preferably, in step (6), the conventional Fourier finite difference migration imaging conditions are:

[0032] (8);

[0033] In the formula For the imaging results, For the total number of guns, For frequency band range, This is the complex transpose of the source wave field. To suppress the influence of low-frequency noise and inter-layer crosstalk noise on the final inversion imaging results in conventional imaging methods, a tilt-adaptive diagonal diffusion operator was designed to operate under conventional imaging conditions, representing the wavefield at the detector point.

[0034] (9);

[0035] In the formula The initial imaging result is obtained after processing with a tilt-adaptive diagonal diffusion operator. A diagonal diffusion operator that is adaptive to tilt angle. Denotes matrix convolution, where The specific expression is:

[0036] (10);

[0037] In the formula This is the tilt angle attribute. For the number of iterations, Where is the diffusion coefficient. This is the tilt modulation function for the tilt angle variable.

[0038] Preferably, in step (7), by introducing a point spread function, the blurring effect of the seismic migration imaging results is simplified to a local convolutional form:

[0039] (11);

[0040] in The radius chosen for the point spread function is half the spacing between scattering points in step (2). To simplify the calculation, the convolution is converted into a product in the wavenumber domain:

[0041] (12);

[0042] In the formula Indicates Fourier transform, Representing the wavenumber, the initial imaging result can be deblurred by performing deconvolution of the imaging result and the point spread function in two-dimensional wavenumber, and then transforming the result back to the spatial domain.

[0043] (13);

[0044] in Indicates the inverse Fourier transform. It is a positive real number used to ensure numerical stability.

[0045] The present invention has the following beneficial effects:

[0046] Compared with existing conventional Fourier finite difference migration imaging techniques, this method can effectively improve imaging resolution and deep amplitude, especially to make up for the shortcomings of the single-pass wave method in imaging regions with drastic changes in lateral velocity and steep structures. Compared with the mainstream least squares method that improves the effect of conventional Fourier finite difference migration imaging techniques, this method has higher computational efficiency and occupies less memory. Attached Figure Description

[0047] Figure 1 This is a schematic diagram of the process of the present invention;

[0048] Figure 2 The local point spread function of the Marmousi model is obtained by the scattering point perturbation method;

[0049] Figure 3 The effect of Gaussian filtering on suppressing artifacts in the point spread function that are far from the center point is shown in the following order: point spread function before Gaussian filtering suppression, and point spread function after Gaussian filtering suppression.

[0050] Figure 4 The inclination property of the work area obtained by the plane wave destruction method;

[0051] Figure 5 Smoothed velocity model used for Marmousi's accurate velocity model and offset method;

[0052] Figure 6 The performance of different imaging methods for the Marmousi model is compared as follows: conventional Fourier finite difference method, reflectivity model as a control, conventional Fourier finite difference deconvolution imaging method based on point spread function, and Fourier finite difference deconvolution imaging method based on tilt angle adaptive noise suppression.

[0053] Figure 7 for Figure 6 Comparison of two imaging details between the conventional Fourier finite difference deconvolution imaging method and the Fourier finite difference deconvolution imaging method based on tilt-adaptive noise suppression;

[0054] Figure 8 This is a velocity model of an actual land work area and the corresponding local point spread function;

[0055] Figure 9 This is the actual seismic gun record for the land-based work area;

[0056] Figure 10 The effects of different imaging methods used in this land work area are compared as follows: conventional Fourier finite difference method, least squares Fourier finite difference imaging method with 6 iterations, conventional Fourier finite difference deconvolution imaging method based on point spread function, and Fourier finite difference deconvolution imaging method based on tilt angle adaptive noise suppression.

[0057] Figure 11 for Figure 10 Two magnified images showing details from the imaging results of the four methods. Detailed Implementation

[0058] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings and specific examples:

[0059] Combination Figure 1 A Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression includes...

[0060] (1) Input velocity field, frequency domain source wavelet, and observation system parameters.

[0061] (2) Construct the velocity field of equally spaced scattering points based on the velocity field, and use the Fourier finite difference continuation method to obtain the point spread function at these scattering points.

[0062] Conventional migration imaging results Considered as the true underground reflectance After Hessian matrix The blurred image is represented as follows:

[0063] (1);

[0064] To eliminate the blurring effect of the Hessian matrix and obtain more accurate imaging results, the least squares method is usually used to iteratively optimize the initial offset imaging results. However, this process is computationally expensive. The point spread function is equivalent to the image calculated by mapping a column of the Hessian matrix to the imaging domain. Approximating the Hessian matrix using the local point spread function and spatial interpolation, and then performing deblurring, is a more efficient method. The local point spread function can be expressed as:

[0065] (2);

[0066] In the formula For local point spread function, For spatial coordinates, The equation for efficiently calculating the local point spread function using the Fourier finite difference method, where the spatial perturbation distance is:

[0067] (3);

[0068] In the formula, Re represents taking the real part. For frequency band range, Angular frequency, For the source wavelet, The spatial coordinates of the receiver point, The spatial coordinates of the firing point, For the Fourier finite difference depth propagation operator, To represent the transpose, the depth extension operator of the Fourier finite difference method is specifically expressed as:

[0069] (4);

[0070] In the formula For two-dimensional wave field values, The horizontal position Vertical depth position The imaginary unit, For depth Background wave propagation speed, This represents the wave propagation speed at the corresponding spatial location. and represents the matching parameters corresponding to the order of the Fourier finite difference, where By setting scattering points at equal intervals in the initial velocity model, and taking the velocity of the point as 1.5 times that of the initial model, the scattering point velocity model is obtained. Then, the imaging value at the location of these scattering points is obtained by using equation (3), which is approximately the point spread function of that point.

[0071] (3) Apply Gaussian filtering to the obtained partial point spread function to suppress artifacts far from the center point.

[0072] The Gaussian filtering expression applied to the local point spread function is:

[0073] (5);

[0074] In the formula, These are the two-dimensional local point spread function values ​​after Gaussian filtering. is the standard deviation of the Gaussian kernel function; the value used in the numerical examples is 3.

[0075] (4) The point spread function of the entire work area is obtained by applying the inverse distance weighted interpolation method.

[0076] The formula for the inverse distance weighted interpolation method is:

[0077] (6);

[0078] The formula selects the four known points closest to the interpolation points. Perform interpolation. For these four known Spatial location, The point spread function value for the entire region is obtained after inverse distance weighted interpolation. This represents the spatial location of the initially calculated partial point spread function. The weighting coefficients are expressed as follows:

[0079] (7);

[0080] In the formula The value can be any positive real number. The larger the value, the more it emphasizes the role of the nearest control point. The value used in the example is 2.

[0081] (5) Use the plane wave destruction method to determine the inclination properties of the work area.

[0082] The plane wave destruction method is a highly efficient technique for extracting seismic dip attributes. Its basic principle is based on the assumption that local seismic data can be composed of a series of superimposed plane waves with different dip angles. By designing a predictive filter to eliminate these in-phase axes, the dip angle corresponding to the minimum residual energy is the local true dip angle. Taking two-dimensional data as an example, the main steps include: 1) extracting data at the target point with a time window; 2) constructing a plane wave predictive filter containing dip angle parameters; 3) scanning different dip angles and solving for the optimal dip angle value that minimizes the prediction error; 4) outputting the dip angle attribute at that point, and calculating the dip angle data volume for the entire seismic area point by point.

[0083] (6) Design a tilt-adaptive diagonal diffusion Fourier finite difference migration method to obtain high-quality initial images.

[0084] Conventional Fourier finite difference migration imaging conditions:

[0085] (8);

[0086] In the formula For the imaging results, The total number of guns, For frequency band range, This is the complex transpose of the source wave field. To suppress the influence of low-frequency noise and inter-layer crosstalk noise on the final inversion imaging results in conventional imaging methods, a tilt-adaptive diagonal diffusion operator was designed to operate under conventional imaging conditions, representing the wavefield at the detector point.

[0087] (9);

[0088] In the formula The initial imaging result is obtained after processing with a tilt-adaptive diagonal diffusion operator. A diagonal diffusion operator that is adaptive to tilt angle. Denotes matrix convolution, where The specific expression is:

[0089] (10);

[0090] In the formula This is the tilt angle attribute. For the number of iterations, Where is the diffusion coefficient. This is the tilt modulation function for the tilt angle variable.

[0091] (7) Perform multidimensional deconvolution on the point spread function of the initial image and the entire work area in the wavenumber domain to output the final offset imaging result.

[0092] By introducing a point spread function, the blurring effect of seismic migration imaging results is simplified to a local convolutional form:

[0093] (11);

[0094] in The radius chosen for the point spread function is half the spacing between scattering points in step (2). To simplify the calculation, the convolution is converted into a product in the wavenumber domain:

[0095] (12);

[0096] In the formula Indicates Fourier transform, Representing the wavenumber, the initial imaging result can be deblurred by performing deconvolution of the imaging result and the point spread function in two-dimensional wavenumber, and then transforming the result back to the spatial domain.

[0097] (13);

[0098] in Indicates the inverse Fourier transform. It is a positive real number used to ensure numerical stability.

[0099] (8) Output the final imaging results.

[0100] (a) Analysis of the test results of the Marmousi model

[0101] This invention presents a Fourier finite difference deconvolution imaging method based on tilt-adaptive noise suppression, which has been applied to the Marmousi model and actual land data, achieving ideal computational results. Figure 1 A flowchart illustrating the present invention is shown. Figure 2 This paper presents partial point spread functions (PSFs) of the Marmousi model obtained through the scattering point perturbation method. Magnified images at three locations reveal different shape characteristics of the PSFs at various points within the work area. In the shallow region, with uniform illumination and simple geological structure, the obtained PSFs exhibit good energy focusing and high resolution. In the complex central region, due to drastic velocity variations and complex geological structures, the PSFs show severe deformation, with dramatic changes in amplitude, energy, and structure between adjacent PSFs. In the deep region, due to complex propagation paths and offset limitations, the obtained PSFs exhibit weak energy focusing, resulting in a significant decrease in overall resolution. Figure 3The effect of Gaussian filtering on the artifacts of the Marmousi model point spread function far from the center point is shown. It can be seen that the energy of each point spread function far from the center point is significantly suppressed. Figure 4 The dip angle properties of the Marmousi model work area obtained by the plane wave destruction method. Figure 5 The smoothed velocity model used for the Marmousi accurate velocity model and offset method.

[0102] Figure 6 The first and second images in the diagram show the results of conventional Fourier finite-difference migration imaging using the Marmousi model and the true reflectance model, respectively. Comparing the two images reveals that while Fourier finite-difference migration can accurately locate flat-surface regions, it suffers from low imaging axis amplitude energy and insufficient continuity in deep regions and areas with drastic velocity changes in the middle. Traditional deconvolution imaging methods are characterized by enhanced noise. To verify the difference in the effectiveness of deconvolution between conventional Fourier finite-difference imaging and Fourier finite-difference imaging based on tilt-adaptive noise suppression, this paper examines the results of conventional Fourier finite-difference imaging and Fourier finite-difference imaging based on tilt-adaptive noise suppression. Figure 6 The third and fourth images in the diagram show the imaging results of the conventional Fourier finite difference deconvolution imaging method and the method of this patent, respectively. Figure 6 Compared to the first image, the two images obtained by introducing PSF deblurring show significantly better characterization of the structure. The continuity of the imaging axis in the complex central region is significantly improved, and the overall amplitude is also enhanced. However, we can see... Figure 6 The third image reveals a lot of crosstalk noise and low-frequency noise enhanced by the point spread function, which greatly reduces the final quality of deconvolution imaging. Figure 6 The fourth figure shows the imaging effect obtained by using the Fourier finite difference method with tilt-adaptive noise suppression as input. It can be seen that the method of this patent can not only preserve the optimization of Fourier finite difference with migration imaging results by point spread function deconvolution to the greatest extent, but also effectively suppress the interference of noise components. Figure 7 Showing Figure 6 The third picture and Figure 6 The magnified view of two details in the fourth image more clearly shows the improvement of the inversion imaging results by suppressing the noise of the initial migration result in advance.

[0103] (b) Analysis of the effectiveness of actual land data testing

[0104] Figure 8 This is a velocity model of an actual land-based work area and corresponding partial point spread function images. Figure 9 This is the actual seismic shot record for this land-based work area. The shot record contains common shot gathers for 264 shots, each received by 240 receivers at a sampling interval of 40m. The data underwent preprocessing including bad path removal and direct wave static correction.

[0105] Figure 10 The effects of different imaging methods on this land work area are compared in the following order: conventional Fourier finite difference method, least squares Fourier finite difference imaging method with 6 iterations, conventional Fourier finite difference deconvolution imaging method based on point spread function, and Fourier finite difference deconvolution imaging method based on tilt angle adaptive noise suppression. Figure 10 The first image shows that Fourier finite-difference migration can accurately image the overall structural contours of the actual data area, clearly depicting the shallow sedimentary layers and the right-side flat-bed region. However, similar to the conclusions of the model tests mentioned earlier, the imaging results obtained using only the Fourier finite-difference method have relatively low overall resolution, and the imaging energy of steep structural boundaries is weak. Figure 10 A comparison shows that both the data domain least squares migration algorithm and the deconvolution method can effectively improve the resolution and deep amplitude energy of the imaging results, and can make up for the shortcomings of the conventional Fourier finite difference one-way wave algorithm in imaging high and steep structures. Figure 11 for Figure 10 Two magnified details of the imaging results show that in the deconvolution imaging results obtained through Fourier finite difference (FTD), least squares Fourier finite difference (LSF) migration, and tilt-free adaptive noise suppression, interlayer crosstalk noise with multiple channels intersecting the effective axis severely affects the accuracy of the imaging results. While LSF and point spread function deconvolution methods improve the initial imaging effect, they also further enhance the relative intensity of noise; the noise interference in their images is more severe than that of traditional Fourier finite difference migration. Compared with the previous three methods, the tilt-adaptive noise suppression-based Fourier finite difference deconvolution imaging results proposed in this patent effectively suppress noise, resulting in more accurate overall structural information and significantly better imaging performance than the other three methods.

[0106] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression, characterized in that, include (1) Input velocity field, frequency domain source wavelet, and observation system parameters; (2) Construct the velocity field of equally spaced scattering points based on the velocity field, and use the Fourier finite difference continuation method to obtain the point spread function at these scattering points; (3) Apply Gaussian filtering to the obtained partial point spread function to suppress artifacts far from the center point; (4) The point spread function of the entire work area is obtained by applying the inverse distance weighted interpolation method; (5) Use the plane wave failure method to determine the dip angle property of the work area; (6) Design a tilt-adaptive diagonal diffusion Fourier finite difference migration method to obtain high-quality initial images; (7) Perform multidimensional deconvolution on the point spread function of the initial image and the entire work area in the wavenumber domain to output the final offset imaging result.

2. The Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression as described in claim 1, characterized in that, In step (2), the conventional migration imaging results Considered as the true underground reflectance After Hessian matrix The blurred image is represented as follows: （1）； The point spread function is equivalent to the image calculated by mapping a column of the Hessian matrix to the imaging domain. Approximating the Hessian matrix using the local point spread function and spatial interpolation, followed by deblurring, is a more efficient method. The local point spread function can be expressed as: （2）； In the formula For local point spread function, For spatial coordinates, The equation for efficiently calculating the local point spread function using the Fourier finite difference method, where the spatial perturbation distance is: （3）； In the formula, Re represents taking the real part. For frequency band range, Angular frequency, For the source wavelet, The spatial coordinates of the receiver point, The spatial coordinates of the firing point, For the Fourier finite difference depth propagation operator, To represent the transpose, the depth extension operator of the Fourier finite difference method is specifically expressed as: （4）； In the formula For two-dimensional wave field values, The horizontal position Vertical depth position The imaginary unit, For depth Background wave propagation speed, This represents the wave propagation speed at the corresponding spatial location. and represents the matching parameters corresponding to the order of the Fourier finite difference, where By setting scattering points at equal intervals in the initial velocity model, and taking the velocity of the point as 1.5 times that of the initial model, the scattering point velocity model is obtained. Then, the imaging value at the location of these scattering points is obtained by using equation (3), which is approximately the point spread function of that point.

3. The Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression as described in claim 1, characterized in that, In step (3), the Gaussian filter expression applied to the local point spread function is: （5）； In the formula, These are the two-dimensional local point spread function values ​​after Gaussian filtering. is the standard deviation of the Gaussian kernel function.

4. The Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression as described in claim 1, characterized in that, In step (4), the formula for the inverse distance weighted interpolation method is: （6）； The formula selects the four known points closest to the interpolation points. Perform interpolation. For these four known Spatial location, The point spread function value for the entire region is obtained after inverse distance weighted interpolation. This represents the spatial location of the initially calculated partial point spread function. The interpolation weighting coefficients are expressed as follows: （7）； In the formula The value can be any positive real number; the larger the value, the more emphasis is placed on the role of the nearest control point.

5. The Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression as described in claim 1, characterized in that, In step (6), the conventional Fourier finite difference migration imaging conditions are as follows: （8）； In the formula For the imaging results, The total number of guns, For frequency band range, This is the complex transpose of the source wave field. To suppress the influence of low-frequency noise and inter-layer crosstalk noise on the final inversion imaging results in conventional imaging methods, a tilt-adaptive diagonal diffusion operator was designed to operate under conventional imaging conditions, representing the wavefield at the detector point. （9）； In the formula The initial imaging result is obtained after processing with a tilt-adaptive diagonal diffusion operator. A diagonal diffusion operator that is adaptive to tilt angle. Denotes matrix convolution, where The specific expression is: （10）； In the formula This is the tilt angle attribute. For the number of iterations, Where is the diffusion coefficient. This is the tilt modulation function for the tilt angle variable.

6. The Fourier finite-difference deconvolution imaging method based on tilt-adaptive noise suppression as described in claim 1, characterized in that, In step (7), by introducing a point spread function, the blurring effect of the seismic migration imaging results is simplified to a local convolutional form: （11）； in The radius chosen for the point spread function is half the spacing between scattering points in step (2). To simplify the calculation, the convolution is converted into a product in the wavenumber domain: （12）； In the formula Indicates Fourier transform, Representing the wavenumber, the initial imaging result can be deblurred by performing deconvolution of the imaging result and the point spread function in two-dimensional wavenumber, and then transforming the result back to the spatial domain. （13）； in Indicates the inverse Fourier transform. It is a positive real number used to ensure numerical stability.

Images