Fast spatial convolution method for point spread function

Abstract

The application provides a point spread function fast spatial convolution method, which comprises the following steps: step 1, calculating an imaging profile and a point spread function field; step 2, determining the spatial position coordinates of the center position of the point spread function by circulation; step 3, determining an interpolation method and calculating the interpolation coefficients of the local imaging profile which contributes to the point spread function; step 4, performing wave number domain fast spatial convolution; step 5, multiplying the profile by the corresponding interpolation coefficients after filtering; and step 6, accumulating the filtering results of the local imaging profile to obtain a final filtered profile. The point spread function fast spatial convolution method can greatly improve the calculation efficiency of the point spread function spatial convolution filtering process and reduce the calculation cost based on the linear accumulation characteristics of the interpolation method.

Description

Technical Field

[0001] This invention relates to the field of seismic data processing technology for petroleum geophysical exploration, and in particular to a fast spatial convolution method for point spread functions. Background Technology

[0002] As a key technology in oil and gas exploration, migration imaging plays a vital role in both theoretical research and practical applications of seismic exploration, and is increasingly important in imaging increasingly complex and refined subsurface media. However, traditional migration imaging methods are merely the conjugate transpose of the Born forward operator, which can only produce fuzzy structural imaging results for band-limited seismic data acquired by finite observation systems. Furthermore, complex subsurface structures and irregular spatial sampling of seismic data can lead to migration artifacts and non-uniform imaging illumination, severely affecting the reliability of imaging amplitude.

[0003] Least square migration based on linearized seismic inversion theory is an effective approach to solving the aforementioned problems. This technique treats seismic migration imaging as a linear inverse problem and uses optimization methods for iterative solutions. Theoretically, it can eliminate the adverse effects of irregular acquisition and band-limited wavelet on imaging results, improving imaging resolution and illumination. However, the theoretical advantages of this technique have not translated into practical applications, primarily due to its prohibitively high computational cost. Classical data-domain least square migration requires a conventional forward modeling and migration operation for each iteration, and the complete inversion process often requires tens of times the computational load of conventional migration. Therefore, the enormous computational cost severely restricts the practicality of this method. Imaging-domain least square migration requires the computation and storage of a massive Hessian matrix. Even with optimization techniques such as data encoding, the computational and storage overhead remains expensive, making it difficult to apply to large-scale 3D real-world seismic data processing.

[0004] Chinese patent application CN202010476294.0 discloses a depth-domain imaging simulation method and system. The method includes: dividing an initial density model into windows and setting density perturbation points within each window; determining a density model at the perturbation points based on the initial density model under the influence of these perturbation points; obtaining point scattering wavefield data based on the perturbation point density model and the initial density model; performing migration processing on the point scattering wavefield data to obtain a point spread function; and obtaining a depth-domain simulated image based on the point spread function and a reflectivity model. This invention can obtain a stable and continuous simulated image that matches the actual migration imaging.

[0005] Chinese patent application CN201210290287.7 discloses a point spread function (PSF) estimation method based on a Kalman model for image restoration. This invention first establishes a symmetric full-plane Kalman state equation and observation equation for PSF estimation; secondly, it determines the observation matrix of the PSF model based on selected stationary regions of the image; finally, it performs Kalman iteration based on the stationary regions to estimate the PSF of the image. This invention effectively improves upon the shortcomings of inaccurate and non-adaptive PSF estimation, enabling adaptive PSF estimation based on the image, resulting in more accurate PSF values. This provides a more accurate PSF for subsequent image restoration.

[0006] Chinese patent application CN201910427576.9 discloses a point spread function reconstruction method, comprising: (1) inputting a target point spread function; (2) sampling the input target point spread function; (3) estimating the initial phase and initial amplitude of each sampled point as a dipole; (4) calculating the radiation field of the dipole on the pupil surface; (5) obtaining the amplitude and phase distribution on the pupil surface through the radiation field; (6) generating a point spread function through the obtained amplitude and phase; (7) comparing the point spread function with the target point spread function, and if the requirements are met, outputting the obtained amplitude and phase distribution; if the requirements are not met, updating the initial phase and initial amplitude, and returning to step (4). This invention directly obtains the phase and amplitude distribution on the pupil surface from the target point spread function, greatly improving computational efficiency. Furthermore, since the light radiated by the dipole propagates in random directions, this method can be extended to multiple microscopic systems with multiple pupil surfaces.

[0007] The existing technologies described above are significantly different from the present invention and have failed to solve the technical problem we want to address. Therefore, we have invented a new method for fast spatial convolution of point spread functions. Summary of the Invention

[0008] The purpose of this invention is to provide a fast spatial convolution method for point spread functions that improves the computational efficiency of least-squares offset in the imaging domain.

[0009] The objective of this invention can be achieved through the following technical measures: a fast spatial convolution method for point spread functions, which includes:

[0010] Step 1: Calculate the imaging profile and point spread function field;

[0011] Step 2: Iterate through the point spread function to determine the spatial coordinates of its center position;

[0012] Step 3: Determine the interpolation method and calculate the interpolation coefficients for the local imaging profiles that contribute to the point spread function;

[0013] Step 4: Perform fast spatial convolution in the wavenumber domain;

[0014] Step 5: After filtering, multiply the profile by the corresponding interpolation coefficients;

[0015] Step 6: Accumulate the filtering results of the local imaging profiles to obtain the final filtered profile.

[0016] The objective of this invention can also be achieved through the following technical measures:

[0017] In step 1, the migration imaging profile and the point spread function field located in the sparse space are calculated using seismic migration and point spread function.

[0018] In step 2, the calculation loop begins, the spread function of each point is read, and the spatial location of its center point is determined.

[0019] In step 3, the local imaging range in which the spread function contributes to the point is determined according to the selected interpolation method, as well as the interpolation coefficients for each imaging sample point.

[0020] In step 4, the product of the point spread function and the wavenumber spectrum of the local imaging profile is calculated using the fast Fourier transform.

[0021] In step 4, the formula for calculating the product of the point spread function and the wavenumber spectrum of the local imaging profile is:

[0022] M(k)=H(k)R(k)

[0023] Where H(k) and R(k) are the wavenumber spectra of the point spread function and the local imaging profile, respectively, and M(k) is the wavenumber spectrum of the local imaging profile after filtering with the point spread function.

[0024] In step 5, the wavenumber spectrum is inversely transformed to the spatial domain to obtain the spatially convolutionally filtered imaging profile, which is then multiplied by the corresponding interpolation coefficients.

[0025] In step 5, the formula for multiplying the filtered profile by the corresponding interpolation coefficients is as follows:

[0026] m(x) = cof(x) * iFFT(M(k))

[0027] Where m(x) is the spatial domain filtered imaging profile after applying interpolation coefficients, cof(x) is the interpolation coefficient at spatial location x, and iFFT represents the inverse Fourier transform.

[0028] In step 6, repeat steps (2) to (5) above, and accumulate the image profiles after convolution filtering to obtain the final point spread function convolution filtering result.

[0029] In step 6, the formula for calculating the point spread function spatial filtering imaging profile is:

[0030] m′(x)=m′(x)+m(x)

[0031] Where m′(x) is the final point spread function spatial filtering imaging profile; m(x) is the spatial domain filtering imaging profile after applying interpolation coefficients.

[0032] The fast spatial convolution method for point spread function (PSF) in this invention first uses FFT for fast spatial convolution filtering, then applies appropriate interpolation coefficients and sums the filtering results to obtain the same computational result as spatial convolution filtering. This invention provides a fast PSF spatial convolution filtering method that, based on the linear accumulation characteristic of interpolation methods, uses FFT to convolve the PSF with a local imaging profile, and then applies interpolation coefficients to calculate and sum the results to obtain the convolution filtering result. Compared to traditional methods, this method can significantly improve the computational efficiency of the PSF spatial convolution filtering process while maintaining computational accuracy and reducing computational costs. Attached Figure Description

[0033] Figure 1 This is a flowchart of a specific embodiment of the point spread function fast spatial convolution method of the present invention;

[0034] Figure 2 This is a schematic diagram of the reflectance profile corresponding to the Marmousi model in a specific embodiment of the present invention;

[0035] Figure 3 This is a schematic diagram of the point spread function field located in the sparse space of the Marmousi model in a specific embodiment of the present invention;

[0036] Figure 4 This is a schematic diagram of the result of spatial convolution filtering using a traditional point spread function in a specific embodiment of the present invention;

[0037] Figure 5 This is a schematic diagram of the result of spatial convolution filtering using the fast point spread function given in this invention in a specific embodiment of the invention. Detailed Implementation

[0038] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0039] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0040] The point spread function fast spatial convolution method of the present invention includes the following steps:

[0041] (1) The migration imaging profile and the point spread function field located in the sparse space were obtained using the seismic migration and point spread function calculation program;

[0042] (2) Iterate through the point spread function to determine the spatial coordinates of its center position;

[0043] (3) Based on the selected interpolation method, determine the local imaging range in which the spread function contributes to the point, and the interpolation coefficients corresponding to each imaging sample point;

[0044] (4) Calculate the product of the point spread function and the wavenumber spectrum of the local imaging profile using the Fast Fourier Transform (FFT);

[0045] (5) Transform the wavenumber spectrum inversely to the spatial domain to obtain the spatially convolutionally filtered imaging profile, and multiply it by the corresponding interpolation coefficients;

[0046] (6) Repeat steps (2) to (5) above, and accumulate the image profiles after convolution filtering to obtain the final point spread function convolution filtering result.

[0047] The following are several specific embodiments of the application of the present invention.

[0048] Example 1

[0049] In a specific embodiment 1 of the present invention, Figure 1 This is a flowchart of an embodiment of the point spread function spatial convolution filtering method described in this invention. The specific implementation steps include:

[0050] (1) Calculate the imaging profile and point spread function field

[0051] The migration imaging profile and the point spread function field located in sparse spatial locations were obtained using a seismic migration and point spread function calculation program.

[0052] (2) Start the calculation loop, read the spread function of each point, and determine the spatial location of its center point.

[0053] By iterating through the point spread function, the spatial coordinates of its center position are determined.

[0054] (3) Determine the interpolation method and calculate the interpolation coefficients for the local imaging profiles that contribute to the point spread function.

[0055] Based on the selected interpolation method, determine the local imaging range in which the spread function contributes to the point, and the corresponding interpolation coefficients for each imaging sample point;

[0056] (4) Perform fast spatial convolution in the wavenumber domain

[0057] The product of the point spread function and the wavenumber spectrum of the local imaging profile is calculated using the Fast Fourier Transform (FFT).

[0058] (5) Multiply the filtered profile by the corresponding interpolation coefficient.

[0059] The wavenumber spectrum is inversely transformed to the spatial domain to obtain the spatially convolutionally filtered imaging profile, which is then multiplied by the corresponding interpolation coefficients.

[0060] (6) Accumulate the filtering results of the local imaging profiles to obtain the final filtered profile.

[0061] Repeat steps (2) to (5) above, and accumulate the image profiles after convolution filtering to obtain the final point spread function convolution filtering result.

[0062] Example 2:

[0063] In a specific embodiment 2 of the present invention, the invention is applied. Figures 2-5 The application effects of the point spread function spatial convolution filtering method described in this invention are as follows:

[0064] 1) The reflectivity profile was calculated using the Marmousi model velocity field, and the results are as follows: Figure 2 As shown;

[0065] 2) Using the point spread function calculation method, the point spread function of the Marmousi model located in the sparse space is calculated, and the results are as follows: Figure 3 As shown;

[0066] 3) Using the conventional point spread function spatial convolution method, the convolution filter profile of the point spread function and the Marmousi model reflectivity profile is obtained. The result is as follows: Figure 4 As shown, it took approximately 3.6 seconds;

[0067] 4) Using this invention, the convolutional filter profile of the point spread function and the Marmousi model reflectivity is obtained, and the results are as follows: Figure 5 As shown, it takes approximately 0.5 seconds;

[0068] Comparative analysis revealed that this invention can significantly improve the computational efficiency of the point spread function spatial convolution filtering process, reducing the computational cost by more than 7 times compared to traditional methods.

[0069] Example 3:

[0070] In a specific embodiment 3 of the present invention, the point spread function spatial convolution filtering method described in the present invention specifically includes the following steps:

[0071] (1) Calculate the imaging profile and point spread function field

[0072] The migration imaging profile and the point spread function field located in sparse spatial locations are obtained using a seismic migration and point spread function calculation program; the calculation of the imaging profile and point spread function can be implemented using conventional algorithms.

[0073] (2) Start the calculation loop, read the spread function of each point, and determine the spatial location of its center point.

[0074] By iterating through the point spread function, the spatial coordinates of its center position are determined.

[0075] (3) Determine the interpolation method and calculate the interpolation coefficients for the local imaging profiles that contribute to the point spread function.

[0076] Based on the selected interpolation method, determine the local imaging range in which the spread function contributes to the point, and the corresponding interpolation coefficients for each imaging sample point; the selection of the interpolation method depends on the required accuracy, but generally conventional linear interpolation can meet the accuracy requirements.

[0077] (4) Perform fast spatial convolution in the wavenumber domain

[0078] The product of the point spread function and the wavenumber spectrum of the local imaging profile is calculated using the Fast Fourier Transform (FFT).

[0079] M(k)=H(k)R(k)

[0080] Where H(k) and R(k) are the wavenumber spectra of the point spread function and the local imaging profile, respectively, and M(k) is the wavenumber spectrum of the local imaging profile after filtering with the point spread function.

[0081] (5) Multiply the filtered profile by the corresponding interpolation coefficient.

[0082] The wavenumber spectrum is inversely transformed to the spatial domain to obtain the spatially convolutionally filtered imaging profile, which is then multiplied by the corresponding interpolation coefficients.

[0083] m(x) = cof(x) * iFFT(M(k))

[0084] Where m(x) is the spatial domain filtered imaging profile after applying interpolation coefficients, cof(x) is the interpolation coefficient at spatial location x, and iFFT represents the inverse Fourier transform.

[0085] (6) Accumulate the filtering results of the local imaging profiles to obtain the final filtered profile.

[0086] Repeat steps (2) to (5) above, and accumulate the image profiles after convolution filtering to obtain the final point spread function convolution filtering result.

[0087] m′(x)=m′(x)+m(x)

[0088] Where m′(x) is the final point spread function spatial filtering imaging profile.

[0089] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0090] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. A fast spatial convolution method using point spread functions, characterized in that, The fast spatial convolution method for the point spread function includes: Step 1: Calculate the imaging profile and point spread function field; Step 2: Iterate through the point spread function to determine the spatial coordinates of its center position; Step 3: Determine the interpolation method and calculate the interpolation coefficients for the local imaging profiles that contribute to the point spread function; Step 4: Perform fast spatial convolution in the wavenumber domain; Step 5: After filtering, multiply the profile by the corresponding interpolation coefficients; Step 6: Accumulate the filtering results of the local imaging profiles to obtain the final filtered profile; In step 4, the product of the point spread function and the wavenumber spectrum of the local imaging profile is calculated using the fast Fourier transform. In step 4, the formula for calculating the product of the point spread function and the wavenumber spectrum of the local imaging profile is: in, , These are the wavenumber spectra of the point spread function and the local imaging profile, respectively. The wavenumber spectrum of the local imaging profile after point spread function filtering; In step 5, the wavenumber spectrum is inversely transformed to the spatial domain to obtain the spatially convolutionally filtered imaging profile, which is then multiplied by the corresponding interpolation coefficients. In step 5, the formula for multiplying the filtered profile by the corresponding interpolation coefficients is as follows: in, It is the spatial domain filtered imaging profile after applying interpolation coefficients. For spatial location Interpolation coefficients at the location, Represents the inverse Fourier transform; In step 6, repeat steps 2-5 above and accumulate the image profiles after convolution filtering to obtain the final point spread function convolution filtering result.

2. The point spread function fast spatial convolution method according to claim 1, characterized in that, In step 1, the migration imaging profile and the point spread function field located in the sparse space are calculated using seismic migration and point spread function.

3. The point spread function fast spatial convolution method according to claim 1, characterized in that, In step 2, the calculation loop begins, the spread function of each point is read, and the spatial location of its center point is determined.

4. The point spread function fast spatial convolution method according to claim 1, characterized in that, In step 3, the local imaging range in which the spread function contributes to the point is determined according to the selected interpolation method, as well as the interpolation coefficients for each imaging sample point.

5. The fast spatial convolution method for point spread function according to claim 1, characterized in that, In step 6, the formula for calculating the point spread function spatial filtering imaging profile is: in, This is the final point spread function spatial filtering imaging profile; It is the spatial domain filtered imaging profile after applying interpolation coefficients.

Images