Discrete method, system, terminal and medium for angle-domain generalized Radon transform

Abstract

The present application relates to a discretization method, system, terminal, and storage medium for generalized Radon transform in the angle domain. The method includes: reading seismic trace data, setting the coordinates of imaging points of the seismic trace data, and reading the travel time, amplitude, and direction table of each imaging point coordinate; calculating the scattering angle and discrete interval; calculating the sampling point θ0 closest to the scattering angle and its partition boundary #imgabs0# and calculating the distance from the scattering angle to the partition boundary θ bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the discrete element splitting algorithm to split the partition boundary θ bound The nearby discrete cells are split into two cells, and the areas of the two cells after the split are calculated according to the discrete interval; the two cells after the split are used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ. This application can obtain a continuous and smooth angle domain inversion result.

Description

Technical Field

[0001] The present application belongs to the technical field of angle domain information extraction, and in particular relates to a discrete method, system, terminal and storage medium for angle domain generalized Radon transform. Background Art

[0002] The inventors proposed and established the angle-domain generalized Radon transform (AD-GRT) inversion method (Li Xuelei et al., 2020), which can achieve multi-parameter inversion of acoustic waves and elastic media in a non-iterative manner. The theoretical framework of the angle-domain generalized Radon transform is:

[0003] Based on the Born approximation, the wave field of the acoustic medium can be expressed as a generalized Radon transform (GRT) form (frequency domain):

[0004] p′(r,s,ω)=ω 2 ∫dx1dx3A s A r exp[iω(φ s +φ r )]κ 0 F(ω)f(x,θ) (1)

[0005] Among them, p′(r,s,ω) is the acoustic pressure scattering wave field, A s and φ s Represent the amplitude and travel time at the source end, A r and φ r are the amplitude and travel time at the receiving end, respectively, and F(ω) is the source wavelet. θ = θ(s, x, r) is the scattering angle, which is the angle between the s-to-x-ray and the r-to-x-ray at x, and satisfies:

[0006] θ=α s -α r (2)

[0007] Among them, α s =α s (x) is the direction angle from s to the x-ray at x, α r =α r (x) is the direction angle from r to the x-ray at x.

[0008] f(x,θ) is an angle domain model:

[0009]

[0010] σ 0 (x) and κ 0 (x) are the background models of the density inverse and the compressibility, respectively, and σ′ and κ′ are the corresponding perturbation parameters.

[0011] The GRT integral transform is equivalent to the integral transform of the angle domain model f(x,θ) to the acoustic pressure scattering wave field p′(r,s,t), where s and r are only distributed on the boundary line of the 2D disturbance region. In theory, both f(x,θ) and p′(r,s,t) have three-dimensional distributions, so the corresponding inverse transform from p′(r,s,t) to f(x,θ) can also be constructed. By establishing the corresponding angle domain generalized inverse Radon transform (AD-GRT), the inverse integral transform from p′(r,s,t) to f(x,θ) is supported. The inverse integral transform expression is as follows:

[0012]

[0013] Among them, J s (y) and J r (y) is the Jacobian determinant and can be expressed as and α s (y) and α r (y) is the direction angle between the mappable source s and the receiving point r. s (y) and J r (y) can be represented by the geometric diffusion from y to s and r:

[0014]

[0015] However, δ(θ-θ0) in equation (4) cannot be directly used in the calculation and needs to be replaced by other approximate interpolation or summation algorithms. In addition, the discrete sampling points of θ0 do not directly correspond to the discrete distribution of the source s and the receiving point r, which brings certain difficulties to the accurate numerical calculation of AD-GRT. Therefore, it is necessary to design an accurate discrete algorithm for AD-GRT based on the conventional angle domain discrete algorithm. The inverse integral transform expression (4) can be approximately expressed as:

[0016]

[0017] Where Δθ is the sampling interval of the scattering angle θ0. In the traditional angle domain discrete algorithm, if If the scattering angle θ corresponding to the imaging point y falls within the interval [θ0-Δθ / 2,θ0+Δθ / 2), then the integral term in Equation (4) is accumulated at the θ0 sampling point of f(y,θ0) and finally multiplied by 1 / Δθ to achieve averaging. However, due to the discrete integer nature of the algorithm, this conventional algorithm results in very uneven calculations of f(y,θ0) for different θ0 values, prone to sawtooth-like fluctuations along the θ0 distribution. This fluctuation is particularly severe when the source and receiving points are sparsely distributed, or when the θ0 values are densely sampled, making it difficult to extract an accurate distribution relationship for f(y,θ0), rendering the angle-domain inversion results invalid. Even if some smoothing factors can mitigate the fluctuations, significant fluctuations still exist, and the smoothing factors can result in additional information loss. Summary of the Invention

[0018] The present application provides a discrete method, system, terminal and storage medium for angle-domain generalized Radon transform, aiming to solve at least one of the above-mentioned technical problems in the prior art to a certain extent.

[0019] In order to solve the above problems, this application provides the following technical solutions:

[0020] A discretization method for angle-domain generalized Radon transform, comprising:

[0021] Read seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate;

[0022] Calculate the scattering angle and discrete interval based on the travel time, amplitude and direction table corresponding to each imaging point coordinate;

[0023] Calculate the nearest sampling point θ0 and its partition boundary of the scattering angle Δθ is the scattering angle interval, and the scattering angle to the partition boundary θ is calculated bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the discrete element splitting algorithm to split the partition boundary θ bound The nearby discrete cells are split into two cells, and the areas of the two cells after the split are calculated according to the discrete interval;

[0024] The two split units are used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ respectively.

[0025] The technical solution adopted in the embodiment of the present application further includes: before reading the seismic trace data, the following steps are further included:

[0026] Set the discrete interval Δs of the shot points, the discrete interval Δr of the receiving points, and the scattering angle interval Δθ;

[0027] The shot point coordinates and the corresponding ray travel time, amplitude and direction table are set, and the receiving point coordinates and the corresponding ray travel time, amplitude and direction table are set.

[0028] The technical solution adopted in the embodiment of the present application also includes: the travel time, amplitude and direction table for reading the coordinates of each imaging point is specifically:

[0029] The direction table includes the direction angle α s and α r , where α s =α s (x) is the direction angle from s to the x-ray at x, α r =α r (x) is the direction angle from r to the x-ray at x.

[0030] The technical solution adopted in the embodiment of the present application also includes: the calculation of the scattering angle and discrete interval according to the travel time, amplitude and direction table corresponding to the coordinates of each imaging point is specifically as follows:

[0031] The scattering angle θ mid =α s -α r ;

[0032] The discrete interval d s =J s Δs and d r =J r Δr, where J s and J r is the Jacobian determinant, Δs is the discrete interval of the shot points, and Δr is the discrete interval of the receiving points.

[0033] The technical solution adopted in the embodiment of the present application also includes: the discrete unit splitting algorithm includes a rectangular splitting algorithm and a parallel splitting algorithm.

[0034] The technical solution adopted in the embodiment of the present application also includes: the rectangle splitting algorithm is specifically:

[0035] Assume d s =Δα s , d r =Δα r , l=|θ bound -θ mid |, Δα s and Δα r are the ray direction intervals between rays s and r propagating to point y, Δα s ≈J s Δs,Δα r ≈J r Δr;

[0036] if Then it represents the partition boundary θ bound Passing through the discrete unit, solve the two split units S1 and S2:

[0037] If d s ≥d r :

[0038]

[0039] If d s <d r :

[0040]

[0041] S2=d s d r -S1

[0042] The technical solution adopted in the embodiment of the present application also includes: the parallel splitting algorithm is specifically:

[0043] if represents the partition boundary θ bound Passing through the discrete element, solve S1 and S2:

[0044]

[0045]

[0046] S1+S2=d s d r .

[0047] Another technical solution adopted by the embodiment of the present application is: a discrete system of angle-domain generalized Radon transform, comprising:

[0048] Data reading module: used to read seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate;

[0049] Scattering angle calculation module: used to calculate the scattering angle and discrete interval based on the travel time, amplitude and direction table corresponding to the coordinates of each imaging point;

[0050] Discrete unit splitting module: used to calculate the sampling point θ0 closest to the scattering angle and its partition boundary Δθ is the scattering angle interval, and the scattering angle to the partition boundary θ is calculated bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the discrete element splitting algorithm to split the partition boundary θ boundThe nearby discrete cells are split into two cells, and the areas of the two cells after the split are calculated according to the discrete interval;

[0051] Imaging module: used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ using the two split units respectively.

[0052] Another technical solution adopted by the embodiment of the present application is: a terminal, the terminal including a processor and a memory coupled to the processor, wherein:

[0053] The memory stores program instructions for implementing a discrete method of the angle domain generalized Radon transform;

[0054] The processor is configured to execute the program instructions stored in the memory to control discretization of angle-domain generalized Radon transform.

[0055] Another technical solution adopted by the embodiment of the present application is: a storage medium storing program instructions executable by a processor, wherein the program instructions are used to execute the discrete method of the angle domain generalized Radon transform.

[0056] Compared with the prior art, the beneficial effects of the embodiments of the present application are as follows: the discrete method, system, terminal and storage medium of the angle domain generalized Radon transform of the embodiments of the present application combine the theoretical framework and physical meaning of AD-GRT transformation to design two angle domain discrete difference algorithms that conform to the AD-GRT transformation theory. By quantitatively splitting the discrete unit area, the discrete unit after discretization is rationally split into two units, and the accurate calculation of each discrete unit is achieved. While ensuring the efficiency and convenience of the algorithm, the smooth continuity of the angle domain amplitude is ensured, and the amplitude oscillation problem of the conventional angle domain discrete algorithm is solved. The embodiments of the present application can obtain continuous and smooth angle domain inversion results, and realize efficient and stable angle domain information extraction. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 is a flowchart of a discrete method of angle-domain generalized Radon transform according to an embodiment of the present application;

[0058] Figure 2 Schematic diagram of the distribution of rectangular discrete units and scattering angle contour lines according to an embodiment of the present application;

[0059] Figure 3 Schematic diagram of the distribution of parallelogram discrete units and scattering angle contour lines according to an embodiment of the present application;

[0060] Figure 4 a is a schematic diagram of a horizontal single interface model. Figure 4 (b) is a schematic diagram of the synthetic single-shot seismic record;

[0061] Figure 5 Schematic diagram of the inversion results of the angle domain model f(x,θ) solved using the AD-GRT inversion method;

[0062] Figure 6 Schematic diagram of the angle domain inversion value distribution on the horizontal line z = 1000m, where (a), (b), and (c) are the true values and inverted values using the traditional discretization method, rectangular splitting algorithm, and parallel splitting algorithm, respectively;

[0063] Figure 7 Schematic diagram of the discrete system structure of the angle domain generalized Radon transform in an embodiment of the present application;

[0064] Figure 8 This is a schematic diagram of the terminal structure of an embodiment of the present application;

[0065] Figure 9 A schematic diagram of the structure of the storage medium of an embodiment of the present application. DETAILED DESCRIPTION

[0066] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0067] See also Figure 1 , is a flow chart of a discrete method for the angle domain generalized Radon transform according to an embodiment of the present application. The discrete method for the angle domain generalized Radon transform according to an embodiment of the present application comprises the following steps:

[0068] S1: Set the operating environment parameters such as the discrete interval Δs of the shot points, the discrete interval Δr of the receiving points, and the scattering angle interval Δθ;

[0069] S2: Set the shot point coordinates and the corresponding ray travel time, amplitude and direction table, and set the receiving point coordinates and the corresponding ray travel time, amplitude and direction table;

[0070] Among them, the direction table includes the direction angle α s and α r , where α s =α s (x) is the direction angle from s to the x-ray at x, α r =α r (x) is the direction angle from r to the x-ray at x.

[0071] S3: read seismic trace data;

[0072] S4: Set the imaging point coordinates of the seismic trace data and read the travel time, amplitude and direction table corresponding to each imaging point coordinate;

[0073] S5: Calculate the scattering angle θ based on the travel time, amplitude and direction table corresponding to each imaging point coordinate mid =α s -α r and its discrete interval d s =J s Δs and d r =J r Δr, where J s and J r is the Jacobian determinant;

[0074] S6: Calculate the scattering angle θ mid The nearest scattering angle sampling point θ0 and its partition boundary And calculate the scattering angle θ mid To the partition boundary θ bound The distance l=|θ bound -θ mid |;

[0075] S7: According to the distance l=|θ bound -θ mid | Determine the partition boundary θ bound Whether to split the scattering angle θ mid If the discrete surface is, the partition boundary θ is split by the discrete element splitting algorithm. bound The nearby discrete unit is split into two units S1 and S2, and the area size of the two units after the split is calculated;

[0076] In this step, in order to solve the discretization problem of the traditional angle domain algorithm, the embodiment of the present application further analyzes the integral summation related to the scattering angle θ and designs a reasonable numerical solution. The discrete infinitesimal element satisfies the relationship:

[0077] Δα s ≈J s Δs,Δα r ≈J r Δr (7)

[0078] Where Δs and Δr are the discrete intervals between the shot point and the receiver point, respectively, s and Δα r The ray direction intervals of rays s and r propagating to y are respectively, which can be obtained by the Jacobi coefficient J s and J r Calculated. Therefore, a discrete sampling point (s, r) corresponds to a discrete unit in y, and the area of the discrete unit is Δα s ×Δα rWhen the sampling point is very close to the partition boundary θ0±Δθ / 2, a part of the discrete unit should belong to the partition near the sampling point θ0. Based on the above, the embodiment of the present application sets the partition boundary θ bound The nearby discrete cells are split into two cells, and the area sizes S1 and S2 of the two cells after the split are solved.

[0079] Specifically, such as Figure 2 The following is a schematic diagram of the distribution of rectangular discrete units and scattering angle contour lines in the embodiment of the present application. The specific rectangular splitting algorithm is: each (α s ,α r ) The discrete unit of the sampling point is represented by a rectangle with an area of Δα s ×Δα r , the sampling point is located in the middle of the length and width, θ b o und =θ0±Δθ / 2 is the boundary of θ0 partition, θ mid is (α s ,α r ) The θ value corresponding to the sampling point. When the partition boundary θ bound When passing through a discrete unit, the discrete unit size is split into two parts, S1 and S2, and replaced by J in formula (4) s J r The ΔsΔr part is added to the inversion results of f(y,θ0) in and around θ0.

[0080] Assume d s =Δα s , d r =Δα r , l=|θ bound -θ mid |; if Then the partition boundary θ bound Passing through this discrete unit, it is necessary to solve the two split units S1 and S2. The solution method is as follows:

[0081] If d s ≥d r :

[0082]

[0083] If d s <d r :

[0084]

[0085] Based on S1, S2 is:

[0086] S2=d s d r -S1. (10)

[0087] The parallel splitting algorithm is as follows: the scattering angle θ satisfies If you set Δα s =0(or α s is a constant), then the integral range of θ in formula (2) can be obtained by r Individual control. Figure 3 As shown in the figure, it is a schematic diagram of the distribution of parallelogram discrete units and scattering angle contour lines. s ,α r ) The discrete unit of the sampling point is represented by a parallelogram (i.e., two sides of which are at α s The other two sides are parallel to the θ isovalue line), and the discrete unit splitting at this time is only related to α r Discrete situations.

[0088] if Then θ bound Passing through the discrete unit, S1 and S2 need to be solved, and S1 and S2 satisfy:

[0089]

[0090]

[0091] Here, S1+S2=d s d r This method can also be used through α s Control discrete units individually and swap all d in the formula r and d s The parallel splitting algorithm is simple and convenient to calculate, and is suitable for future high-dimensional discrete computing.

[0092] S8: Use two units S1 and S2 to calculate the imaging accumulation at the sampling point θ0±Δθ and the scattering angle sampling point θ0 respectively;

[0093] S9: Determine whether the imaging point of the seismic trace data has been cycled through. If not, obtain the next imaging point and re-execute S4 to S8; otherwise, execute S10;

[0094] S10: Determine whether there are any remaining receiving points. If so, set the next receiving point and re-execute S2 to S9; otherwise, execute S11.

[0095] S11: Determine whether there are any remaining shot points. If there are any remaining shot points, set the next shot point and re-execute S2 to S10; otherwise, end.

[0096] Data experiments:

[0097] In order to verify the feasibility and effectiveness of the embodiment of the present application, the following simple model data is used to verify the calculation effect of the embodiment of the present application. Specifically: the horizontal single interface model is as follows: Figure 4 As shown in Figure a, the horizontal interface is at a depth of z = 1000 m, and the acoustic wave velocities above and below the interface are 2000 m / s and 2100 m / s, respectively. The model density is set according to the Gardner relationship ρ = 0.31 × c0.25 (velocity units: m / s, density units: g / cm³). The model grid spacing is 5 m in both the horizontal and vertical directions. Shots and receivers are distributed on the boundary line z = 0 m, with a shot and receiver spacing of 10 m. The maximum offset from shot to receiver is 2000 m. The model seismic records are synthesized using the finite difference method of the acoustic wave equation to verify the validity of the accurate wavefield data inversion results. The background velocity and density models used for AD-GRT are obtained by smoothing the original model. Figure 4 (b) is a schematic diagram of a synthetic single-shot seismic record, showing a single-shot record with a source x = 2000 m. Here, only the main reflection information is shown, and the direct wave information is excluded.

[0098] Figure 5 The inversion results of the angle-domain model f(x,θ) (x = 2000 m) obtained using the AD-GRT inversion method are presented. To eliminate amplitude oscillations caused by filtering, the source wavelet is set to F(ω) = 1, and the wavelet inverse filtering operation is ignored. The main inverted information is distributed near the horizontal line z = 1000 m, and the phase axis has a clear horizontal distribution feature. At the same time, there is significant tilt noise at the edge of the inversion result (see area A), which is related to data truncation at the maximum offset boundary of the seismic record. The record boundary is characterized by data discontinuity, and the discontinuity does not conform to the local approximation assumption of the AD-GRT.

[0099] In order to verify the effectiveness of the discrete unit splitting method in numerical calculation of the embodiment of the present application, the numerical calculation inversion results of the traditional discrete method, the rectangular splitting algorithm and the parallel splitting algorithm are respectively used below. Figure 6 As shown in FIG, a schematic diagram of the angle domain inversion value distribution on the horizontal line z = 1000 m, where (a), (b), and (c) are the true values and inverted values using the traditional discrete method, rectangular splitting algorithm, and parallel splitting algorithm, respectively. Figure 6The horizontal axis direction is represented by the cosine value of the angle to show the linear tilt distribution characteristics of the angle domain model. The true value is the difference (jump value) of the disturbance value above and below the interface z = 1000m, and for the convenience of comparison, the inversion value is multiplied by a constant. Experimental results show that the traditional discrete method is prone to form obvious sawtooth oscillations in the extraction of angle domain information, especially in a small angle range. The rectangular splitting algorithm and the parallel splitting algorithm proposed in the embodiments of the present application can both provide a very continuous and smooth amplitude distribution effect, and the inverted amplitude distribution and the model true value distribution can be highly fitted within the effective angle range.

[0100] Based on the above, the discrete method of the angle domain generalized Radon transform in the embodiment of the present application combines the theoretical framework and physical meaning of the AD-GRT transform to design two angle domain discrete difference algorithms that conform to the AD-GRT transform theory. By quantitatively splitting the discrete unit area, the discrete unit after discretization is reasonably split into two units, and the accurate calculation of each discrete unit is achieved. While ensuring the efficiency and convenience of the algorithm, it ensures the smooth continuity of the angle domain amplitude and solves the amplitude oscillation problem of the conventional angle domain discrete algorithm. The embodiment of the present application can obtain a continuous and smooth angle domain inversion result and realize efficient and stable angle domain information extraction.

[0101] See also Figure 7 , is a schematic diagram of the discrete system structure of the angle domain generalized Radon transform in the embodiment of the present application. The discrete system 40 of the angle domain generalized Radon transform in the embodiment of the present application includes:

[0102] Data reading module 41: used to read seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate;

[0103] Scattering angle calculation module 42: used to calculate the scattering angle and discrete interval according to the travel time, amplitude and direction table corresponding to each imaging point coordinate;

[0104] Discrete unit splitting module 43: used to calculate the nearest sampling point θ0 of the scattering angle and its partition boundary Δθ is the scattering angle interval, and the scattering angle to the partition boundary θ is calculated bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the discrete element splitting algorithm to split the partition boundary θ bound The nearby discrete cells are split into two cells, and the areas of the two cells after the split are calculated based on the discrete interval;

[0105] Imaging module 44 is used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ using the two split units respectively.

[0106] See also Figure 8 , is a schematic diagram of the terminal structure of an embodiment of the present application. The terminal 50 includes a processor 51 and a memory 52 coupled to the processor 51.

[0107] The memory 52 stores program instructions for implementing the above-mentioned discrete method of angle-domain generalized Radon transform.

[0108] The processor 51 is configured to execute program instructions stored in the memory 52 to control the discretization of the angle-domain generalized Radon transform.

[0109] The processor 51 may also be referred to as a CPU (Central Processing Unit). The processor 51 may be an integrated circuit chip having signal processing capabilities. The processor 51 may also be a general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic device, a discrete gate or transistor logic device, or a discrete hardware component. The general-purpose processor may be a microprocessor or any conventional processor.

[0110] See also Figure 9 , which is a structural diagram of the storage medium of an embodiment of the present application. The storage medium of the embodiment of the present application stores a program file 61 that can implement all the above methods, wherein the program file 61 can be stored in the above storage medium in the form of a software product, including a number of instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) or a processor (processor) to execute all or part of the steps of the methods of each embodiment of the present invention. The aforementioned storage medium includes: various media that can store program codes, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, or terminal devices such as a computer, a server, a mobile phone, and a tablet.

[0111] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application is not limited to the embodiments shown herein, but is intended to encompass the broadest scope consistent with the principles and novel features disclosed herein.

Claims

1. A discretization method for generalized Radon transform in angle domain, characterized in that: include: Read seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate; Calculate the scattering angle and discrete interval based on the travel time, amplitude and direction table corresponding to each imaging point coordinate; Calculate the nearest sampling point θ0 and its partition boundary of the scattering angle Δθ is the scattering angle interval, and the scattering angle to the partition boundary θ is calculated bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the rectangular splitting algorithm to split the partition boundary θ bound The nearby discrete unit is split into two units, and the areas S1 and S2 of the two units after the split are calculated according to the discrete interval; The two split units are used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ respectively; The rectangle splitting algorithm is specifically as follows: Let d s = Δα s , d r = Δα r , l = |θ bound - θ mid |, Δα s and Δα r are respectively the intervals of the directions at which s and r propagate to y, Δα s ≈ J s Δs, Δα r ≈ J r Δr; if Then it represents the partition boundary θ bound Passing through the discrete unit, the areas S1 and S2 of the two units after splitting are solved: If d s ≥d r : If d s <d r : S2=d s d r -S1; Before reading the seismic trace data, the method further includes: Set the discrete interval Δs of the shot points, the discrete interval Δr of the receiving points, and the scattering angle interval Δθ; Setting the shot point coordinates and the corresponding ray travel time, amplitude, and direction table, and setting the receiving point coordinates and the corresponding ray travel time, amplitude, and direction table; The travel time, amplitude and direction table for reading the coordinates of each imaging point is specifically: The direction table includes the direction angle α s (x) and α r (x), where α s (x) is the direction angle of the ray from s to x at x, α r (x) is the direction angle of the ray from r to x at x; The calculation of the scattering angle and discrete interval according to the travel time, amplitude and direction table corresponding to each imaging point coordinate is specifically as follows: The scattering angle θ mid =α s (x)-α r (x); The discrete interval d s =J s Δs and d r =J r Δr, where J s and J r is the Jacobian determinant, Δs is the discrete interval of the shot points, and Δr is the discrete interval of the receiving points.

2. A discretization method for generalized Radon transform in angle domain, characterized in that: include: Read seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate; Calculate the scattering angle and discrete interval based on the travel time, amplitude and direction table corresponding to each imaging point coordinate; Calculate the nearest sampling point θ0 and its partition boundary of the scattering angle Δθ is the scattering angle interval, and the scattering angle to the partition boundary θ is calculated bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the parallel splitting algorithm to split the partition boundary θ bound The nearby discrete unit is split into two units, and the areas S1 and S2 of the two units after the split are calculated according to the discrete interval; The two split units are used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ respectively; The parallel splitting algorithm is specifically as follows: if represents the partition boundary θ bound Passing through the discrete unit, the areas S1 and S2 of the two units are solved: S1+S2=d s d r ; Before reading the seismic trace data, the method further includes: Set the discrete interval Δs of the shot points, the discrete interval Δr of the receiving points, and the scattering angle interval Δθ; Setting the shot point coordinates and the corresponding ray travel time, amplitude, and direction table, and setting the receiving point coordinates and the corresponding ray travel time, amplitude, and direction table; The travel time, amplitude and direction table for reading the coordinates of each imaging point is specifically: The direction table includes the direction angle α s (x) and α r (x), where α s (x) is the direction angle of the ray from s to x at x, α r (x) is the direction angle of the ray from r to x at x; The calculation of the scattering angle and discrete interval according to the travel time, amplitude and direction table corresponding to each imaging point coordinate is specifically as follows: The scattering angle θ mid =α s (x)-α r (x); The discrete interval d s =J s Δs and d r =J r Δr, where J s and J r is the Jacobian determinant, Δs is the discrete interval of the shot points, and Δr is the discrete interval of the receiving points.

3. A discrete system for implementing the angle domain generalized Radon transform discretization method of angle domain generalized Radon transform according to claim 1 or 2, characterized in that: include: Data reading module: used to read seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate; Scattering angle calculation module: used to calculate the scattering angle and discrete interval based on the travel time, amplitude and direction table corresponding to the coordinates of each imaging point; Discrete unit splitting module: used to calculate the sampling point θ0 closest to the scattering angle and its partition boundary Δθ is the scattering angle interval, and the scattering angle to the partition boundary θ is calculated bound The distance l is used to determine the partition boundary θ bound Whether to split the scattering angle discrete surface, if so, use the discrete element splitting algorithm to split the partition boundary θ bound Splitting the nearby discrete cells into two cells, and calculating the areas of the two cells after the splitting according to the discrete interval, wherein the discrete cell splitting algorithm includes a rectangular splitting algorithm and a parallel splitting algorithm; Imaging module: used to calculate the imaging accumulation at the sampling point θ0 and θ0±Δθ using the two split units respectively.

4. A terminal, characterized in that: The terminal includes a processor and a memory coupled to the processor, wherein: The memory stores program instructions for implementing the discrete method of the angle domain generalized Radon transform according to any one of claims 1-2; The processor is configured to execute the program instructions stored in the memory to control discretization of angle-domain generalized Radon transform.

5. A storage medium, characterized in that: Program instructions executable by a processor are stored, and the program instructions are used to execute the discretization method of the angle domain generalized Radon transform according to any one of claims 1 to 2.

Images