Scattered wave imaging method and device

Abstract

The present invention discloses a scattered wave imaging method and device, wherein the method comprises: calculating the kernel function of the angle gather generated by seismic data; calculating the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather; performing L1 norm and L2 norm sparse inversion on the angle gather according to the generalized Fourier slice theorem interpolation operator of the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather; obtaining the scattered wave imaging result according to the three-dimensional model representing the two-dimensional array of the angle gather. The present invention can realize high-resolution imaging of scattered waves of the angle gather; and at the same time improve the calculation efficiency of scattered wave imaging.

Description

Technical Field

[0001] The present invention relates to the technical field of geophysical oil and gas exploration, and in particular to a scattered 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] Seismic exploration is a method that uses the differences in elasticity and density of underground media to observe and analyze the response of the earth to artificially excited seismic waves to infer the underground structural morphology and rock properties. It has become the most important means of oil and natural gas resource exploration and is also widely used in coalfield and engineering geological surveys, regional geological research, and crustal research. At present, oil and gas resource exploration mainly uses reflection wave seismic exploration, which consists of three parts: seismic data acquisition, processing, and interpretation. Among them, seismic data migration imaging is the key core link of seismic data processing.

[0004] 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 scattered waves. The use of scattered waves is expected to improve the imaging resolution of seismic data.

[0005] Within the wavelength scale of seismic waves, the velocity and density inhomogeneities of underground media will produce scattering phenomena, such as the reflection of rough interfaces, which is usually called diffuse reflection. The energy of reflected waves is concentrated in a certain direction, while the energy of scattered waves is scattered in all directions. Therefore, the energy of scattered waves is weaker than that of reflected waves, which makes scattered wave imaging much more difficult than reflection wave imaging. Prestack seismic data migration imaging gathers are usually depth-offset distance domain gathers, which are easy to convert into depth-scattering angle domain gathers. When the migration speed is appropriate, the reflected waves appear as horizontal straight line events on these two gathers, while the scattered waves appear as convex curves. After superposition, only the reflected waves can be strengthened and the scattered waves can be suppressed. If the depth-illumination dip angle domain gathers are generated during the migration process, the reflected waves appear as concave smile curve events, while the scattered waves appear as horizontal straight line events. Therefore, it is easier to separate the scattered waves on the depth-illumination dip angle domain gathers (hereinafter referred to as angle gathers), and then scattered wave imaging can be achieved through superposition.

[0006] At present, Radon transform is usually used to separate scattered waves. However, due to the influence of factors such as aperture, the resolution of conventional Radon transform is low, which seriously affects the effect of scattered wave imaging. The Radon transform with sparse constraints of preconditioned conjugate gradient (PCG) has high resolution, but the calculation amount increases by hundreds of times by iteratively solving the Radon transform equation constrained by the model, and the calculation efficiency is low.

[0007] In summary, at present, when performing scattered wave imaging, the resolution is usually low, which seriously affects the effect of scattered wave imaging. At the same time, due to the large amount of calculation, the calculation efficiency of scattered wave imaging is seriously affected. Summary of the invention

[0008] An embodiment of the present invention provides a scattered wave imaging method for improving the resolution and calculation efficiency of scattered wave imaging. The method includes:

[0009] Compute kernel functions for angle gathers generated from seismic data;

[0010] According to the kernel function of the angle gather, a generalized Fourier slice theorem interpolation operator of the angle gather is calculated; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather;

[0011] According to the generalized Fourier slice theorem interpolation operator of the angle gather, the L1 norm and L2 norm sparse inversion of the angle gather is performed to obtain a three-dimensional model representing the two-dimensional array of the angle gather;

[0012] The scattered wave imaging results are obtained based on the three-dimensional model that characterizes the two-dimensional array of angle gathers.

[0013] The embodiment of the present invention further provides a scattered wave imaging device for improving the resolution and calculation efficiency of scattered wave imaging, the device comprising:

[0014] A kernel function calculation module is used to calculate the kernel function of the angle gather generated by seismic data;

[0015] An interpolation operator calculation module is used to calculate a generalized Fourier slice theorem interpolation operator of the angle gather according to a kernel function of the angle gather; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather;

[0016] A sparse inversion module is used to perform L1 norm and L2 norm sparse inversion on the angular gathers according to the generalized Fourier slice theorem interpolation operator of the angular gathers to obtain a three-dimensional model representing the two-dimensional array of the angular gathers;

[0017] The imaging module is used to obtain scattered wave imaging results based on a three-dimensional model representing a two-dimensional array of angular gathers.

[0018] 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 scattered wave imaging method when executing the computer program.

[0019] 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 scattered wave imaging method.

[0020] In an embodiment of the present invention, a kernel function of an angle gather generated by seismic data is calculated; based on the kernel function of the angle gather, a generalized Fourier slice theorem interpolation operator of the angle gather is calculated; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather; based on the generalized Fourier slice theorem interpolation operator of the angle gather, L1 norm and L2 norm sparse inversion are performed on the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather; based on the three-dimensional model representing the two-dimensional array of the angle gather , and the scattered wave imaging result is obtained. Thus, the generalized Fourier slice theorem and the L1 norm and L2 norm sparse inversion can be combined through the generalized Fourier slice theorem interpolation operator. There is no need to use the conventional Radon transform under the existing technology to avoid being affected by factors such as aperture, thereby avoiding the problem of resolution reduction and achieving high-resolution imaging of scattered waves of angular gathers. At the same time, there is no need to iteratively solve the Radon transform equation constrained by the model under the existing technology, thus avoiding the problem of large amount of calculation and improving the computational efficiency of scattered wave imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] 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:

[0022] Figure 1 A schematic flow chart of an example of a scattered wave imaging method in an embodiment of the present invention;

[0023] Figure 2 It is a kernel function surface diagram in an embodiment of the present invention;

[0024] Figure 3 is a position diagram of the interpolation operator of the generalized Fourier slice theorem in an embodiment of the present invention;

[0025] Figure 4 This is a trend diagram of the inner loop residual of the sparse inversion in an embodiment of the present invention;

[0026] Figure 5 is a trend diagram of the L0 norm of the outer loop of the sparse inversion in an embodiment of the present invention;

[0027] Figure 6 This is a comparison diagram of computing efficiency in an embodiment of the present invention;

[0028] Figure 7 is an angle trace set diagram input in an embodiment of the present invention;

[0029] Figure 8 is a cross-sectional view of pz in an embodiment of the present invention and Comparison of cross-sections;

[0030] Fig. 9 A comparison diagram of scattered wave imaging effects in an embodiment of the present invention;

[0031] Fig.10 Schematic diagram of a scattered wave imaging method according to an embodiment of the present invention;

[0032] Fig.11 Schematic diagram of the structure of a scattered wave imaging device in an embodiment of the present invention. DETAILED DESCRIPTION

[0033] 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.

[0034] Scattering wave imaging technology can be used to identify the locations of lithology mutations such as ancient river channels, lithology pinch-out points, faults, carbonate caves, and igneous rocks, which is helpful for the exploration of carbonate rocks, igneous rocks, small fault blocks, and lithology oil and gas reservoirs. However, the inventors have found that: at present, when performing scattered wave imaging, the resolution is usually low, which seriously affects the effect of scattered wave imaging. At the same time, due to the large amount of calculation, the calculation efficiency of scattered wave imaging is seriously affected. The present invention can be used to solve the above problems and improve the resolution and calculation efficiency of scattered wave imaging. A scattered wave imaging method provided in an embodiment of the present invention is as follows: Fig.10 As shown, the following steps may be included:

[0035] Step 1001: Calculate the kernel function of the angle gather generated by seismic data;

[0036] Step 1002: Calculate the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather;

[0037] Step 1003: performing L1 norm and L2 norm sparse inversion on the angle gathers according to the generalized Fourier slice theorem interpolation operator of the angle gathers to obtain a three-dimensional model representing the two-dimensional array of the angle gathers;

[0038] Step 1004: Obtain scattered wave imaging results based on the three-dimensional model that characterizes the two-dimensional array of angle gathers.

[0039] In specific implementation, the kernel function of the angle gather generated by seismic data is first calculated.

[0040] In the embodiment, the kernel function of the angle gather generated by the seismic data is calculated according to the following formula:

[0041]

[0042] in, represents the kernel function of the angular gather; θ is the illumination inclination angle of the angular gather, in radians, θ = θ min +iθ·dθ,θ min is the minimum illumination inclination angle of the angle gather, in radians, dθ is the sampling interval of the illumination inclination angle, in radians; iθ is the sampling sequence number of the illumination inclination angle, iθ=0,1,…,nθ-1, nθ is the total number of samples of the illumination inclination angle; is the formation dip angle of the angle gather, in radians, is the minimum formation dip angle of the angle gather, in radians, is the sampling interval of formation dip, in radians, is the sampling number of the formation dip, is the total number of samples of formation dip.

[0043] In the above embodiment, by calculating the kernel function of the angle gather, it is helpful to calculate the generalized Fourier slice theorem interpolation operator of the angle gather in the subsequent step to realize scattered wave imaging.

[0044] In a specific implementation, after calculating the kernel function of the angle gather generated by seismic data, the generalized Fourier slice theorem interpolation operator of the angle gather is calculated according to the kernel function of the angle gather; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather.

[0045] In an embodiment, calculating the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather can include: discretizing the wavenumber corresponding to the curvature in the angle gather according to the kernel function of the angle gather; calculating the sequence number of the wavenumber corresponding to the depth in the angle gather according to the wavenumber corresponding to the curvature in the discretized angle gather; calculating the linear interpolation operator coefficient according to the sequence number of the wavenumber corresponding to the depth in the angle gather; and calculating the generalized Fourier slice theorem interpolation operator of the angle gather according to the linear interpolation operator coefficient.

[0046] In the embodiment, the generalized Fourier slice theorem interpolation operator of the angle gather is calculated according to the linear interpolation operator coefficients according to the following formula:

[0047]

[0048] Where S represents the generalized Fourier slice theorem interpolation operator of the angular gather, dimensionless; δ is the coefficient of the linear interpolation operator, δ = (k p -fkp) / dkp-ikp;k p It is the wave number corresponding to the event axis curvature in the angular channel gather, with the unit of 1 / km, referred to as the curvature wave number; is the kernel function of the angle gather, dimensionless; k z k is the wave number corresponding to the depth of the angular channel concentration, the unit is 1 / km, referred to as the depth wave number; z =fkz+ikz·dkz; fkz is the depth wave number k z The minimum value of the depth is expressed in 1 / km; dkz is the depth wave number k z The sampling interval is in 1 / km; ikz is the depth wave number k z The sampling sequence number is ikz = 0, 1, ..., nkz-1; nkz is the depth wave number k z The total number of samples; ikp is the curvature wave number k p The sampling number, represents the integer rounded downward; fkp is the curvature wave number k p The minimum value of the unit is 1 / km, dkp is the curvature wave number k p The sampling interval is 1 / km, ikp = 0, 1, ..., nkp-1, nkp is the curvature wave number k p The total number of samples.

[0049] In the above embodiment, an algorithm for calculating the generalized Fourier slice theorem interpolation operator for the angle gather is provided, which is helpful for calculating the scattered wave imaging result in the subsequent steps.

[0050] In specific implementation, after calculating the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather, the L1 norm and L2 norm sparse inversion are performed on the angle gather according to the generalized Fourier slice theorem interpolation operator of the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather.

[0051] In an embodiment, according to the generalized Fourier slice theorem interpolation operator of the angle gather, L1 norm and L2 norm sparse inversion is performed on the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather, including: according to the generalized Fourier slice theorem interpolation operator of the angle gather, L1 norm and L2 norm optimization problem expressions are established; the L1 norm and L2 norm optimization problem expressions are converted into L1 norm and L2 norm linear equation groups; the L1 norm and L2 norm linear equation groups are solved by the conjugate gradient algorithm to obtain a three-dimensional model representing the two-dimensional array of the angle gather.

[0052] In the embodiment, the L1 norm and L2 norm optimization problem expressions are established according to the following formulas:

[0053]

[0054] Where d is a vector of two-dimensional arrays in the z-θ domain representing the angle gather; m is a vector of two-dimensional arrays d representing the angle gather. A vector of arrays of domain 3D models; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; W m is the weight matrix of the three-dimensional model; represents the L2 norm; represents the L1 norm; k is the preset number of iterations in the process of solving the L1 norm and L2 norm optimization problem expressions, m (k) represents the three-dimensional model m obtained in the kth iteration, m (k+1) represents the three-dimensional model m obtained in the k+1th iteration; represents the three-dimensional model m that minimizes the value of {·}; is the sampling number of the formation dip, is the total number of samples of formation dip.

[0055] In the embodiment, the L1 norm and L2 norm optimization problem expressions can be converted into L1 norm and L2 norm linear equations according to the following formula:

[0056]

[0057] Where A is the cross-correlation matrix of the generalized Fourier slice theorem interpolation operator S, b is the vector obtained by applying the interpolation operator S to d, X is the vector composed of the three-dimensional model m, For the three-dimensional model Quantity; is the interpolation operator S of the generalized Fourier slice theorem The matrix consists of components, that is, where i = 0, 1, ..., nθ, j = 0, 1, ..., nkp, and is the generalized Fourier slice theorem interpolation operator The superscript T indicates the matrix transpose. For the matrix The transpose of W m is the weight matrix of the three-dimensional model. The weight matrix is ​​a diagonal matrix. Except for the diagonal elements, all other elements are zero. The i-th diagonal element of the weight matrix is is the mean value of the three-dimensional model m; m i is the value of the i-th element of the three-dimensional model m; Among them, nz, np and are the depth z, event curvature p and formation dip of the three-dimensional model m respectively. Total number of upsamples;

[0058] is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; F -1 It is an operator for performing two-dimensional inverse Fourier transform on the depth wavenumber kz coordinates and curvature wavenumber kp coordinates in the three-dimensional model in the wavenumber domain.

[0059] In the embodiment, the conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers, including:

[0060] According to the residual and L0 norm between the input angle gather and the pre-calculated angle gather, the conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

[0061] In the above embodiment, by establishing L1 norm and L2 norm optimization problem expressions according to the generalized Fourier slice theorem interpolation operator of the angle gather, the optimal solution of the three-dimensional model can be obtained;

[0062] The L1 norm and L2 norm optimization problem expressions are transformed into L1 norm and L2 norm linear equations, and the optimal solution of the three-dimensional model can be calculated by numerical calculation methods;

[0063] The conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model that represents the two-dimensional array of angle gathers. The optimal solution of the three-dimensional model can be obtained quickly and stably.

[0064] In specific implementation, after performing L1 norm and L2 norm sparse inversion on the angle gathers according to the generalized Fourier slice theorem interpolation operator of the angle gathers to obtain a three-dimensional model representing the two-dimensional array of the angle gathers, the scattered wave imaging result is obtained according to the three-dimensional model representing the two-dimensional array of the angle gathers.

[0065] In the embodiment, the scattered wave imaging result can be obtained according to the following formula based on the three-dimensional model representing the two-dimensional array of angle gathers:

[0066]

[0067] Wherein, I(iz) is the scattered wave imaging result; iz is the depth domain sampling number, iz = 0, 1, ..., nz-1; nz is the total number of depth domain sampling; dz is the depth domain sampling interval, the unit is km; iθ is the sampling number of the illumination inclination, iθ = 0, 1, ..., nθ-1, nθ is the total number of illumination inclination sampling; is the sampling number of the formation dip, is the total number of samples of formation dip; ip is the sampling number of event curvature, ip = 0, 1, ..., np-1; np is the total number of samples of event curvature; p min is the minimum curvature value specified by the user, in km; ip1 is the serial number of the lower curvature value used in scattered wave imaging, ip1 = (p a -p min ) / dp; ip2 is the serial number of the upper curvature value used in scattered wave imaging, ip2 = (p b -p min ) / dp;p a and p bis the preset curvature p value range, p a <0 and p b >0, all units are km; The kernel function that describes the angle gather, I(iz) represents the obtained scattered wave imaging result, iz=0,1,…,nz-1.

[0068] In the above embodiment, a method for obtaining scattered wave imaging results is provided, which avoids the problem of resolution reduction and can achieve high-resolution imaging of scattered waves of angle gathers; at the same time, it avoids the problem of large amount of calculation and improves the calculation efficiency of scattered wave imaging.

[0069] A specific embodiment is given below to illustrate the specific application of the method of the present invention. In this embodiment, Figure 1 As shown, the following steps may be included:

[0070] S101. Kernel function calculation: that is, calculating the kernel function of the angle gather generated by seismic data

[0071] S102, calculation of the generalized Fourier slice theorem interpolation operator: that is, calculating the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather;

[0072] S103, L1-L2 norm sparse inversion: that is, according to the generalized Fourier slice theorem interpolation operator of the angle gather, L1 norm and L2 norm sparse inversion are performed on the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather;

[0073] S104, scattered wave imaging: that is, obtaining scattered wave imaging results according to the three-dimensional model representing the two-dimensional array of angle gathers;

[0074] S105, all angle gathers are calculated through S103 and S104 to complete the scattered wave imaging of the entire work area: that is, the above steps are performed on each angle gather generated by the seismic data.

[0075] S101-S105 in this embodiment are described in detail below:

[0076] S101, kernel function calculation

[0077] Given an angle gather d(iz,iθ), (iz=0,1,…,nz-1;iθ=0,1,…,nθ-1), where nz is the number of sampling points in the depth domain of the angle gather, and nθ is the number of illumination inclination angles. The depth coordinate of the angle gather is z=iz·dz, where dz is the depth domain sampling interval (in km). The illumination inclination angle coordinate of the angle gather is θ=θ min +iθ·dθ, where θ min is the minimum illumination inclination (in radians), and dθ is the illumination inclination sampling interval (in radians). in The kernel function calculation formula is:

[0078]

[0079] For example, the number of sampling points nz of each seismic channel is 601, and the depth sampling interval is 5 meters. The number of illumination inclination channels nθ is 51, and the minimum illumination inclination angle θ min = -50°, the illumination angle sampling interval dθ = 2°. Figure 2 As shown, Figure 2 The kernel function surface diagram in the embodiment of the present invention is shown, indicating that the kernel function surface is shaped like a butterfly wing.

[0080] S102, Generalized Fourier Slice Theorem Interpolation Operator Calculation

[0081] The generalized Fourier slice theorem can be expressed as

[0082]

[0083] in is the kernel function, k z is the wave number corresponding to the depth (unit: 1 / km), k p is the wave number corresponding to the event curvature p (unit: 1 / km). z Discretized into k z =fkz+ikz·dkz, where fkz is k z minimum value (unit: 1 / km), dkz is the sampling interval of the depth wave number (unit: 1 / km), ikz=0,1,…,nkz-1, nkz is the total number of sampling of the depth wave number. p Discretized into k p =fkp+ikp·dkp, where fkp is k p minimum value (unit: 1 / km), dkp is the sampling interval of the event curvature wave number (unit: 1 / km), ikp=0,1,…,nkp-1,nkp is the total number of samplings of the event curvature wave number. In the embodiment of the present invention, fkz=0 (unit: 1 / km), nkz=309, dkz=2.04 (unit: 1 / km), fkp=0 (unit: 1 / km), nkp=130, dkp=1.036 (unit: 1 / km).

[0084] So using a four-dimensional array represents the generalized Fourier slice theorem interpolation operator, and the calculation steps are as follows:

[0085] (a) Calculation

[0086] (b) Calculation in Indicates rounding down to an integer.

[0087] (c) Calculate the linear interpolation operator coefficients

[0088] (d) Calculate the element value of the four-dimensional array S

[0089]

[0090] Where ikz=0,1,…,nkz-1;iθ=0,1,…,nθ-1; like Figure 3 As shown, Figure 3 is the interpolation operator position diagram of the generalized Fourier slice theorem in the embodiment of the present invention, and the interpolation operator position ikp is the depth wave number ikz, the illumination dip angle θ and the formation dip angle This figure corresponds to ikz = 154. When the maximum ikp = 65, the interpolation operator positions greater than 65 will be clipped.

[0091] S103, L1-L2 norm sparse inversion

[0092] Input an angle gather and use the generalized Fourier slice theorem to establish the L1-L2 norm optimization problem expression:

[0093] Where d is a vector of two-dimensional arrays in the z-θ domain representing the angle gather; m is a vector of two-dimensional arrays d representing the angle gather. A vector of arrays of domain 3D models; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; W m is the weight matrix of the three-dimensional model; represents the L2 norm; represents the L1 norm; k is the number of iterations in the optimization problem solving process, for example, m (k) represents the three-dimensional model m obtained in the kth iteration, m (k+1) represents the three-dimensional model m obtained in the k+1th iteration; represents the three-dimensional model m that minimizes the value of {·}; is the sampling number of the formation dip, is the total number of samples of formation dip. This optimization problem is transformed into the following linear equations:

[0094]

[0095] Where A is the cross-correlation matrix of the generalized Fourier slice theorem interpolation operator S, b is the vector obtained by applying the interpolation operator S to d, X is a vector consisting of the three-dimensional model m to be determined, For the three-dimensional model Quantity; is the interpolation operator S of the generalized Fourier slice theorem The matrix consists of components, that is, where i = 0, 1, ..., nθ, j = 0, 1, ..., nkp, and is the generalized Fourier slice theorem interpolation operator The superscript T indicates the matrix transpose. For the matrix The transpose of W m is the weight matrix of the three-dimensional model. The weight matrix is ​​a diagonal matrix. Except for the diagonal elements, all other elements are zero. The i-th diagonal element of the weight matrix is is the mean value of the three-dimensional model m; m i is the value of the i-th element of the three-dimensional model m; Among them, nz, np and are the depth z, event curvature p and formation dip of the three-dimensional model m respectively. Total number of upsamples; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; F -1 is an operator for performing two-dimensional Fourier inverse transform on the depth wave number kz coordinate and the curvature wave number kp coordinate in the three-dimensional model in the wave number domain. The above linear equations are solved using the preconditioned conjugate gradient (PCG) algorithm, where the residual between the input angle gather and the pre-calculated angle gather (abbreviated as residual, the same below) is used as the iteration stop criterion in the inner loop, where the residual between the input angle gather and the pre-calculated angle gather can be expressed as The outer loop uses the L0 norm as the iteration stop criterion. This means that the inner loop guarantees the convergence of the L2 norm, and the outer loop guarantees the convergence of the L1 norm, thus achieving sparse inversion. Figure 4 This is a residual trend diagram of the inner loop of sparse inversion in an embodiment of the present invention. The diagram shows that the residual of the inner loop iteration has an obvious downward trend. After 5 iterations, the residual has dropped by 95%. Generally, the requirement can be met after 2-3 iterations. Figure 5 This is a trend diagram of the L0 norm of the sparse inversion outer loop in the embodiment of the present invention. The smaller the L0 norm of the model, the sparser the model, indicating a higher resolution. This figure shows that the outer loop iteration significantly improves the resolution of the model.

[0096] S104, Scattering Wave Imaging

[0097] Computational scattered wave imaging:

[0098]

[0099] Where, I(iz) is the scattered wave imaging result; iz is the depth domain sampling number, iz = 0, 1, ..., nz-1; nz is the total number of depth domain sampling; dz is the depth domain sampling interval (km); iθ is the sampling number of the illumination inclination, iθ = 0, 1, ..., nθ-1, nθ is the total number of illumination inclination sampling; is the sampling number of the formation dip, is the total number of samples of formation dip; ip is the sampling number of event curvature, ip = 0, 1, ..., np-1; np is the total number of samples of event curvature; p min is the minimum curvature value specified by the user (km); ip1 is the serial number of the lower curvature value used in scattered wave imaging, ip1 = (p a -p min ) / dp; ip2 is the serial number of the upper curvature value used in scattered wave imaging, ip2 = (p b -p min ) / dp;p a and p b The curvature p value range specified by the user requires p a <0 and p b >0, their units are all km; is the kernel function of the above angle trace set. The final output is the imaging trace I(iz)(iz=0,1,…,nz-1).

[0100] S105. All the angle gathers are calculated through S103 and S104, thus completing the scattered wave imaging of the entire work area.

[0101] Compared with the prior art, the embodiment of the present invention combines the generalized Fourier slice theorem and the L1-L2 norm sparse inversion to achieve high-resolution imaging of angular gather scattered waves, thereby improving both the resolution and the computational efficiency.

[0102] Figure 2 is a kernel function surface diagram in the embodiment of the present invention, and the illumination inclination angle θ and the inclination angle used are The range is [-50°, 50°], and the discrete interval is 2°. Figure 2 It shows that the kernel function surface is shaped like butterfly wings.

[0103] Figure 3is the interpolation operator position diagram of the generalized Fourier slice theorem in the embodiment of the present invention, and the interpolation operator position ikp is the depth wave number kz, the illumination inclination angle θ and the inclination angle This figure corresponds to ikz = 154. When the maximum ikp = 65, the interpolation operator positions greater than 65 will be clipped.

[0104] Figure 4 This is a residual trend diagram of the inner loop of sparse inversion in an embodiment of the present invention. The diagram shows that the residual of the inner loop iteration has an obvious downward trend. After 5 iterations, the residual has dropped by 95%. Generally, the requirement can be met after 2-3 iterations.

[0105] Figure 5 This is a trend diagram of the L0 norm of the outer loop of sparse inversion in an embodiment of the present invention. The smaller the L0 norm of the model, the sparser the model, indicating a higher resolution. Figure 5 It shows that the outer loop iteration significantly improves the resolution of the model.

[0106] Figure 6 This is a comparison chart of computational efficiency in embodiments of the present invention. The PCG algorithm for one angle gather takes 210 seconds, while the generalized Fourier slice theorem (GFST)-PCG algorithm of the present invention only takes 4 seconds, and the computational efficiency is improved by 52 times.

[0107] Figure 7 is the angle gather map input in the embodiment of the present invention, the gather has 51 traces, each trace 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 reflected waves and one scattered wave, where the reflected wave event axis is presented as a concave smile curve, and the scattered wave event axis is presented as a horizontal straight line.

[0108] Figure 8 The pz cross-sectional view (left) and Sectional view (right), left The reflected wave and scattered wave are effectively separated, and only the scattered wave phase axis exists in the right figure (p=0km). Fig. 9 It is a comparison diagram of scattered wave imaging effects in the embodiments of the present invention.

[0109] Fig. 9 The left picture shows the images of reflected waves and scattered waves. It can be seen that this method highlights the scattered waves at a depth of 2 km, while the reflected waves at depths of 1 km and 1.5 km are also imaged, but the intensity of the scattered wave imaging is significantly greater than that of the reflected wave. Fig. 9The right figure is a scattered wave imaging diagram. It can be seen that the scattered wave at a depth of 2 km is clearly imaged, while the two reflected waves above it are almost not imaged. From the waveform of the scattered wave, the scattered wave imaging basically maintains the waveform shape. This embodiment shows that the present invention improves the resolution and calculation efficiency of the angular gather scattered wave imaging.

[0110] In summary, Figure 6 The computational efficiency of the embodiment of the present invention is positively demonstrated. Compared with the PCG method, the computational efficiency of the method of the embodiment of the present invention is increased by 52 times. Figure 8 and Fig. 9 The technical effects of the embodiments of the present invention are demonstrated. These figures illustrate the effectiveness and practicability of the method provided by the embodiments of the present invention in achieving high-resolution rapid imaging of angular gather scattered waves.

[0111] In an embodiment of the present invention, a kernel function of an angle gather generated from seismic data is calculated; based on the kernel function of the angle gather, a generalized Fourier slice theorem interpolation operator of the angle gather is calculated; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather; based on the generalized Fourier slice theorem interpolation operator of the angle gather, L1 norm and L2 norm sparse inversion are performed on the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather; based on the three-dimensional model representing the two-dimensional array of the angle gather, a The scattered wave imaging result is obtained, and then the generalized Fourier slice theorem and the L1 norm and L2 norm sparse inversion can be combined through the generalized Fourier slice theorem interpolation operator. There is no need to use the conventional Radon transform under the existing technology to avoid being affected by factors such as aperture, thereby avoiding the problem of resolution reduction and achieving high-resolution imaging of scattered waves of angular gathers. At the same time, there is no need to iteratively solve the Radon transform equation constrained by the model under the existing technology, thus avoiding the problem of large amount of calculation and improving the computational efficiency of scattered wave imaging.

[0112] The present invention also provides a scattered wave imaging device, such as the following embodiment. Since the principle of solving the problem by the device is similar to that of the scattered wave imaging method, the implementation of the device can refer to the implementation of the scattered wave imaging method, and the repeated parts will not be repeated.

[0113] An embodiment of the present invention provides a scattered wave imaging device, such as Fig.11 As shown, it may include:

[0114] Kernel function calculation module 01, used to calculate the kernel function of the angle gather generated by seismic data;

[0115] The interpolation operator calculation module 02 is used to calculate the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather;

[0116] The sparse inversion module 03 is used to perform L1 norm and L2 norm sparse inversion on the angle gather according to the generalized Fourier slice theorem interpolation operator of the angle gather, and obtain a three-dimensional model representing the two-dimensional array of the angle gather;

[0117] The imaging module 04 is used to obtain the scattered wave imaging result according to the three-dimensional model representing the two-dimensional array of angle gathers.

[0118] In one embodiment, the kernel function calculation module is specifically used to calculate the kernel function of the angle gather generated by the seismic data according to the following formula:

[0119]

[0120] in, represents the kernel function of the angular gather; θ is the illumination inclination angle of the angular gather, in radians, θ = θ min +iθ·dθ,θ min is the minimum illumination inclination angle of the angle gather, in radians, dθ is the sampling interval of the illumination inclination angle, in radians; iθ is the sampling sequence number of the illumination inclination angle, iθ=0,1,…,nθ-1, nθ is the total number of samples of the illumination inclination angle; is the formation dip angle of the angle gather, in radians, is the minimum formation dip angle of the angle gather, in radians, is the sampling interval of formation dip, in radians, is the sampling number of the formation dip, is the total number of samples of formation dip.

[0121] In one embodiment, the interpolation operator calculation module is specifically used to calculate the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather, including:

[0122] According to the kernel function of the angle gather, the wave number corresponding to the curvature in the angle gather is discretized;

[0123] According to the wave number corresponding to the curvature in the discretized angle gather, the sequence number of the wave number corresponding to the depth in the angle gather is calculated;

[0124] Calculate the linear interpolation operator coefficients according to the sequence number of the wave number corresponding to the depth in the angle gather;

[0125] Based on the linear interpolation operator coefficients, the generalized Fourier slice theorem interpolation operator of the angle gather is calculated.

[0126] In one embodiment, the interpolation operator calculation module is specifically used to calculate the generalized Fourier slice theorem interpolation operator of the angle gather according to the linear interpolation operator coefficients according to the following formula:

[0127]

[0128] Where S represents the generalized Fourier slice theorem interpolation operator of the angular gather, dimensionless; δ is the coefficient of the linear interpolation operator, δ = (k p -fkp) / dkp-ikp;k p It is the wave number corresponding to the event axis curvature in the angular channel gather, with the unit of 1 / km, referred to as the curvature wave number; is the kernel function of the angle gather, dimensionless; k z k is the wave number corresponding to the depth of the angular channel concentration, the unit is 1 / km, referred to as the depth wave number; z =fkz+ikz·dkz; fkz is the depth wave number k z The minimum value of the depth is expressed in 1 / km; dkz is the depth wave number k z The sampling interval is in 1 / km; ikz is the depth wave number k z The sampling sequence number is ikz = 0, 1, ..., nkz-1; nkz is the depth wave number k z The total number of samples; ikp is the curvature wave number k p The sampling number, represents the integer rounded downward; fkp is the curvature wave number k p The minimum value of the unit is 1 / km, dkp is the curvature wave number k p The sampling interval is 1 / km, ikp = 0, 1, ..., nkp-1, nkp is the curvature wave number k p The total number of samples.

[0129] In one embodiment, the sparse inversion module is specifically used to:

[0130] According to the generalized Fourier slice theorem interpolation operator of the angular gather, the L1 norm and L2 norm optimization problem expressions are established;

[0131] Convert the L1 norm and L2 norm optimization problem expressions into L1 norm and L2 norm linear equations;

[0132] The conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

[0133] In one embodiment, the sparse inversion module is specifically used to establish L1 norm and L2 norm optimization problem expressions according to the following formula:

[0134]

[0135] Where d is a vector of two-dimensional arrays in the z-θ domain representing the angle gather; m is a vector of two-dimensional arrays d representing the angle gather. A vector of arrays of domain 3D models; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; W m is the weight matrix of the three-dimensional model; represents the L2 norm; represents the L1 norm; k is the preset number of iterations in the process of solving the L1 norm and L2 norm optimization problem expressions, m (k) represents the three-dimensional model m obtained in the kth iteration, m (k+1) represents the three-dimensional model m obtained in the k+1th iteration; represents the three-dimensional model m that minimizes the value of {·}; is the sampling number of the formation dip, is the total number of samples of formation dip.

[0136] In one embodiment, the sparse inversion module is specifically used to convert the L1 norm and L2 norm optimization problem expressions into L1 norm and L2 norm linear equations according to the following formula:

[0137]

[0138] Where A is the cross-correlation matrix of the generalized Fourier slice theorem interpolation operator S, b is the vector obtained by applying the interpolation operator S to d, X is the vector composed of the three-dimensional model m, For the three-dimensional model Quantity; is the interpolation operator S of the generalized Fourier slice theorem The matrix consists of components, that is, where i = 0, 1, ..., nθ, j = 0, 1, ..., nkp, and is the generalized Fourier slice theorem interpolation operator The superscript T indicates the matrix transpose. For the matrix The transpose of W m is the weight matrix of the three-dimensional model. The weight matrix is ​​a diagonal matrix. Except for the diagonal elements, all other elements are zero. The i-th diagonal element of the weight matrix is is the mean value of the three-dimensional model m; m i is the value of the i-th element of the three-dimensional model m; Among them, nz, np and are the depth z, event curvature p and formation dip of the three-dimensional model m respectively. Total number of upsamples; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; F -1 is the depth wave number k in the three-dimensional model in the wave number domain z Coordinates and curvature wave number k p Operator that performs the two-dimensional inverse Fourier transform of coordinates.

[0139] In one embodiment, the sparse inversion module is specifically used to solve the L1 norm and L2 norm linear equations using a conjugate gradient algorithm to obtain a three-dimensional model representing the two-dimensional array of angle gathers, including:

[0140] According to the residual and L0 norm between the input angle gather and the pre-calculated angle gather, the conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

[0141] In one embodiment, the imaging module is specifically used to obtain the scattered wave imaging result according to the following formula based on the three-dimensional model representing the two-dimensional array of angle gathers:

[0142]

[0143] Wherein, I(iz) is the scattered wave imaging result; iz is the depth domain sampling number, iz = 0, 1, ..., nz-1; nz is the total number of depth domain sampling; dz is the depth domain sampling interval, the unit is km; iθ is the sampling number of the illumination inclination, iθ = 0, 1, ..., nθ-1, nθ is the total number of illumination inclination sampling; is the sampling number of the formation dip, is the total number of samples of formation dip; ip is the sampling number of event curvature, ip = 0, 1, ..., np-1; np is the total number of samples of event curvature; p min is the minimum curvature value specified by the user, in km; ip1 is the serial number of the lower curvature value used in scattered wave imaging, ip1 = (p a-p min ) / dp; ip2 is the serial number of the upper curvature value used in scattered wave imaging, ip2 = (p b -p min ) / dp;p a and p b is the preset curvature p value range, p a <0 and p b >0, all units are km; The kernel function that describes the angle gather, I(iz) represents the obtained scattered wave imaging result, iz=0,1,…,nz-1.

[0144] 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 scattered wave imaging method when executing the computer program.

[0145] 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 scattered wave imaging method.

[0146] In an embodiment of the present invention, a kernel function of an angle gather generated by seismic data is calculated; based on the kernel function of the angle gather, a generalized Fourier slice theorem interpolation operator of the angle gather is calculated; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather; based on the generalized Fourier slice theorem interpolation operator of the angle gather, L1 norm and L2 norm sparse inversion are performed on the angle gather to obtain a three-dimensional model representing the two-dimensional array of the angle gather; based on the three-dimensional model representing the two-dimensional array of the angle gather , and the scattered wave imaging result is obtained. Thus, the generalized Fourier slice theorem and the L1 norm and L2 norm sparse inversion can be combined through the generalized Fourier slice theorem interpolation operator. There is no need to use the conventional Radon transform under the existing technology to avoid being affected by factors such as aperture, thereby avoiding the problem of resolution reduction and achieving high-resolution imaging of scattered waves of angular gathers. At the same time, there is no need to iteratively solve the Radon transform equation constrained by the model under the existing technology, thus avoiding the problem of large amount of calculation and improving the computational efficiency of scattered wave imaging.

[0147] 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.

[0148] 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.

[0149] 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.

[0150] 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.

[0151] 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 scattered wave imaging method, characterized in that: include: Compute kernel functions for angle gathers generated from seismic data; According to the kernel function of the angle gather, a generalized Fourier slice theorem interpolation operator of the angle gather is calculated; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather; According to the generalized Fourier slice theorem interpolation operator of the angle gather, the L1 norm and L2 norm sparse inversion of the angle gather is performed to obtain a three-dimensional model representing the two-dimensional array of the angle gather; The scattered wave imaging result is obtained based on the three-dimensional model representing the two-dimensional array of angle gathers; According to the kernel function of the angle gather, the generalized Fourier slice theorem interpolation operator of the angle gather is calculated, including: According to the kernel function of the angle gather, the wave number corresponding to the curvature in the angle gather is discretized; According to the wave number corresponding to the curvature in the discretized angle gather, the sequence number of the wave number corresponding to the depth in the angle gather is calculated; Calculate the linear interpolation operator coefficients according to the sequence number of the wave number corresponding to the depth in the angle gather; According to the linear interpolation operator coefficients, the generalized Fourier slice theorem interpolation operator of the angle gather is calculated; The generalized Fourier slice theorem interpolation operator for the angle gather is calculated from the linear interpolation operator coefficients according to the following formula: Where S represents the generalized Fourier slice theorem interpolation operator of the angular gather, dimensionless; δ is the coefficient of the linear interpolation operator, δ = (k p -fkp) / dkp-ikp;k p is the wave number corresponding to the event curvature p in the angular channel gather, with the unit of 1 / km, referred to as the curvature wave number; is the kernel function of the angle gather, dimensionless; k z k is the wave number corresponding to the depth of the angular channel concentration, the unit is 1 / km, referred to as the depth wave number; z =fkz+ikz·dkz; fkz is the depth wave number k z The minimum value of the depth is expressed in 1 / km; dkz is the depth wave number k z The sampling interval is in 1 / km; ikz is the depth wave number k z The sampling sequence number is ikz = 0, 1, ..., nkz-1; nkz is the depth wave number k z The total number of samples; ikp is the curvature wave number k p The sampling number, represents the integer rounded downward; fkp is the curvature wave number k p The minimum value of the unit is 1 / km, dkp is the curvature wave number k p The sampling interval is in 1 / km, ikp = 0, 1, ..., nkp-1, nkp is the curvature wave number k p The total number of samples of the illumination inclination angle θ; iθ is the sampling number of the illumination inclination angle θ, iθ=0,1,…,nθ-1, nθ is the total number of samples of the illumination inclination angle θ; is the formation dip The sampling number, is the formation dip The total number of samples.

2. The method according to claim 1, characterized in that The kernel function of the angle gather generated from seismic data is calculated as follows: in, represents the kernel function of the angular gather; θ is the illumination inclination angle of the angular gather, in radians, θ = θ min +iθ·dθ,θ min is the minimum illumination inclination angle of the angle gather, in radians, dθ is the sampling interval of the illumination inclination angle, in radians; iθ is the sampling sequence number of the illumination inclination angle, iθ = 0, 1, ..., nθ-1, nθ is the total number of samples of the illumination inclination angle; is the formation dip angle of the angle gather, in radians, is the minimum formation dip angle of the angle gather, in radians, is the sampling interval of formation dip, in radians, is the sampling number of the formation dip, is the total number of samples of formation dip, 3. The method according to claim 1, characterized in that According to the generalized Fourier slice theorem interpolation operator of the angle gather, the L1 norm and L2 norm sparse inversion of the angle gather is performed to obtain a three-dimensional model representing the two-dimensional array of the angle gather, including: According to the generalized Fourier slice theorem interpolation operator of the angular gather, the L1 norm and L2 norm optimization problem expressions are established; Convert the L1 norm and L2 norm optimization problem expressions into L1 norm and L2 norm linear equations; The conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

4. The method according to claim 3, characterized in that The L1 norm and L2 norm optimization problem expressions are established as follows: Where d is a vector of two-dimensional arrays in the z-θ domain representing the angle gather; m is a vector representing the two-dimensional array d of the angle gather. A vector of arrays of domain 3D models; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; Wm is the weight matrix of the three-dimensional model; represents the L2 norm; represents the L1 norm; k is the preset number of iterations in the process of solving the L1 norm and L2 norm optimization problem expressions, m (k) represents the three-dimensional model m obtained in the kth iteration, m (k+1) represents the three-dimensional model m obtained in the k+1th iteration; represents the three-dimensional model m that minimizes the value of {·}; is the sampling number of the formation dip, is the total number of samples of formation dip.

5. The method according to claim 4, characterized in that The L1 norm and L2 norm optimization problem expressions are converted into L1 norm and L2 norm linear equations according to the following formula: Where A is the cross-correlation matrix of the generalized Fourier slice theorem interpolation operator S, b is the vector obtained by applying the interpolation operator S to d, X is the vector composed of the three-dimensional model m, For the three-dimensional model Quantity; is the interpolation operator S of the generalized Fourier slice theorem The matrix consists of components, that is, where i = 0, 1, ..., nθ, j = 0, 1, ..., nkp, and is the generalized Fourier slice theorem interpolation operator The superscript T indicates the matrix transpose. For the matrix The transpose of W m is the weight matrix of the three-dimensional model. The weight matrix is ​​a diagonal matrix. Except for the diagonal elements, all other elements are zero. The i-th diagonal element of the weight matrix is is the mean value of the three-dimensional model m; m i is the value of the i-th element of the three-dimensional model m; Among them, nz, np and are the depth z, event curvature p and formation dip of the three-dimensional model m respectively. The total number of upsampling; F is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; F -1 is the depth wave number k in the three-dimensional model in the wave number domain z Coordinates and curvature wave number k p Operator that performs the two-dimensional inverse Fourier transform of coordinates.

6. The method according to claim 3, characterized in that The conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers, including: According to the residual and L0 norm between the input angle gather and the pre-calculated angle gather, the conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

7. The method according to claim 1, characterized in that The scattered wave imaging result is obtained according to the following formula based on the three-dimensional model representing the two-dimensional array of angle gathers: Wherein, I(iz) is the scattered wave imaging result; iz is the depth domain sampling number, iz = 0, 1, ..., nz-1; nz is the total number of depth domain sampling; dz is the depth domain sampling interval, the unit is km; iθ is the sampling number of the illumination inclination, iθ = 0, 1, ..., nθ-1, nθ is the total number of illumination inclination sampling; is the sampling number of the formation dip, is the total number of samples of formation dip; ip is the sampling number of event curvature p, ip = 0, 1, ..., np-1; np is the total number of samples of event curvature p; p min is the minimum curvature value specified by the user, in km; ip1 is the serial number of the lower curvature value used in scattered wave imaging, ip1 = (p a -p min ) / dp; ip2 is the serial number of the upper curvature value used in scattered wave imaging, ip2 = (p b -p min ) / dp;p a and p b is the preset curvature p value range, p a <0 and p b >0, the unit is km; The kernel function that describes the angle gather, I(iz) represents the obtained scattered wave imaging result, iz = 0, 1, ..., nz-1.

8. A scattered wave imaging device, characterized in that: include: A kernel function calculation module is used to calculate the kernel function of the angle gather generated by seismic data; An interpolation operator calculation module is used to calculate a generalized Fourier slice theorem interpolation operator of the angle gather according to a kernel function of the angle gather; the generalized Fourier slice theorem interpolation operator is an array, and the parameters of the array include the illumination inclination angle of the angle gather and the wave number corresponding to the depth in the angle gather; A sparse inversion module is used to perform L1 norm and L2 norm sparse inversion on the angular gathers according to the generalized Fourier slice theorem interpolation operator of the angular gathers to obtain a three-dimensional model representing the two-dimensional array of the angular gathers; An imaging module, used to obtain scattered wave imaging results based on a three-dimensional model representing a two-dimensional array of angular gathers; The interpolation operator calculation module is specifically used to calculate the generalized Fourier slice theorem interpolation operator of the angle gather according to the kernel function of the angle gather, including: According to the kernel function of the angle gather, the wave number corresponding to the curvature in the angle gather is discretized; According to the wave number corresponding to the curvature in the discretized angle gather, the sequence number of the wave number corresponding to the depth in the angle gather is calculated; Calculate the linear interpolation operator coefficients according to the sequence number of the wave number corresponding to the depth in the angle gather; According to the linear interpolation operator coefficients, the generalized Fourier slice theorem interpolation operator of the angle gather is calculated; The interpolation operator calculation module is specifically used to calculate the generalized Fourier slice theorem interpolation operator of the angle gather according to the linear interpolation operator coefficients according to the following formula: Where S represents the generalized Fourier slice theorem interpolation operator of the angular gather, dimensionless; δ is the coefficient of the linear interpolation operator, δ = (k p -fkp) / dkp-ikp;k p is the wave number corresponding to the event curvature p in the angular channel gather, with the unit of 1 / km, referred to as the curvature wave number; is the kernel function of the angle gather, dimensionless; k z k is the wave number corresponding to the depth of the angular channel concentration, the unit is 1 / km, referred to as the depth wave number; z =fkz+ikz·dkz; fkz is the depth wave number k z The minimum value of the depth is expressed in 1 / km; dkz is the depth wave number k z The sampling interval is in 1 / km; ikz is the depth wave number k z The sampling sequence number is ikz = 0, 1, ..., nkz-1; nkz is the depth wave number k z The total number of samples; ikp is the curvature wave number k p The sampling number, represents the integer rounded downward; fkp is the curvature wave number k p The minimum value of the unit is 1 / km, dkp is the curvature wave number k p The sampling interval is in 1 / km, ikp = 0, 1, ..., nkp-1, nkp is the curvature wave number k p The total number of samples of the illumination inclination angle θ; iθ is the sampling number of the illumination inclination angle θ, iθ=0,1,…,nθ-1, nθ is the total number of samples of the illumination inclination angle θ; is the formation dip The sampling number, is the formation dip The total number of samples.

9. The device according to claim 8, characterized in that The kernel function calculation module is specifically used to calculate the kernel function of the angle gather generated by seismic data according to the following formula: in, represents the kernel function of the angular gather; θ is the illumination inclination angle of the angular gather, in radians, θ = θ min +iθ·dθ,θ min is the minimum illumination inclination angle of the angle gather, in radians, dθ is the sampling interval of the illumination inclination angle, in radians; iθ is the sampling sequence number of the illumination inclination angle, iθ=0,1,…,nθ-1, nθ is the total number of samples of the illumination inclination angle; is the formation dip angle of the angle gather, in radians, is the minimum formation dip angle of the angle gather, in radians, is the sampling interval of formation dip, in radians, is the sampling number of the formation dip, is the total number of samples of formation dip, 10. The device according to claim 8, characterized in that Sparse inversion module, specifically used for: According to the generalized Fourier slice theorem interpolation operator of the angular gather, the L1 norm and L2 norm optimization problem expressions are established; Convert the L1 norm and L2 norm optimization problem expressions into L1 norm and L2 norm linear equations; The conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

11. The device according to claim 10, characterized in that The sparse inversion module is specifically used to establish the L1 norm and L2 norm optimization problem expressions according to the following formula: Where d is a vector of two-dimensional arrays in the z-θ domain representing the angle gather; m is a vector representing the two-dimensional array d of the angle gather. A vector of arrays of domain 3D models; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; Wm is the weight matrix of the three-dimensional model; represents the L2 norm; represents the L1 norm; k is the preset number of iterations in the process of solving the L1 norm and L2 norm optimization problem expressions, m (k) represents the three-dimensional model m obtained in the kth iteration, m (k+1) represents the three-dimensional model m obtained in the k+1th iteration; represents the three-dimensional model m that minimizes the value of {·}; is the sampling number of the formation dip, is the total number of samples of formation dip.

12. The device according to claim 11, characterized in that Sparse inversion module, specifically used for: Where A is the cross-correlation matrix of the generalized Fourier slice theorem interpolation operator S, b is the vector obtained by applying the interpolation operator S to d, X is the vector composed of the three-dimensional model m, For the three-dimensional model Quantity; is the interpolation operator S of the generalized Fourier slice theorem The matrix consists of components, that is, where i = 0, 1, ..., nθ, j = 0, 1, ..., nkp, and is the generalized Fourier slice theorem interpolation operator The superscript T indicates the matrix transpose. For the matrix The transpose of W m is the weight matrix of the three-dimensional model. The weight matrix is ​​a diagonal matrix. Except for the diagonal elements, all other elements are zero. The i-th diagonal element of the weight matrix is is the mean value of the three-dimensional model m; m i is the value of the i-th element of the three-dimensional model m; Among them, nz, np and are the depth z, event curvature p and formation dip of the three-dimensional model m respectively. Total number of upsamples; is the operator for Fourier transform of the depth z coordinate in the angular gather; F is the operator for two-dimensional Fourier transform of the depth z coordinate and the curvature p coordinate in the angular gather; F -1 is the depth wave number k in the three-dimensional model in the wave number domain z Coordinates and curvature wave number k p Operator that performs the two-dimensional inverse Fourier transform of coordinates.

13. The device according to claim 10, characterized in that The sparse inversion module is specifically used to solve the L1 norm and L2 norm linear equations using the conjugate gradient algorithm to obtain a three-dimensional model representing the two-dimensional array of angle gathers, including: According to the residual and L0 norm between the input angle gather and the pre-calculated angle gather, the conjugate gradient algorithm is used to solve the L1 norm and L2 norm linear equations to obtain a three-dimensional model representing the two-dimensional array of angle gathers.

14. The device according to claim 8, characterized in that The imaging module is specifically used to obtain the scattered wave imaging result according to the following formula based on the three-dimensional model representing the two-dimensional array of angle gathers: Wherein, I(iz) is the scattered wave imaging result; iz is the depth domain sampling number, iz = 0, 1, ..., nz-1; nz is the total number of depth domain sampling; dz is the depth domain sampling interval, the unit is km; iθ is the sampling number of the illumination inclination angle θ, iθ = 0, 1, ..., nθ-1, nθ is the total number of sampling of the illumination inclination angle θ; is the formation dip The sampling number, is the formation dip The total number of samples; ip is the sampling number of the event axis curvature p, ip = 0, 1, ..., np-1; np is the total number of samples of the event axis curvature p; p min is the minimum curvature value specified by the user, in km; ip1 is the serial number of the lower curvature value used in scattered wave imaging, ip1 = (p a -p min ) / dp; ip2 is the serial number of the upper curvature value used in scattered wave imaging, ip2 = (p b -p min ) / dp;p a and p b is the preset curvature p value range, p a <0 and p b >0, the unit is km; The kernel function that describes the angle gather, I(iz) represents the obtained scattered wave imaging result, iz = 0, 1, ..., nz-1.

15. 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 7 is implemented.

16. 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 7.

Images