Diffraction wave imaging method and device

Abstract

The present invention discloses a diffraction wave imaging method and device, wherein the method comprises: calculating the imaging value and illumination times of each seismic trace at each imaging point underground according to the seismic recording data, shot point position and receiving point position of each seismic trace in the seismic data; calculating the depth domain angle trace set corresponding to each seismic trace according to the imaging value and illumination times of each seismic trace at each imaging point underground; performing sparse constrained Radon transform on multiple depth domain angle trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle trace set into three-dimensional model data; superimposing multiple Radon domain three-dimensional models to obtain diffraction wave imaging data. The present invention can improve the resolution and fidelity of diffraction wave imaging.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data processing, and in particular to a diffraction wave imaging method and device. Background Art

[0002] This section is intended to provide a background or context to the embodiments of the invention recited in the claims. No admission is made that the description herein is prior art by inclusion in this section.

[0003] The primary task of seismic exploration is to improve the ability of seismic data to distinguish underground structures and lithology. The main object of reflection wave seismic exploration is the underground large set (large-scale) layered lateral distribution of strata, which plays an irreplaceable and important role in the conventional structural oil and gas exploration stage. With the expansion of the field of oil and gas exploration, the exploration has shifted from structural exploration to lithological exploration, from medium and shallow exploration to medium and deep exploration, from sedimentary oil and gas reservoir exploration to carbonate and igneous oil and gas reservoir exploration, and even tight oil and gas reservoir exploration. This places higher and higher requirements on the accuracy of seismic exploration, and improving the resolution of seismic imaging has become a major demand for oil and gas exploration. For non-structural oil and gas reservoirs such as carbonate rocks, igneous rocks, small fault blocks and lithology, their common characteristics are small-scale inhomogeneous bodies or irregular surfaces, which are difficult to form effective reflection waves, and generate a large number of diffraction waves or scattered waves. Diffraction wave imaging technology has been a hot topic of research in the past decade and is expected to improve the imaging resolution of seismic data.

[0004] At present, diffraction wave imaging technology is mainly realized by imaging first and then separating, that is, diffraction wave separation is performed on the imaged depth domain gathers. Especially in the depth domain angular gathers, the reflection wave presents a smile curve, while the diffraction wave presents a linear feature, which provides a theoretical basis for diffraction wave imaging in complex structural areas.

[0005] When separating diffraction waves on the depth domain angle gather, the two-dimensional Radon transform is often used, but this transform cannot reflect the formation dip information, and thus cannot determine the vertex position of the reflection wave hyperbola, which affects the diffraction wave separation effect and reduces the resolution of diffraction wave imaging. If three-dimensional Radon transform is used, especially three-dimensional high-resolution Radon transform, although the quality of diffraction wave imaging is improved, when using the split Bregman iterative algorithm, wavelet deconvolution must be performed first. Due to the uncertainty of the wavelet, the uncertainty of the diffraction wave imaging is increased, and the fidelity of the diffraction wave imaging is reduced. Summary of the invention

[0006] An embodiment of the present invention provides a diffraction wave imaging method for improving the resolution and fidelity of diffraction wave imaging. The method includes:

[0007] According to the seismic record data, shot point position and receiving point position of each seismic channel in the seismic data, the imaging value and illumination times of each imaging point of each seismic channel in the underground are calculated;

[0008] According to the imaging value and illumination times of each seismic trace at each imaging point underground, the depth domain angle gather corresponding to each seismic trace is calculated;

[0009] Performing sparse constrained Radon transform on multiple depth domain angle gathers to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle gathers into three-dimensional model data;

[0010] Multiple Ladon domain three-dimensional models are superimposed to obtain diffraction wave imaging data.

[0011] The embodiment of the present invention further provides a diffraction wave imaging device for improving the resolution and fidelity of diffraction wave imaging, the device comprising:

[0012] An imaging value and illumination times calculation module is used to calculate the imaging value and illumination times of each imaging point of each seismic channel in the underground according to the seismic record data, the shot point position and the receiving point position of each seismic channel in the seismic data;

[0013] A depth domain angle gather calculation module is used to calculate the depth domain angle gather corresponding to each seismic trace according to the imaging value and illumination times of each imaging point of each seismic trace in the underground;

[0014] A sparse constrained Radon transform module is used to perform sparse constrained Radon transform on multiple depth domain angle gathers to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle gathers into three-dimensional model data;

[0015] The model superposition module is used to superimpose multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data.

[0016] An embodiment of the present invention further provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the above-mentioned diffraction wave imaging method when executing the computer program.

[0017] An embodiment of the present invention further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program for executing the above-mentioned diffraction wave imaging method.

[0018] In the embodiment of the present invention, the imaging value and illumination times of each seismic trace at each imaging point underground are calculated according to the seismic recording data, the shot point position and the receiving point position of each seismic trace in the seismic data; the depth domain angle trace set corresponding to each seismic trace is calculated according to the imaging value and illumination times of each seismic trace at each imaging point underground; the sparse constrained Radon transform is performed on multiple depth domain angle trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle trace set into three-dimensional model data; the multiple Radon domain three-dimensional models are performed The diffraction wave imaging data is obtained by superposition, so that the linear characteristics of the diffraction wave can be concentrated by calculating the imaging value and the number of illuminations, which can avoid the problem that the formation dip information cannot be retained in the prior art, thereby affecting the diffraction wave separation effect, and improving the resolution of the diffraction wave imaging; the waveform characteristics of the diffraction wave are effectively maintained through the sparse constrained Radon transform, and since the wavelet deconvolution processing is no longer needed, the problem that the uncertainty of the diffraction wave imaging is increased due to the uncertainty of the wavelet in the prior art, thereby reducing the fidelity of the diffraction wave imaging is avoided, and the fidelity of the diffraction wave imaging is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the prior art descriptions. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work. In the drawings:

[0020] Figure 1 is an example flow chart of a diffraction wave imaging method in an embodiment of the present invention;

[0021] Figure 2 A schematic diagram of the process of generating a depth domain angle gather in an embodiment of the present invention;

[0022] Figure 3 Schematic diagram of the process of sparse constrained Radon transform in an embodiment of the present invention;

[0023] Figure 4 Schematic diagram of the process of diffraction wave imaging in an embodiment of the present invention;

[0024] Figure 5 A schematic diagram of a depth domain angle gather map generated in an embodiment of the present invention;

[0025] Figure 6 It is a schematic diagram of a depth domain angle gather map after Radon forward and inverse transformation in an embodiment of the present invention;

[0026] Figure 7 Schematic diagram of a separated diffraction wave depth domain angular gather diagram in an embodiment of the present invention;

[0027] Figure 8 is a comparison diagram of the reflected wave and diffraction wave imaging and the diffraction wave imaging in an embodiment of the present invention;

[0028] Fig. 9 is the conventional Radon transform zp gather diagram and Comparison chart of the gather graph;

[0029] Fig.10 is the sparse constrained Radon transform zp gather map and Comparison chart of the gather graph;

[0030] Fig.11 Schematic diagram of a diffraction wave imaging method in an embodiment of the present invention;

[0031] Fig.12 Schematic diagram of the structure of a diffraction wave imaging device in an embodiment of the present invention. DETAILED DESCRIPTION

[0032] To make the purpose, technical solution and advantages of the embodiments of the present invention more clear, the embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings. Here, the exemplary embodiments of the present invention and their descriptions are used to explain the present invention, but are not intended to limit the present invention.

[0033] Diffraction wave imaging technology has roughly gone through three stages of development. The focus of research in the first decade of this century was separation before imaging, that is, first separating the time domain diffraction wave field from the recorded seismic data, and then performing migration imaging. This method works well in simple structural areas, but it is difficult to identify the diffraction wave field in complex structural areas due to the complexity of the wave field. The second stage is scattering target imaging, which does not require separation of the diffraction wave field first, but applies an aperture control factor during the imaging process. This factor is complementary to the reflection wave imaging aperture control factor. Since the optimal migration aperture of reflection wave imaging itself is difficult to determine in complex structural areas, scattering target imaging technology is also difficult to adapt to complex structural areas. The third stage is imaging first and then separation, that is, performing diffraction wave separation on the depth domain gather of imaging. Especially on the depth domain angle gather, the reflection wave presents a smile curve, while the diffraction wave presents a linear feature, which provides a theoretical basis for diffraction wave imaging in complex structural areas.

[0034] The inventor found that the diffraction wave imaging technology can be used to identify the positions of lithology mutations such as ancient river channels, lithology pinch-out points, faults, carbonate caves, igneous rocks, etc. However, due to the weak energy of diffraction waves, the effective application of diffraction wave imaging technology in complex tectonic areas is limited. When separating diffraction waves on depth domain angle gathers, two-dimensional Radon transform is commonly used, but this transform cannot reflect the stratum dip information, and then the vertex position of the reflection wave hyperbola cannot be determined, which affects the diffraction wave separation effect and reduces the resolution of diffraction wave imaging. If three-dimensional Radon transform is used, especially three-dimensional high-resolution Radon transform, although the diffraction wave imaging quality is improved, when using the split Bregman iterative algorithm, wavelet deconvolution must be performed first, because the uncertainty of the wavelet increases the uncertainty of the diffraction wave imaging, and the fidelity of the diffraction wave imaging is reduced.

[0035] The diffraction wave imaging method provided in the embodiment of the present invention can be used to solve the above problems and improve the resolution and fidelity of diffraction wave imaging. Fig.11 As shown, it may include:

[0036] Step 1101: Calculate the imaging value and illumination times of each imaging point in the underground of each seismic channel according to the seismic record data, shot point position and receiving point position of each seismic channel in the seismic data;

[0037] Step 1102: Calculate the depth domain angle gather corresponding to each seismic trace according to the imaging value and illumination times of each seismic trace at each imaging point underground;

[0038] Step 1103: performing sparse constrained Radon transform on multiple depth-domain angle gathers to obtain multiple Radon-domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth-domain angle gathers into three-dimensional model data;

[0039] Step 1104: Overlay multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data.

[0040] In the embodiment of the present invention, the imaging value and illumination times of each seismic trace at each imaging point underground are calculated according to the seismic recording data, the shot point position and the receiving point position of each seismic trace in the seismic data; the depth domain angle trace set corresponding to each seismic trace is calculated according to the imaging value and illumination times of each seismic trace at each imaging point underground; the sparse constrained Radon transform is performed on multiple depth domain angle trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle trace set into three-dimensional model data; the multiple Radon domain three-dimensional models are performed The diffraction wave imaging data is obtained by superposition, so that the linear characteristics of the diffraction wave can be concentrated by calculating the imaging value and the number of illuminations, which can avoid the problem that the formation dip information cannot be retained in the prior art, thereby affecting the diffraction wave separation effect, and improving the resolution of the diffraction wave imaging; the waveform characteristics of the diffraction wave are effectively maintained through the sparse constrained Radon transform, and since the wavelet deconvolution processing is no longer needed, the problem that the uncertainty of the diffraction wave imaging is increased due to the uncertainty of the wavelet in the prior art, thereby reducing the fidelity of the diffraction wave imaging is avoided, and the fidelity of the diffraction wave imaging is improved.

[0041] In specific implementation, firstly, the imaging value and illumination times of each underground imaging point of each seismic channel are calculated according to the seismic record data, shot point position and receiving point position of each seismic channel in the seismic data.

[0042] In an embodiment, the diffraction wave imaging method provided by the embodiment of the present invention may also include: calculating the illumination inclination angle of each seismic trace in the seismic data based on the travel time from the shot point and the receiving point in each seismic trace in the seismic data to each underground imaging point, the ray azimuth, the ray polar angle and the Jacobian determinant element value.

[0043] Calculating the imaging value and illumination times of each seismic channel at each underground imaging point based on the seismic record data, shot point location and receiving point location of each seismic channel in the seismic data may include: calculating the imaging value and illumination times of each seismic channel at each underground imaging point based on the illumination inclination angle of each seismic channel in the seismic data, seismic record data, shot point location, receiving point location, travel time from the shot point and the receiving point to each underground imaging point, ray azimuth, ray polar angle and Jacobian determinant element value.

[0044] In an embodiment, the illumination inclination angle of each seismic trace in the seismic data may be calculated according to the following formula:

[0045] When the lighting inclination angle is the true inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is:

[0046]

[0047] When the lighting inclination angle is the apparent inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is:

[0048]

[0049] When the lighting inclination angle is the true inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is:

[0050]

[0051] When the lighting inclination angle is the apparent inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is:

[0052]

[0053] Where, pxs = sinα SM cosβ SM , pys = sinα SM sinβ SM ,

[0054] pxg=sinα GM cosβ GM , pyg = sinα GM sinβ GM ;

[0055] pxs, pys, pxg, and pyg are all ray parameters. pxs and pys are the ray parameters of the x-direction and y-direction of the ray from the shot point S to the imaging point M; pxg and pyg are the ray parameters of the x-direction and y-direction of the ray from the receiving point R to the imaging point M; α represents the scattering angle of each seismic trace, and cos 2α=cosα SM cosα GM +sinα SM sinα GM cos(β SM -β GM ), α SM and β SM are the polar angle and azimuth angle of the ray from the shot point S to the imaging point M at the imaging point M; α GM and β GM They are the polar angle and azimuth angle of the ray from the receiving point R to the imaging point M at the imaging point M respectively.

[0056] In the above embodiment, the travel time and other information are read from the travel time and other tables, and then the next imaging can be continued through the Kirchhoff integral imaging theory, and the depth domain angle gather is output after all the tracks are imaged; the embodiment of the present invention utilizes the linear characteristics of the diffraction wave in the depth domain illumination dip angle gather, and effectively maintains the waveform characteristics of the diffraction wave through the three-dimensional sparse constrained Radon transform, thereby improving the resolution and fidelity of the diffraction wave imaging.

[0057] In specific implementation, after calculating the imaging value and illumination times of each underground imaging point of each seismic channel based on the seismic record data, shot point position and receiving point position of each seismic channel in the seismic data, the depth domain angular track set corresponding to each seismic channel is calculated based on the imaging value and illumination times of each underground imaging point of each seismic channel.

[0058] In the embodiment, the depth domain angle gather corresponding to each seismic trace is calculated according to the following formula:

[0059]

[0060] Wherein, I(ix,iy,iz;θ) represents the set of nx×ny depth domain angular trace sets corresponding to each seismic trace; img(ix,iy,iz;k;θ) represents the imaging value of the k-th seismic trace at each underground imaging point; ni(ix,iy,iz;k;θ) represents the number of illuminations of the k-th seismic trace at each underground imaging point, ix, iy and iz are the numbers of the underground imaging points on the x, y and z coordinate axes respectively; nx, ny and nz are the total number of discrete points of the underground imaging points on the x, y and z coordinate axes respectively.

[0061] In the above embodiment, by calculating the depth domain angle gathers corresponding to each seismic trace, it is helpful to perform sparse constrained Radon transform on multiple depth domain angle gathers in the subsequent steps, which can concentrate the linear characteristics of the diffraction waves, avoid the problem in the prior art that the formation dip information cannot be retained, which affects the diffraction wave separation effect, and improve the resolution of diffraction wave imaging.

[0062] In specific implementation, after calculating the depth domain angular trace set corresponding to each seismic trace based on the imaging value and illumination times of each imaging point underground of each seismic trace, sparse constrained Radon transform is performed on multiple depth domain angular trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angular trace set into three-dimensional model data.

[0063] In an embodiment, sparse constrained Radon transform is performed on multiple depth-domain angle gathers to obtain multiple Radon domain three-dimensional models, which may include: for each depth-domain angle gather, performing Fourier transform, Radon transform and inverse Fourier transform on the depth-domain angle gather to obtain an initial value of the Radon domain three-dimensional model; calculating the diagonal elements of the weight matrix of the Radon domain three-dimensional model based on the initial value of the Radon domain three-dimensional model; and performing L1 norm sparse constraint on the initial value of the Radon domain three-dimensional model based on the diagonal elements of the weight matrix of the Radon domain three-dimensional model to obtain a Radon domain three-dimensional model corresponding to the depth-domain angle gather.

[0064] In the above embodiment, the waveform characteristics of the diffraction wave are effectively maintained by the sparse constrained Radon transform. Since the wavelet deconvolution process is no longer required, the problem of increasing the uncertainty of the diffraction wave imaging due to the uncertainty of the wavelet and reducing the fidelity of the diffraction wave imaging in the prior art is avoided, thereby improving the fidelity of the diffraction wave imaging.

[0065] During specific implementation, for each depth domain angle gather, Fourier transform, Radon transform and inverse Fourier transform are performed on the depth domain angle gather to obtain the initial value of the Radon domain three-dimensional model, including: for each depth domain angle gather, Fourier transform is performed on the depth domain angle gather to obtain the depth wavenumber domain two-dimensional array of the depth domain angle gather; Radon transform and inverse Fourier transform are performed on the depth wavenumber domain two-dimensional array to obtain the initial value of the Radon domain three-dimensional model.

[0066] In the embodiment, the diagonal elements of the weight matrix of the Ladon domain three-dimensional model may be calculated according to the initial value of the Ladon domain three-dimensional model according to the following formula:

[0067]

[0068] Among them, [W m ] ii The model weight matrix W represents the three-dimensional model of the Radon domain m The i-th diagonal element of the weight matrix in is denoted by b is the mean value of m in the three-dimensional model of the Ladon domain, Can be written as m i , Indicates that the sequence numbers in the three-dimensional model m of the Ladon domain are iz, ip and The three-dimensional model value of the Ladon domain.

[0069] In the above embodiment, the diagonal elements of the weight matrix of the Radon domain three-dimensional model are calculated according to the initial values ​​of the Radon domain three-dimensional model, which is helpful for calculating the diffraction wave imaging data in the subsequent steps.

[0070] In a specific implementation, after sparse constrained Radon transform is performed on multiple depth-domain angular gathers to obtain multiple Radon-domain three-dimensional models, the multiple Radon-domain three-dimensional models are superimposed to obtain diffraction wave imaging data.

[0071] In an embodiment, superimposing multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data may include: extracting model trace data with a curvature within a preset numerical range in each Ladon domain three-dimensional model; performing Fourier transform processing on the depth coordinates of the seismic trace in each Ladon domain three-dimensional model to obtain depth wavenumber domain data of multiple Ladon domain three-dimensional models; performing Radon inverse transform calculation and depth domain Fourier inverse transform calculation on the extracted model trace data with a curvature within a preset numerical range in each Ladon domain three-dimensional model and the depth wavenumber domain data of multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data corresponding to each seismic trace; and superimposing the diffraction wave imaging data corresponding to each seismic trace to obtain the diffraction wave imaging data.

[0072] In the embodiment of the present invention, the imaging value and illumination times of each seismic trace at each imaging point underground are calculated according to the seismic recording data, the shot point position and the receiving point position of each seismic trace in the seismic data; the depth domain angle trace set corresponding to each seismic trace is calculated according to the imaging value and illumination times of each seismic trace at each imaging point underground; the sparse constrained Radon transform is performed on multiple depth domain angle trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle trace set into three-dimensional model data; the multiple Radon domain three-dimensional models are performed The diffraction wave imaging data is obtained by superposition, so that the linear characteristics of the diffraction wave can be concentrated by calculating the imaging value and the number of illuminations, which can avoid the problem that the formation dip information cannot be retained in the prior art, thereby affecting the diffraction wave separation effect, and improving the resolution of the diffraction wave imaging; the waveform characteristics of the diffraction wave are effectively maintained through the sparse constrained Radon transform, and since the wavelet deconvolution processing is no longer needed, the problem that the uncertainty of the diffraction wave imaging is increased due to the uncertainty of the wavelet in the prior art, thereby reducing the fidelity of the diffraction wave imaging is avoided, and the fidelity of the diffraction wave imaging is improved.

[0073] A specific embodiment is given below to illustrate the specific application of the method of the present invention. Figure 1 As shown, it is a general flow chart provided by an embodiment of the present invention, and the method capable of realizing three-dimensional high-resolution Radon transform diffraction wave fidelity imaging includes:

[0074] S101, depth domain angle gather generation: that is, according to the seismic record data, shot point position and receiving point position of each seismic channel in the seismic data, the imaging value and illumination times of each seismic channel at each imaging point underground are calculated; according to the imaging value and illumination times of each seismic channel at each imaging point underground, the depth domain angle gather corresponding to each seismic channel is calculated;

[0075] S102, sparse constrained Radon transform: that is, performing sparse constrained Radon transform on multiple depth domain angle gathers to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle gathers into three-dimensional model data;

[0076] S103, diffraction wave imaging: that is, superimposing multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data.

[0077] Combine the following Figure 2-Figure 10 This embodiment is described in detail:

[0078] S101, depth domain angle gather generation:

[0079] Figure 2 A flow chart of depth domain angle gather generation is shown, and the method comprises:

[0080] S201. Input a seismic record data u(i), (it=0,1,…,nt-1), where nt is the number of sampling points of the record. At the same time, input the shot point position S of the record. k (x s ,y s ,z s ) and the receiving point location G k (x g ,y g ,z g ), where k is the sequence number of the channel. In the embodiment of the present invention, the number of sampling points nt per channel is 3500 and the sampling interval is 2 milliseconds.

[0081] S202, read the gun point position S from the travel time table calculated in advance k and receiving point location G k The information such as travel time, ray azimuth, ray polar angle, Jacobian determinant elements, etc. to each underground imaging point position M(ix,iy,iz) (ix=0,1,…,nx-1;iy=0,1,…,ny-1;iz=0,1,…,nz-1), wherein nx,ny,nz are the total number of discrete points in the x,y,z directions of the three-dimensional imaging area, respectively. The grid size used for the travel time table and other information in the embodiment of the present invention is 50 meters × 100 meters × 20 meters, and the number of grid points is nxv=241,nyv=91,nzv=401. The imaging grid size is 12.5 meters × 25 meters × 5 meters, and the number of imaging grid points is nx=500,ny=60,nz=601.

[0082] S203. According to the seismic Boschhof migration imaging theory, the imaging value img(ix,iy,iz;k;θ) and the illumination times ni(ix,iy,iz;k;θ) of each underground imaging point are calculated, where θ is the illumination inclination angle. When the illumination inclination angle is the true inclination angle or the apparent inclination angle of the illumination surface in the x direction, the true inclination angle expression is:

[0083]

[0084] The expression of apparent inclination is:

[0085]

[0086] When the lighting inclination angle is the true inclination angle or apparent inclination angle of the lighting surface in the y direction, the true inclination angle expression is:

[0087]

[0088] The expression of apparent inclination is:

[0089]

[0090] Where, pxs = sinα SM cosβ SM , pys = sinα SM sinβ SM ,

[0091] pxg=sinα GM cosβ GM , pyg = sinα GM sinβ GM ;

[0092] pxs, pys, pxg, and pyg are all ray parameters. Specifically, pxs and pys are the ray parameters of the x and y directions of the ray from the shot point S to the imaging point M; pxg and pyg are the ray parameters of the x and y directions of the ray from the receiving point R to the imaging point M. α represents the scattering angle of each seismic trace, and cos 2α = cosα SM cosα GM +sinα SM sinα GM cos(β SM -β GM ), α SM and β SM are the polar angle and azimuth angle of the ray from the shot point S to the imaging point M at the imaging point M; α GM and β GM They are the polar angle and azimuth angle of the ray from the receiving point R to the imaging point M at the imaging point M respectively.

[0093] S204, repeating steps S201-S203 until all seismic traces have completed migration imaging.

[0094] S205, Calculate the depth domain angle gather

[0095] I(ix,iy,iz;θ)=∑ k img(ix,iy,iz;k;θ) / ∑ kni(ix,iy,iz;k;θ). In this way, there are nx×ny depth domain angle gathers, which are abbreviated as I ij (z, θ), (i=0,1,…,nx-1;j=0,1,…,ny-1). In the embodiment of the present invention, nx=500,ny=60.

[0096] S102, Sparse Constrained Radon Transform

[0097] For each depth domain angle gather I ij (z,θ), (i=0,1,…,nx-1;j=0,1,…,ny-1), perform sparse constrained Radon transform, which is a special Radon transform, which transforms two-dimensional data into three-dimensional model data, and applies L1 norm sparse constraint to the Radon domain model, using L0 norm as the sparsity criterion of the model. For the convenience of description, this depth domain angle gather is represented as d(z,θ), discretized into d(iz,iθ), (iz=0,1,…,nz-1;iθ=0,1,…,nθ-1), where z=z min +iz·dz;θ=θ min +iθ·dθ,z min is the minimum depth (m), dz is the depth increment (m), θ min is the minimum lighting inclination angle, dθ is the lighting inclination angle increment. In the embodiment of the present invention, nz=601, dz=5 meters, z min =0 meters, nθ=51, dθ=2°, θ min =-50°. Figure 3 A sparse constrained Radon transform flow chart is shown, and the method includes:

[0098] S301, outer loop iteration initialization

[0099] Use the adjoint Radon transform model as the initial value of the Radon domain model m, where m represents a three-dimensional array The curvature of the smile curve p = p min +ip·dp,p min is the minimum curvature, dp is the increment of discrete curvature, where curvature has the dimension of depth and the unit is km; the position of the vertex of the smile curve is the minimum position of the vertex, is the increment of the vertex position, The range and increment of is generally the same as the discrete range and increment of the lighting inclination angle θ, that is, In this embodiment, p min =-2.47, p max =0.51, dp=0.046, np=65. The following are the steps for calculating the adjoint Radon transform:

[0100] (a) Fourier transform of depth: Input a depth domain angle gather d(iz,iθ), (iz=0,1,…,nz-1;iθ=0,1,…,nθ-1), calculate in is the Fourier transform operator of the depth z coordinate, and d1 is the two-dimensional array d1(ikz,iθ) in the depth wavenumber domain (ikz=0,1,…,nkz-1;iθ=0,1,…,nθ-1). z =ikz·dkz, wave number increment dkz=1 / (nzfft·dz), where nzfft is the number of depth discrete points required for Fourier transform, and dz is the depth sampling interval.

[0101] (b) Radon transform in the depth wavenumber domain: For each depth wavenumber ikz, calculate d2 = L T d1, where T represents the matrix transpose and L is the Radon transform kernel function The decision matrix is ​​expressed as:

[0102]

[0103] where k z =ikz·dkz;p ip =p min +ip·dp,(ip=0,1,…,np-1);θ iθ =θ min +iθ·dθ,(iθ=0,1,…,nθ-1); Radon transform kernel function for:

[0104]

[0105] When θ and In the range of ±57.2°, It can be approximated as:

[0106]

[0107] (c) Perform inverse Fourier transform on the depth and calculate

[0108] (d) Take d3 as the initial value of the Radon domain model m of the outer loop.

[0109] S302, outer loop iteration: used to update the model weight matrix, implement L1 norm sparsity constraint, and improve resolution.

[0110] S401, calculate model weight matrix W m The diagonal elements of :

[0111]

[0112] Among them, [W m ] ii The model weight matrix W represents the three-dimensional model of the Radon domain m The i-th diagonal element of the weight matrix in . In order to calculate [W m ] ii , which is expressed as and Where b is the mean value of m in the three-dimensional model of the Ladon domain, Indicates that the sequence numbers in the three-dimensional model m of the Ladon domain are iz, ip and The three-dimensional model value of the Ladon domain.

[0113] S402, PCG inner loop iteration: used to update the Radon domain model

[0114] (a) Inner loop initialization: Assume that the initial value of the inner loop model is zero, that is, m0 = 0, the initial residual of the inner loop is r0 = d3, and the weighted residual is The conjugate gradient is p0=z0;

[0115] (b) Inner loop iteration (the number of iterations is k=0, 1, ..., imax-1, where imax is the maximum number of iterations of the inner loop), in this embodiment, imax=5:

[0116] (b.1) Calculate the model iteration step in I is the identity matrix, γ is the pre-whitening parameter;

[0117] (b.2) Update model m k+1 =m k +α k p k ;

[0118] (b.3) Update residual

[0119] (b.4) Calculate the weighted residual

[0120] (b.5) Calculate the conjugate gradient iteration step size

[0121] (b.6) Update the conjugate gradient p k+1 =z k+1 +β k p k .

[0122] S403, using the residual as the inner loop iteration stopping criterion. When the preset condition is met, the preset condition in this embodiment can be when When ieps is the parameter specified by the user, it is usually 0.001.

[0123] S303, using the L0 norm as the outer loop iteration stop criterion, that is, when the L0 norm of the calculated model is less than oesp times the L0 norm of the initial model, the outer loop iteration is stopped, where oeps is a parameter specified by the user, generally 0.001.

[0124] S304, output the Radon domain model m. ij Depth domain angle gather I ij Radon domain model of (z,θ), where i=0,1,…,nx-1; j=0,1,…,ny-1.

[0125] S103, Diffraction Wave Imaging

[0126] Three-dimensional Radon domain model ij The model traces with curvature p = 0 contain the diffraction wave imaging energy. These traces can be directly superimposed to obtain the diffraction wave imaging.

[0127]

[0128] Where i=0,1,…,nx-1; j=0,1,…,ny-1, ip0=-p min / dp. In this embodiment, ip0=54.

[0129] Taking into account the influence of practical factors such as imaging speed, the present invention also provides another method for calculating diffraction wave imaging, which is to take the three-dimensional Radon domain model m ij The model track near the mean curvature p = 0 is subjected to inverse Radon transformation to obtain the illumination inclination-depth domain data, which are then superimposed to obtain the diffraction wave imaging. Figure 4 A flow chart of diffraction wave imaging is shown, and the method comprises:

[0130] S501, extracting the three-dimensional Ladon domain model m ij Model path near the mean curvature p = 0 in p a <0 and p b >0 is the curvature p value range specified by the user. In this embodiment, p a =-0.2, p b =0.2.

[0131] S502, perform Fourier transform on the depth coordinates of each channel to obtain depth wave number domain data

[0132] S503, Radon inverse transform calculation:

[0133] S504, yes Perform inverse Fourier transform in the depth domain to obtain d ij (z,θ), where z = z min +iz·dz;θ=θ min +iθ·dθ.

[0134] S505, superimpose and output diffraction wave imaging trace:

[0135]

[0136] Figure 2 A flow chart for generating a depth domain angle gather is shown, and the method comprises: inputting a seismic record, reading travel time and other information from a travel time and other table, Kirchhoff integral imaging, continuing with the next imaging, and outputting a depth domain angle gather after all traces have been imaged.

[0137] Figure 3 The sparse constrained Radon transform flow chart is shown, and the method includes: inputting a depth domain angle track set to initialize the outer loop iteration, the outer loop iteration updates the model weight matrix, calculates the diagonal elements of the model weight matrix, the PCG inner loop iteration updates the Radon domain model, uses the residual as the inner loop iteration stopping criterion, uses the L0 norm as the outer loop iteration stopping criterion, and outputs the Radon domain model.

[0138] Figure 4 A diffraction wave imaging flow chart is shown, and the method includes: extracting model traces near zero curvature in a three-dimensional Radon domain model, Fourier transform in depth domain, inverse Radon transform in depth wavenumber domain, inverse Fourier transform in depth domain, and superimposing and outputting diffraction wave imaging traces.

[0139] Compared with the prior art, the embodiment of the present invention adopts the scheme of imaging first and then separation to solve the difficulty of effectively imaging the diffraction wave with weak energy in complex structural areas, and implements three-dimensional sparse constrained Radon transform on the depth domain angle gather, effectively separating the diffraction wave and imaging it. The embodiment of the present invention utilizes the linear characteristics of the diffraction wave in the depth domain illumination dip angle gather, and effectively maintains the waveform characteristics of the diffraction wave through three-dimensional sparse constrained Radon transform, thereby improving the resolution and fidelity of diffraction wave imaging.

[0140] Figure 5It is a depth domain angle gather map generated in an embodiment of the present invention, which has 51 channels, each channel has 601 sampling points, the depth sampling interval is 5 meters, the minimum illumination inclination angle is -50°, the maximum illumination inclination angle is +50°, and the illumination inclination angle increment is 2°. The depth domain angle gather contains two reflection waves and one diffraction wave, in which the reflection wave event axis is presented as a smile curve, and the diffraction wave event axis is presented as a horizontal straight line.

[0141] Figure 6 is the depth domain angle gather diagram after Radon forward and inverse transformation in the embodiment of the present invention, wherein the parameter used for Radon forward transformation is p min =-2.47, p max =0.51, dp = 0.046. The parameters used for the Radon inverse transform are p a =-2.47, p b =0.51.

[0142] Depend on Figure 5 and Figure 6 It can be seen that the depth domain angle gathers after the Radon forward and inverse transforms in the embodiment of the present invention are almost identical to the input angle gathers, the positive and negative lobes of the waveform remain unchanged, and the amplitudes are almost equal, which indicates that the sparse constrained Radon transform and inverse Radon transform proposed in the present invention are fidelity-preserving.

[0143] Figure 7 is the separated diffraction wave depth domain angular gather diagram in the embodiment of the present invention, wherein the parameter used for the Radon inverse transform is p a =-0.2, p b =0.2. It can be seen that the separated diffraction wave waveform is well maintained, the amplitude and phase remain basically unchanged, and the gather contains almost no reflection wave information. This shows that the diffraction wave is separated with high resolution and the waveform of the diffraction wave is maintained.

[0144] Figure 8 3 is a comparison diagram of the reflected wave and diffraction wave imaging (left) and the diffraction wave imaging (right) in an embodiment of the present invention. The left figure includes the imaging of two reflected waves and one diffraction wave, and the right figure only has one diffraction wave imaging, wherein the amplitude of the diffraction wave is much larger than the residual value of the reflected wave, and the diffraction wave waveform is well maintained.

[0145] Fig. 9 The conventional Radon transform zp gather diagram (left) and The left picture is The zp gather on the right is for p = 0 km It can be seen that the energy distribution range of the reflected wave and diffracted wave of the conventional Radon transform is wide, which is not convenient for high-resolution separation of diffracted waves.

[0146] Fig.10The sparse constrained Radon transform zp gather diagram (left) and The left image corresponds to The zp gather of the right image corresponds to the two reflected wave phase axes in the depth domain angle gather. In the gather, the horizontal event axis corresponds to the diffraction wave event axis in the depth domain angular gather.

[0147] contrast Fig. 9 and Fig.10 , indicating that the sparse constrained Radon transform proposed in the present invention separates the diffraction waves with high resolution. This embodiment shows that the present invention effectively improves the resolution and fidelity of diffraction wave imaging through high-resolution Radon transform of depth domain angular gathers. The embodiment of the present invention specifically relates to a method of migration imaging in the field of seismic data processing in geophysical oil and gas exploration, and is applicable to seismic exploration of carbonate oil and gas reservoirs, seismic exploration of igneous oil and gas reservoirs, seismic exploration of small fault block oil and gas reservoirs, and seismic exploration of lithologic oil and gas reservoirs.

[0148] In summary, Figure 6 and Figure 8 The effect of this embodiment is shown. Figure 9-10 The effect of this embodiment is further illustrated by comparison. These figures illustrate the effectiveness and practicality of the method of the present invention in improving the resolution and fidelity of diffraction wave imaging.

[0149] In the embodiment of the present invention, the imaging value and illumination times of each seismic trace at each imaging point underground are calculated according to the seismic recording data, the shot point position and the receiving point position of each seismic trace in the seismic data; the depth domain angle trace set corresponding to each seismic trace is calculated according to the imaging value and illumination times of each seismic trace at each imaging point underground; the sparse constrained Radon transform is performed on multiple depth domain angle trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle trace set into three-dimensional model data; the multiple Radon domain three-dimensional models are performed The diffraction wave imaging data is obtained by superposition, so that the linear characteristics of the diffraction wave can be concentrated by calculating the imaging value and the number of illuminations, which can avoid the problem that the formation dip information cannot be retained in the prior art, thereby affecting the diffraction wave separation effect, and improving the resolution of the diffraction wave imaging; the waveform characteristics of the diffraction wave are effectively maintained through the sparse constrained Radon transform, and since the wavelet deconvolution processing is no longer needed, the problem that the uncertainty of the diffraction wave imaging is increased due to the uncertainty of the wavelet in the prior art, thereby reducing the fidelity of the diffraction wave imaging is avoided, and the fidelity of the diffraction wave imaging is improved.

[0150] The present invention also provides a device in the embodiment, as described in the following embodiment. Since the principle of the device to solve the problem is similar to the method, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be repeated.

[0151] An embodiment of the present invention provides a diffraction wave imaging device, such as Fig.12 As shown, the device may include:

[0152] The imaging value and illumination times calculation module 01 is used to calculate the imaging value and illumination times of each imaging point of each seismic channel in the underground according to the seismic record data, the shot point position and the receiving point position of each seismic channel in the seismic data;

[0153] The depth domain angle gather calculation module 02 is used to calculate the depth domain angle gather corresponding to each seismic trace according to the imaging value and illumination times of each imaging point of each seismic trace in the underground;

[0154] The sparse constrained Radon transform module 03 is used to perform sparse constrained Radon transform on multiple depth domain angle gathers to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle gathers into three-dimensional model data;

[0155] The model superposition module 04 is used to superimpose multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data.

[0156] In one embodiment, it also includes: an illumination inclination angle calculation module, which is used to calculate the illumination inclination angle of each seismic trace in the seismic data according to the travel time from the shot point and the receiving point in each seismic trace in the seismic data to each imaging point underground, the ray azimuth angle, the ray polar angle and the Jacobian determinant element value;

[0157] The imaging value and illumination times calculation module is specifically used to calculate the imaging value and illumination times of each seismic channel at each underground imaging point based on the illumination inclination of each seismic channel in the seismic data, seismic record data, shot point location, receiving point location, travel time from the shot point and receiving point to each underground imaging point, ray azimuth, ray polar angle and Jacobian determinant element value.

[0158] In one embodiment, the illumination inclination angle calculation module is specifically used to calculate the illumination inclination angle of each seismic trace in the seismic data according to the following formula:

[0159] When the lighting inclination angle is the true inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is:

[0160]

[0161] When the lighting inclination angle is the apparent inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is:

[0162]

[0163] When the lighting inclination angle is the true inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is:

[0164]

[0165] When the lighting inclination angle is the apparent inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is:

[0166]

[0167] Where, pxs = sinα SM cosβ SM , pys = sinα SM sinβ SM ,

[0168] pxg=sinα GM cosβ GM , pyg = sinα GM sinβ GM ;

[0169] pxs, pys, pxg, and pyg are all ray parameters. Specifically, pxs and pys are the ray parameters of the x and y directions of the ray from the shot point S to the imaging point M; pxg and pyg are the ray parameters of the x and y directions of the ray from the receiving point R to the imaging point M. α represents the scattering angle of each seismic trace, and cos 2α = cosα SM cosα GM +sinα SM sinα GM cos(β SM -β GM ), α SM and β SM are the polar angle and azimuth angle of the ray from the shot point S to the imaging point M at the imaging point M; α GM and β GM They are the polar angle and azimuth angle of the ray from the receiving point R to the imaging point M at the imaging point M respectively.

[0170] In one embodiment, the depth domain angle gather calculation module is specifically used to calculate the depth domain angle gather corresponding to each seismic trace according to the following formula:

[0171]

[0172] Wherein, I(ix,iy,iz;θ) represents the set of nx×ny depth domain angular trace sets corresponding to each seismic trace; img(ix,iy,iz;k;θ) represents the imaging value of the k-th seismic trace at each underground imaging point; ni(ix,iy,iz;k;θ) represents the number of illuminations of the k-th seismic trace at each underground imaging point, ix, iy and iz are the numbers of the underground imaging points on the x, y and z coordinate axes respectively; nx, ny and nz are the total number of discrete points of the underground imaging points on the x, y and z coordinate axes respectively.

[0173] In one embodiment, the sparse constrained Radon transform module is specifically used to:

[0174] For each depth domain angle gather, Fourier transform, Radon transform and inverse Fourier transform are performed on the depth domain angle gather to obtain the initial value of the Radon domain three-dimensional model;

[0175] Calculate the diagonal elements of the weight matrix of the Radon domain three-dimensional model according to the initial value of the Radon domain three-dimensional model;

[0176] According to the diagonal elements of the weight matrix of the Radon domain 3D model, the L1 norm sparse constraint is performed on the initial value of the Radon domain 3D model to obtain the Radon domain 3D model corresponding to the depth domain angle gather.

[0177] In one embodiment, the sparse constrained Radon transform module is specifically used to:

[0178] For each depth-domain angle gather, Fourier transform is performed on the depth-domain angle gather to obtain a two-dimensional array of the depth-domain angle gather in the depth wavenumber domain;

[0179] The Radon transform and inverse Fourier transform are performed on the two-dimensional array in the depth wavenumber domain to obtain the initial value of the three-dimensional model in the Radon domain.

[0180] In one embodiment, the sparse constrained Radon transform module is specifically used to calculate the diagonal elements of the weight matrix of the Radon domain three-dimensional model according to the initial value of the Radon domain three-dimensional model according to the following formula:

[0181]

[0182] Among them, [W m ] ii The model weight matrix W represents the three-dimensional model of the Radon domain m The i-th diagonal element of the weight matrix in is denoted by b is the mean value of m in the three-dimensional model of the Ladon domain, Can be written as m i , Indicates that the sequence numbers in the three-dimensional model m of the Ladon domain are iz, ip and The three-dimensional model value of the Ladon domain.

[0183] In one embodiment, the model overlay module is specifically used to:

[0184] Extracting model track data with a curvature within a preset value range in each Ladong domain three-dimensional model;

[0185] Performing Fourier transform processing on the depth coordinates of the seismic traces in each Ladon domain three-dimensional model to obtain depth wavenumber domain data of multiple Ladon domain three-dimensional models;

[0186] Performing Radon inverse transform calculation and depth domain Fourier inverse transform calculation on the extracted model trace data with curvature within a preset numerical range in each Radon domain three-dimensional model and the depth wavenumber domain data of multiple Radon domain three-dimensional models to obtain diffraction wave imaging data corresponding to each seismic trace;

[0187] The diffraction wave imaging data corresponding to each seismic trace is stacked and calculated to obtain the diffraction wave imaging data.

[0188] An embodiment of the present invention further provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the above-mentioned diffraction wave imaging method when executing the computer program.

[0189] An embodiment of the present invention further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program for executing the above-mentioned diffraction wave imaging method.

[0190] In the embodiment of the present invention, the imaging value and illumination times of each seismic trace at each imaging point underground are calculated according to the seismic recording data, the shot point position and the receiving point position of each seismic trace in the seismic data; the depth domain angle trace set corresponding to each seismic trace is calculated according to the imaging value and illumination times of each seismic trace at each imaging point underground; the sparse constrained Radon transform is performed on multiple depth domain angle trace sets to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle trace set into three-dimensional model data; the multiple Radon domain three-dimensional models are performed The diffraction wave imaging data is obtained by superposition, so that the linear characteristics of the diffraction wave can be concentrated by calculating the imaging value and the number of illuminations, which can avoid the problem that the formation dip information cannot be retained in the prior art, thereby affecting the diffraction wave separation effect, and improving the resolution of the diffraction wave imaging; the waveform characteristics of the diffraction wave are effectively maintained through the sparse constrained Radon transform, and since the wavelet deconvolution processing is no longer needed, the problem that the uncertainty of the diffraction wave imaging is increased due to the uncertainty of the wavelet in the prior art, thereby reducing the fidelity of the diffraction wave imaging is avoided, and the fidelity of the diffraction wave imaging is improved.

[0191] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0192] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0193] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0194] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0195] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A diffraction wave imaging method, characterized in that: include: According to the seismic record data, shot point position and receiving point position of each seismic channel in the seismic data, the imaging value and illumination times of each imaging point of each seismic channel in the underground are calculated; According to the imaging value and illumination times of each seismic trace at each imaging point underground, the depth domain angle gather corresponding to each seismic trace is calculated; Perform sparse constrained Radon transform on multiple depth domain angle gathers to obtain multiple Radon domain 3D models; The sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle gather into three-dimensional model data; Multiple Ladon domain three-dimensional models are superimposed to obtain diffraction wave imaging data; Sparse constrained Radon transform is performed on multiple depth domain angle gathers to obtain multiple Radon domain 3D models, including: For each depth domain angle gather, Fourier transform, Radon transform and inverse Fourier transform are performed on the depth domain angle gather to obtain the initial value of the Radon domain three-dimensional model; Calculate the diagonal elements of the weight matrix of the Radon domain three-dimensional model according to the initial value of the Radon domain three-dimensional model; According to the diagonal elements of the weight matrix of the Radon domain 3D model, the initial value of the Radon domain 3D model is subjected to L1 norm sparse constraint, and the Radon domain 3D model corresponding to the depth domain angle gather is obtained; The diagonal elements of the weight matrix of the Radon domain 3D model are calculated according to the initial values ​​of the Radon domain 3D model as follows: Among them, [W m ] ii The model weight matrix W represents the three-dimensional model of the Radon domain m The i-th diagonal element in b is the mean value of m in the three-dimensional model of the Ladon domain, Indicates that the sequence numbers in the three-dimensional model m of the Ladon domain are iz, ip and The three-dimensional model value of the Radon domain can be recorded as m i ; z is the imaging depth, z = z min +iz·dz,(iz=0,1,…,nz-1),z min is the minimum imaging depth, dz is the sampling interval of imaging depth z, nz is the number of sampling points of imaging depth z, and the unit of imaging depth z is km; p is the event axis curvature, p = p min +ip·dp,(ip=0,1,…,np-1),p min is the minimum curvature, dp is the sampling interval of curvature p, np is the number of sampling points of curvature p, and the unit of curvature p is also km; is the formation inclination, is the minimum formation dip, is the formation dip The sampling interval is is the formation dip The number of sampling points, the formation dip The unit is degree.

2. The method according to claim 1, characterized in that Also includes: Calculate the illumination dip angle of each seismic trace in the seismic data according to the travel time from the shot point and the receiving point in each seismic trace to each underground imaging point, the ray azimuth, the ray polar angle and the element value of the Jacobian determinant; According to the seismic record data, shot point location and receiving point location of each seismic channel in the seismic data, the imaging value and illumination times of each imaging point of each seismic channel in the underground are calculated, including: Based on the illumination inclination of each seismic channel in the seismic data, seismic record data, shot point location, receiving point location, travel time from the shot point and receiving point to each underground imaging point, ray azimuth, ray polar angle and Jacobian determinant element values, the imaging value and illumination times of each seismic channel at each underground imaging point are calculated.

3. The method according to claim 2, characterized in that The illumination inclination angle of each seismic channel in the seismic data is calculated according to the following formula: When the lighting inclination angle is the true inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is: When the lighting inclination angle is the apparent inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is: When the lighting inclination angle is the true inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is: When the lighting inclination angle is the apparent inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is: Where, pxs = sinα SM cosβ SM , pys = sinα SM sinβ SM , pxg=sinα GM cosβ GM , pyg = sinα GM sinβ GM ; pxs, pys, pxg, pyg are all ray parameters, pxs and pys are the ray parameters of the x-direction and y-direction of the ray from the shot point s to the imaging point M; pxg and pyg are the ray parameters of the x-direction and y-direction of the ray from the receiving point R to the imaging point M; α represents the scattering angle of each seismic trace, and cos 2α=cosα SM cosα GM +sinα SM sinα GM cos(β SM -β GM ), α SM and β SM are the polar angle and azimuth angle of the ray from the shot point S to the imaging point M at the imaging point M; α GM and β GM They are the polar angle and azimuth angle of the ray from the receiving point R to the imaging point M at the imaging point M respectively.

4. The method according to claim 1, characterized in that The depth domain angle gather corresponding to each seismic trace is calculated according to the following formula: (ix=0,1,…,nx-1; iy=0,1,…,ny-1; iz=0,1,…,nz-1) Where I(ix,iy,iz;θ) represents the set of nx×ny depth domain angle gathers corresponding to each seismic trace; img(ix,iy,iz;k;θ) represents the imaging value of the k-th seismic trace at each underground imaging point; ni(ix,iy,iz;k;θ) represents the illumination times of the k-th seismic trace at each underground imaging point, ix, iy and iz are the numbers of the underground imaging points on the x, y and z axes respectively; nx, ny and nz are the total number of discrete points of the underground imaging points on the x, y and z axes respectively; θ is the illumination inclination; when θ=θ min +iθ·dθ,θ min is the minimum lighting inclination angle, dθ is the sampling interval of the lighting inclination angle θ, nθ is the number of sampling points of the lighting inclination angle θ, and the unit of the lighting inclination angle θ is radian; 5. The method according to claim 1, characterized in that For each depth domain angle gather, Fourier transform, Radon transform and inverse Fourier transform are performed on the depth domain angle gather to obtain the initial value of the Radon domain three-dimensional model, including: For each depth-domain angle gather, Fourier transform is performed on the depth-domain angle gather to obtain a two-dimensional array of the depth-domain angle gather in the depth wavenumber domain; The Radon transform and inverse Fourier transform are performed on the two-dimensional array in the depth wavenumber domain to obtain the initial value of the three-dimensional model in the Radon domain.

6. The method according to claim 1, characterized in that Multiple Ladon domain 3D models are superimposed to obtain diffraction wave imaging data, including: Extracting model track data with a curvature within a preset value range in each Ladong domain three-dimensional model; Performing Fourier transform processing on the depth coordinates of the seismic traces in each Ladon domain three-dimensional model to obtain depth wavenumber domain data of multiple Ladon domain three-dimensional models; Performing Radon inverse transform calculation and depth domain Fourier inverse transform calculation on the extracted model trace data with curvature within a preset numerical range in each Radon domain three-dimensional model and the depth wavenumber domain data of multiple Radon domain three-dimensional models to obtain diffraction wave imaging data corresponding to each seismic trace; The diffraction wave imaging data corresponding to each seismic trace is stacked and calculated to obtain the diffraction wave imaging data.

7. A diffraction wave imaging device, characterized in that: include: An imaging value and illumination times calculation module is used to calculate the imaging value and illumination times of each imaging point of each seismic channel in the underground according to the seismic record data, the shot point position and the receiving point position of each seismic channel in the seismic data; A depth domain angle gather calculation module is used to calculate the depth domain angle gather corresponding to each seismic trace according to the imaging value and illumination times of each imaging point of each seismic trace in the underground; A sparse constrained Radon transform module is used to perform sparse constrained Radon transform on multiple depth domain angle gathers to obtain multiple Radon domain three-dimensional models; the sparse constrained Radon transform is used to transform the two-dimensional data represented by the depth domain angle gathers into three-dimensional model data; A model superposition module is used to superimpose multiple Ladon domain three-dimensional models to obtain diffraction wave imaging data; Sparse constrained Radon transform module, specifically used for: For each depth domain angle gather, Fourier transform, Radon transform and inverse Fourier transform are performed on the depth domain angle gather to obtain the initial value of the Radon domain three-dimensional model; Calculate the diagonal elements of the weight matrix of the Radon domain three-dimensional model according to the initial value of the Radon domain three-dimensional model; According to the diagonal elements of the weight matrix of the Radon domain 3D model, the initial value of the Radon domain 3D model is subjected to L1 norm sparse constraint, and the Radon domain 3D model corresponding to the depth domain angle gather is obtained; Sparse constrained Radon transform module, specifically used for: The diagonal elements of the weight matrix of the Radon domain 3D model are calculated according to the initial values ​​of the Radon domain 3D model as follows: Among them, [W m ] ii The model weight matrix W represents the three-dimensional model of the Radon domain m The i-th diagonal element in b is the mean value of m in the three-dimensional model of the Ladon domain, Indicates that the sequence numbers in the three-dimensional model m of the Ladon domain are iz, ip and The three-dimensional model value of the Radon domain can be recorded as m i ; z is the imaging depth, z = z min +iz·dz,(iz=0,1,…,nz-1),z min is the minimum imaging depth, dz is the sampling interval of imaging depth z, nz is the number of sampling points of imaging depth z, and the unit of imaging depth z is km; p is the event axis curvature, p = p min +ip·dp,(ip=0,1,…,np-1),p min is the minimum curvature, dp is the sampling interval of curvature p, n is the number of sampling points of curvature p, and the unit of curvature p is also km; is the formation inclination, is the minimum formation dip, is the formation dip The sampling interval is is the formation dip The number of sampling points, the formation dip The unit is degree.

8. The device according to claim 7, characterized in that Also includes: The illumination inclination angle calculation module is used to calculate the illumination inclination angle of each seismic channel in the seismic data according to the travel time from the shot point and the receiving point in each seismic channel in the seismic data to each underground imaging point, the ray azimuth angle, the ray polar angle and the Jacobian determinant element value; The imaging value and illumination times calculation module is specifically used to calculate the imaging value and illumination times of each seismic channel at each underground imaging point based on the illumination inclination of each seismic channel in the seismic data, seismic record data, shot point location, receiving point location, travel time from the shot point and receiving point to each underground imaging point, ray azimuth, ray polar angle and Jacobian determinant element value.

9. The device according to claim 8, characterized in that The illumination inclination angle calculation module is specifically used to calculate the illumination inclination angle of each seismic channel in the seismic data according to the following formula: When the lighting inclination angle is the true inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is: When the lighting inclination angle is the apparent inclination angle of the lighting surface in the x direction, the expression for calculating the lighting inclination angle θ is: When the lighting inclination angle is the true inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is: When the lighting inclination angle is the apparent inclination angle of the lighting surface in the y direction, the expression for calculating the lighting inclination angle θ is: where pxs = sinα SM cosβ SM , pys = sinα SM sinβ SM , pxg=sinα GM cosβ GM ,pyg=sinα GM sinβ GM ; pxs, pys, pxg, and pyg are all ray parameters. pxs and pys are the ray parameters of the x-direction and y-direction of the ray from the shot point S to the imaging point M; pxg and pyg are the ray parameters of the x-direction and y-direction of the ray from the receiving point R to the imaging point M; α represents the scattering angle of each seismic trace, and cos 2α=cosα SM cosα GM +sinα SM sinα GM cos(β SM -β GM ), α SM and β SM are the polar angle and azimuth angle of the ray from the shot point S to the imaging point M at the imaging point M; α GM and β GM They are the polar angle and azimuth angle of the ray from the receiving point R to the imaging point M at the imaging point M respectively.

10. The device according to claim 7, characterized in that The depth domain angle gather calculation module is specifically used to calculate the depth domain angle gather corresponding to each seismic trace according to the following formula: (ix=0,1,…,nx-1; iy=0,1,…,ny-1; iz=0,1,…,nz-1) Where I(ix,iy,iz;θ) represents the set of nx×ny depth domain angle gathers corresponding to each seismic trace; img(ix,iy,iz;k;θ) represents the imaging value of the k-th seismic trace at each underground imaging point; ni(ix,iy,iz;k;θ) represents the number of illuminations of the k-th seismic trace at each underground imaging point, ix, iy and iz are the numbers of the underground imaging points on the x, y and z axes respectively; nx, ny and nz are the total number of discrete points of the underground imaging points on the x, y and z axes respectively; θ is the illumination inclination, θ=θ min +iθ·dθ,θ min is the minimum lighting inclination angle, dθ is the sampling interval of the lighting inclination angle θ, nθ is the number of sampling points of the lighting inclination angle θ, and the unit of the lighting inclination angle θ is radian; 11. The device according to claim 7, characterized in that Sparse constrained Radon transform module, specifically used for: For each depth-domain angle gather, Fourier transform is performed on the depth-domain angle gather to obtain a two-dimensional array of the depth-domain angle gather in the depth wavenumber domain; The Radon transform and inverse Fourier transform are performed on the two-dimensional array in the depth wavenumber domain to obtain the initial value of the three-dimensional model in the Radon domain.

12. The device according to claim 7, characterized in that Model overlay module, specifically used for: Extracting model track data with a curvature within a preset value range in each Ladong domain three-dimensional model; Performing Fourier transform processing on the depth coordinates of the seismic traces in each Ladon domain three-dimensional model to obtain depth wavenumber domain data of multiple Ladon domain three-dimensional models; Performing Radon inverse transform calculation and depth domain Fourier inverse transform calculation on the extracted model trace data with curvature within a preset numerical range in each Radon domain three-dimensional model and the depth wavenumber domain data of multiple Radon domain three-dimensional models to obtain diffraction wave imaging data corresponding to each seismic trace; The diffraction wave imaging data corresponding to each seismic trace is stacked and calculated to obtain the diffraction wave imaging data.

13. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the method according to any one of claims 1 to 6 is implemented.

14. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program for executing the method according to any one of claims 1 to 6.

Images