A method, apparatus, storage medium, and electronic device for seismic profile imaging.

Abstract

This invention discloses a method, apparatus, storage medium, and electronic device for seismic profile imaging. The method includes: performing a Radon transform on each first seismic wavefield at the shot point to obtain a first plane wave component and obtaining the value of a first ray parameter, thereby obtaining the incident ray angle corresponding to the first ray parameter; performing a Radon transform on each second seismic wavefield on one side of the receiver point to obtain a second plane wave component and obtaining the value of each second ray parameter, thereby obtaining the scattered ray angle corresponding to each second ray parameter; obtaining each dip angle based on the incident angle and each scattered angle; obtaining the single-shot imaging function for each dip angle through a single-shot imaging calculation formula; iteratively executing the imaging function for all shot points to obtain dip gathers; and superimposing the dip gathers to obtain seismic profile imaging. This method can efficiently obtain accurate dip gathers, thereby obtaining high-resolution, high signal-to-noise ratio seismic profile imaging.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical exploration, and in particular to a method and device for seismic profile imaging, a storage medium and an electronic device. BACKGROUND

[0002] In the depth domain imaging of seismic data, the calculation of dip angle gathers can be used for imaging discontinuous points and velocity analysis, and the dip angle gather is an important gather.

[0003] In the prior art, a method for calculating a dip angle gather is as follows: the travel time of a shot point to an imaging point and the angle of incidence of a wave field at the imaging point are calculated at each imaging point by ray tracing, and the travel time of the imaging point to a receiver point and the propagation angle of a scattered wave field starting from the imaging point are calculated, so that the corresponding dip angle at the imaging point is calculated from the two angles, and the corresponding migration energy is placed on the dip angle. For each imaging point, different combinations of shot points and receiver points will calculate different dip angles, and the dip angle migration results obtained for the same imaging point form the dip angle gather of the imaging point. This method is the source of many methods for calculating dip angle gathers based on ray migration. In the prior art, the ray tracing is replaced by a travel time table of a shot point to an imaging point calculated by an eikonal equation, and a travel time table of a receiver point to the imaging point, and the angles of the incident ray and the scattered ray are calculated from the travel time gradient, so as to calculate the dip angle and obtain the dip angle gather. The common problem of this type of method is that the ray tracing and the eikonal equation cannot overcome the multi-value travel time problem when the structure is complex, and therefore the accuracy of the dip angle gather is affected.

[0004] In the prior art, the angle of incidence of the incident ray from the shot point to the imaging point is obtained by ray tracing, and then the structural dip angle is obtained from the velocity model, and the angle of the scattered ray from the imaging point to the receiver point is obtained by scanning with the two angles. The problem faced by this method is that it is difficult to obtain an accurate structural dip angle from the velocity model in actual work, thereby affecting the calculation accuracy of the dip angle gather.

[0005] In the prior art, the Gaussian beam ray tracing method is used to calculate the angles of the incident ray and the scattered ray, which can adapt to more complex structures to some extent, but is still less adaptable than the wave equation method.

[0006] There is another method for calculating the dip angle gather in the prior art, which firstly obtains a high-frequency shot wave field propagation profile to the entire geological model by using a one-way wave equation simulation method, then obtains the incident ray angle information from the shot point to the imaging point from the profile by using a structure tensor technique, and then performs plane wave decomposition on the receiver wave field in the migration process, and calculates the scattering ray angle from the imaging point to the receiver by using the ray parameter corresponding to each plane wave component. The incident ray angle and the scattering ray angle obtained in the foregoing are used to calculate the dip angle, and the common imaging point dip angle gather after migration is generated. This method needs to perform one-way wave simulation before migration, and the accuracy of the angle calculated by using the structure tensor technique on the single-frequency profile is not high. In the specific implementation process, the plane wave decomposition is performed on the CPU by using the recursive linear Radon transform, and the calculation efficiency is relatively low.

[0007] Therefore, there is a need for a method for obtaining a dip angle gather and imaging a seismic profile with higher efficiency and more accurate angle, so as to more accurately analyze the underground structure in geological exploration and analysis. SUMMARY

[0008] The present application provides a method for imaging a seismic profile, which solves the technical problem of defects in the description of the seismic wave field of the ray-type seismic profile imaging method in the face of multi-value travel time and complex structure, can efficiently and accurately obtain the dip angle gather of the seismic wave field of the complex underground structure, and obtain high-resolution seismic profile imaging.

[0009] The present application provides a method for imaging a seismic profile, which includes:

[0010] Based on the common shot point seismic data, a one-way wave phase shift and interpolation algorithm are used to respectively obtain a first seismic wave field on a preset depth slice extending downward from the shot point and a second seismic wave field on a preset depth slice extending downward from the receiver point on one side of the receiver point;

[0011] A frequency Radon transform is performed on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component of a first preset frequency, a norm square operation is performed on each first plane wave component in a preset incident ray angle range to obtain the numerical value of each first plane wave component and sort it, and a first ray parameter is obtained, the first ray parameter is the ray parameter corresponding to the first plane wave component with the maximum numerical value, and a unique incident ray angle corresponding to the first ray parameter is obtained based on the first ray parameter and a ray parameter calculation formula;

[0012] Performing the frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point respectively to obtain each first plane wave component of a first preset frequency, performing Fourier transform on the first seismic wave field to obtain a third seismic wave field, and performing the frequency Radon transform on the third seismic wave field again to obtain each first plane wave component of the first preset frequency.

[0013] Based on the incident ray angle and each of the scattering ray angles, performing operation of half of a difference value between the incident ray angle and each of the scattering ray angles to obtain each of the dip angle;

[0014] Based on each of the dip angles, the incident ray angle, and each of the scattering ray angles, inputting into a single-shot imaging calculation formula to obtain a single-shot imaging function of each of the dip angles;

[0015] Performing the process of the single-shot imaging function on all shot points to obtain imaging functions of all shot points to obtain a dip angle gather, and stacking the dip angle gather to obtain a seismic profile image.

[0016] In the embodiments of the present application,

[0017] The step of performing the frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point respectively to obtain each first plane wave component of a first preset frequency comprises:

[0018] Performing Fourier transform on the first seismic wave field to obtain a third seismic wave field;

[0019] And performing the frequency Radon transform on the third seismic wave field again to obtain each first plane wave component of the first preset frequency.

[0020] In the embodiments of the present application,

[0021] The calculation formula of the Fourier transform is:

[0022]

[0023] Wherein,

[0024] f(x, t) is a time-space domain seismic wave field,

[0025] F(ω, x) is a frequency-space domain seismic wave field after Fourier transform,

[0026] ω is frequency, t is time, x is a horizontal coordinate of a single-shot imaging point, and i is an imaginary index;

[0027] The calculation formula of the frequency Radon transform is:

[0028] The calculation formula of the frequency Radon transform is:

[0029] wherein,

[0030] G(ω, p) is a plane wave component after Radon transform of frequency,

[0031] F(ω, x) is a frequency-space domain seismic wave field after Fourier transform,

[0032] ω is frequency, x is horizontal coordinate of a single shot imaging point, and i is imaginary index,

[0033] [-δx, +δx] is a range of a local spatial window centered at the imaging point x,

[0034] p k is a ray parameter representing a propagation direction of the kth plane wave component.

[0035] In an embodiment of the present application,

[0036] The step of obtaining a numerical value of each of the first plane wave components and sorting them by performing a square operation of a norm of each of the first plane wave components within a preset incident ray angle range to obtain a first ray parameter, the first ray parameter being the ray parameter corresponding to the first plane wave component with the maximum numerical value, comprises:

[0037] Setting an upper limit value and a lower limit value of the preset incident ray angle;

[0038] Based on the upper limit value, the lower limit value of the preset incident ray angle and the ray parameter calculation formula, obtaining the maximum value and the minimum value of the ray parameter corresponding to the upper limit and the lower limit of the preset incident ray angle;

[0039] Based on the maximum value and the minimum value of the ray parameter, performing average sampling, and obtaining a sampling interval of the ray parameter by dividing the difference between the maximum value and the minimum value of the ray parameter by the sampling number;

[0040] Based on the sampling interval and the minimum value of the ray parameter, obtaining each ray parameter corresponding to the propagation direction of each of the first plane wave components by sequentially multiplying the sampling interval by the product of the sampling number and adding the minimum value of the ray parameter;

[0041] Based on the calculation formula of the Radon transform, the third seismic wave field, the first preset frequency, and the horizontal coordinate of the single shot imaging point, obtaining each of the first plane wave components within the preset incident ray angle range;

[0042] Performing a square operation of a norm of each of the first plane wave components to obtain a numerical value of each of the first plane wave components and sorting them in descending order;

[0043] The ray parameter corresponding to the first plane wave component with the maximum value is the first ray parameter.

[0044] In the embodiments of the present application,

[0045] The ray parameter calculation formula is:

[0046]

[0047] wherein p is the ray parameter, θ is the angle, and v is the plane wave propagation velocity, which is generally obtained by measurement.

[0048] In the embodiments of the present application,

[0049] After the step of obtaining the imaging results of all the shot points in the process of cyclically performing the single-shot imaging results of all the shot points to obtain the dip angle gather, and stacking the dip angle gather to obtain the seismic profile imaging, the method further comprises:

[0050] Performing underground structure analysis based on the seismic imaging profile.

[0051] In the embodiments of the present application,

[0052] The single-shot imaging calculation formula is:

[0053]

[0054] wherein,

[0055] ω is the frequency,

[0056] x is the horizontal coordinate of the single-shot imaging point,

[0057] z is the depth coordinate of the single-shot imaging point,

[0058] θ S is the incident ray angle, θ R is the scattered ray angle, θ D is the dip angle,

[0059] I(x,z,θ D ) is the single-shot imaging function of the imaging point at a preset spatial position in the underground,

[0060] S(x,z,θ S ; ω) is the seismic wave field of all effective frequencies and corresponding unique incident ray angles,

[0061] R(x,z,θ R ; ω) is the seismic wave field of all effective frequencies and corresponding scattered ray angles.

[0062] The present application provides a device for seismic profile imaging, comprising:

[0063] An acquisition seismic wave field module is configured to obtain a first seismic wave field on a preset depth slice along a downward direction from a shot point and a second seismic wave field on a preset depth slice along a downward direction from a receiver point based on common shot seismic data by using a one-way wave phase shift and interpolation algorithm;

[0064] A calculation incident angle module is configured to perform a linear frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component of a first preset frequency, perform a square norm operation on each of the first plane wave components within a preset incident ray angle range to obtain a value of each of the first plane wave components and sort the values, and obtain a first ray parameter, which is a ray parameter corresponding to the first plane wave component with a maximum value, and obtain a unique incident ray angle corresponding to the first ray parameter based on the first ray parameter and a ray parameter calculation formula.

[0065] A calculation scattering angle module is configured to perform the frequency Radon transform on the second seismic wave field of all effective frequencies of the receiver point to obtain each second plane wave component corresponding to the all effective frequencies, obtain a value of each of the second ray parameters under the all effective frequencies by average sampling within a preset scattering ray angle range, and obtain each scattering ray angle corresponding to each of the second ray parameters based on each of the second ray parameters and the ray parameter calculation formula.

[0066] A calculation single-shot dip angle module is configured to perform an operation of half of a difference between the incident ray angle and each of the scattering ray angles based on the incident ray angle and each of the scattering ray angles to obtain each dip angle.

[0067] A generation single-shot imaging module is configured to obtain a single-shot imaging function of each of the dip angles by inputting each of the dip angles, the incident ray angle and each of the scattering ray angles into a single-shot imaging calculation formula.

[0068] A generation seismic profile imaging module is configured to obtain an imaging function of all shot points by performing a process of the single-shot imaging function on all shot points to obtain a dip angle gather, and obtain a seismic profile image by stacking the dip angle gather.

[0069] The application provides a storage medium, which has a computer program stored thereon,

[0070] The program is executed by a processor to implement steps of a method of seismic profile imaging according to any one of the above.

[0071] The application provides an electronic device, which comprises:

[0072] a memory having a computer program stored thereon; and

[0073] a processor for executing the computer program in the memory to implement the steps of the method of seismic profile imaging of any one of the above.

[0074] Compared with the prior art, one or more embodiments of the present application can have the following advantages:

[0075] The present application obtains the corresponding plane wave component by performing Radon transform on the linear frequency of the seismic wave field of the seismic source and the receiver respectively, obtains the unique angle of the incident ray corresponding to the maximum energy plane wave component in each plane wave component of the seismic source, then obtains multiple scattering angles of the receiver, thereby obtaining multiple dip angle angles, applies the imaging condition to the multiple dip angle angles to obtain the single shot imaging function, and stacks the dip angle gathers of all shot points to obtain the seismic profile imaging. The method of the present application can efficiently and accurately obtain the dip angle gather of the seismic wave field of the complex underground structure and obtain high-resolution seismic profile imaging.

[0076] Other features and advantages of the present application will be set forth in the following description, and in part will become apparent to those skilled in the art from the description, or can be learned by practice of the present application. The objects and other advantages of the present application can be realized and achieved by the structure particularly pointed out in the specification, claims and drawings. BRIEF DESCRIPTION OF DRAWINGS

[0077] The accompanying drawings are included to provide a further understanding of the present application, and constitute a part of the specification, illustrate embodiments of the present application, and are used together with the present application to explain the present application, and do not constitute a limitation of the present application. In the drawings:

[0078] Figure 1 is a flowchart of a method of seismic profile imaging of embodiment 1 of the present application;

[0079] Figure 2 is a schematic diagram of the angle of incidence, scattering angle and dip angle of embodiment 1 of the present application;

[0080] Figure 3 is a flowchart of a method of generating a dip angle gather of embodiment 1 of the present application;

[0081] Figure 4 is a schematic diagram of a Sigbee2A seismic model profile of embodiment 1 of the present application;

[0082] Figure 5 is a schematic diagram of a depth domain common imaging point dip angle gather of embodiment 1 of the present application;

[0083] Figure 6 is a schematic diagram of a dip angle gather of embodiment 1 of the present application on a seismic model to form a migration profile;

[0084] Figure 7Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application. DETAILED DESCRIPTION

[0085] In order to make the objects, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the drawings, so that the realization process of how the present application applies technical means to solve technical problems and achieve technical effects can be fully understood and implemented. It should be noted that, as long as there is no conflict, each embodiment in the present application and each feature in each embodiment can be combined with each other, and the technical solutions formed thereby are all within the protection scope of the present application.

[0086] First embodiment

[0087] Figure 1 Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application.

[0088] Figure 2 Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application.

[0089] Figure 3 Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application.

[0090] Figure 4 Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application.

[0091] Figure 5 Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application.

[0092] Figure 6 Fig. 1 is a schematic diagram of a logical framework of a device for seismic profile imaging according to Embodiment 2 of the present application.

[0093] In actual geological media, in addition to having layering, there are many special complex geological structures, such as faults, pinchouts, and the like, which constitute discontinuities (two-dimensional space) or discontinuous lines (three-dimensional space) of strata. When seismic waves propagate to these strata discontinuities (lines), they will be diffracted like light rays passing through a small hole in physical optics, and these discontinuities can be regarded as a new seismic source, from which a new disturbance is generated and propagates to the surrounding elastic space. This disturbance is called a diffracted wave in seismic exploration, and this phenomenon is called diffraction.

[0094] The present embodiment needs to use a graphics processing unit (English: Graphics Processing Unit, abbreviation: GPU), also known as a display core, a visual processor, or a display chip, which is a microprocessor specially used for image and graphics related operations on personal computers, workstations, game consoles, and some mobile devices (such as tablet computers, smart phones, etc.).

[0095] The embodiment provides a method for seismic profile imaging, comprising:

[0096] S100, obtaining a first seismic wave field on a preset depth slice along a downward direction from a shot point and a second seismic wave field along a downward direction from a receiver point based on common shot seismic data by using a one-way wave phase shift and interpolation algorithm.

[0097] The one-way wave migration algorithm is an important branch of seismic migration imaging, and is a method for decomposing a sound wave equation into an upward wave equation and a downward wave equation for wave field continuation calculation.

[0098] Specifically, in the embodiment, the one-way wave phase shift and interpolation algorithm is used to obtain the first seismic wave field on the depth slice along the downward direction from the shot point and the second seismic wave field along the downward direction from the receiver point based on the common shot seismic data. In this process, the one-way wave operator overcomes some shortcomings of the ray-based migration imaging method, such as defects in describing the seismic wave field of complex structures in the face of multi-value travel time. The one-way wave migration technology can better approximate the seismic wave field of complex structures, can finely image the complex structures, and can also ensure a high imaging frequency.

[0099] S110, performing frequency Radon transformation on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component of a first preset frequency, performing norm square operation on each first plane wave component in a preset incident ray angle range to obtain a value of each first plane wave component and sort the value, and obtaining a first ray parameter, and obtaining a unique incident ray angle corresponding to the first ray parameter based on the first ray parameter and a ray parameter calculation formula.

[0100] Specifically, in the embodiment, on the preset depth, the shot wavelet main frequency corresponding to the horizontal wave field component of the shot point source first seismic wave field is obtained, and linear frequency Radon transformation (i.e., Radon transformation) is performed on each horizontal position of the preset depth to obtain each first plane wave component of a first preset frequency. In the process of plane wave decomposition, a GPU acceleration technology is used to improve the calculation efficiency. Each first plane wave component is a function containing a preset frequency and a ray parameter, the ray parameter represents the propagation direction of each plane wave component, after obtaining all first plane wave components at the preset depth point, norm square operation is performed on each first plane wave component in a preset incident ray angle range to obtain a value of each first plane wave component and sort the value, and a first ray parameter is obtained. The first ray parameter is the ray parameter corresponding to the first plane wave component with the maximum value. The incident angle of the imaging point at the preset depth underground is calculated as the unique angle of the incident ray of this shot point at the first preset frequency by using the first ray parameter corresponding to the maximum energy plane wave component and a ray parameter calculation formula.

[0101] In the embodiment, the step of performing the frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point respectively to obtain each first plane wave component of the first preset frequency comprises:

[0102] The first seismic wave field in the time-space domain is converted by Fourier transform to obtain a third seismic wave field in the frequency-space domain;

[0103] The third seismic wave field is converted by the frequency Radon transform to obtain each first plane wave component of the first preset frequency.

[0104] In the embodiment, the calculation formula of the Fourier transform is:

[0105]

[0106] wherein,

[0107] f(x, t) is the first seismic wave field in the time-space domain,

[0108] F(ω, x) is the third seismic wave field in the frequency-space domain after the Fourier transform,

[0109] ω is the frequency, t is the time, x is the horizontal coordinate of the single-shot imaging point, and i is the imaginary index;

[0110] In the embodiment, the calculation formula of the ray parameter-frequency wave field conversion is:

[0111]

[0112] wherein,

[0113] G(ω, p) is the first plane wave component after the frequency Radon transform,

[0114] F(ω, x) is the third seismic wave field in the frequency-space domain after the Fourier transform,

[0115] ω is the frequency, t is the time, x is the horizontal coordinate of the single-shot imaging point, and i is the imaginary index,

[0116] p k is the ray parameter, representing the propagation direction of the kth plane wave component.

[0117] In the embodiment, the calculation formula of the ray parameter is:

[0118]

[0119] wherein, p is the ray parameter, θ is the angle, and v is the plane wave propagation speed, which is generally obtained by measurement.

[0120] Specifically, in the embodiment, the numerical values of the first plane wave components are obtained by performing the square operation of the norm of each first plane wave component in the preset incident ray angle range, and the first ray parameter is obtained by sorting the numerical values, and the step of obtaining the first ray parameter corresponding to the first plane wave component with the maximum numerical value includes:

[0121] 1) setting the upper limit value and the lower limit value of the preset incident ray angle;

[0122] 2) obtaining the maximum value and the minimum value of the ray parameter corresponding to the upper limit value and the lower limit value of the preset incident ray angle based on the upper limit value, the lower limit value and the ray parameter calculation formula of the preset incident ray angle;

[0123] 3) performing average sampling based on the maximum value and the minimum value of the ray parameter, and obtaining the sampling interval of the ray parameter by dividing the difference between the maximum value and the minimum value of the ray parameter by the sampling number;

[0124] 4) obtaining each ray parameter corresponding to the propagation direction of each first plane wave component based on the sampling interval and the minimum value of the ray parameter by multiplying the sampling interval by the product of the sampling number and adding the minimum value of the ray parameter;

[0125] 5) obtaining each first plane wave component in the preset incident ray angle range based on the ray parameter-frequency wave field conversion calculation formula, the third seismic wave field, the first preset frequency and the horizontal coordinate of the single shot imaging point;

[0126] 6) obtaining the numerical values of the first plane wave components by performing the square operation of the norm of each first plane wave component, and sorting the numerical values in descending order;

[0127] 7) selecting the ray parameter corresponding to the first plane wave component with the maximum numerical value as the first ray parameter.

[0128] Specifically, in the embodiment, the upper limit value and the lower limit value of the maximum possible angle of the incident ray angle are first set, and the preset incident ray angle is adjusted according to the actual exploration situation, and the embodiment is preferably [-30 degrees, +30 degrees]. Since the ray parameter p is related to the angle θ by the calculation formula, and the propagation velocity v of the plane wave can be obtained by measurement, the maximum value and the minimum value of the value of the ray parameter p corresponding to the upper limit value and the lower limit value of the preset incident ray angle can be calculated based on the upper limit value and the lower limit value of the preset incident ray angle. The maximum value pmax and the minimum value pmin of the numerical value of the ray parameter p are based on average sampling, and the sampling interval of the ray parameter is obtained by dividing the difference between the maximum value and the minimum value of the ray parameter by the sampling number.

[0129] The sampling interval is calculated as follows: Δp = (p max - pmin) / n, where Δp represents the sampling interval of the ray parameters.

[0130] Based on the sampling interval Δp and the minimum value pmin of the ray parameter, the ray parameters corresponding to the propagation direction of each first plane wave component are obtained by successively multiplying the sampling interval by the product of the number of samplings and adding the minimum value of the ray parameter. Then, the k-th ray parameter p k It can be represented as,

[0131] p k =p min + Δp*k

[0132] From this, the n ray parameters p can also be calculated. k The value of p is n, for example, p∈[p1, p2, ..., p...]. n ].

[0133] Based on the first preset frequency ω0 and related parameter data, the transformed first plane wave is obtained as G(ω0, p). k Based on the calculation formula of frequency Radon transform, the third seismic wave field, the first preset frequency, and the horizontal coordinate x of the single-shot imaging point, the first plane wave components G(ω0, p) within the preset incident ray angle range are obtained. k ).

[0134] Specifically, other parameters and each ray parameter p k Substituting into the formula for calculating the frequency Radon transform, for each first plane wave component G(ω0, p) k Perform norm squaring operation The values ​​of each first plane wave component are obtained and sorted in descending order. The ray parameter corresponding to the first plane wave component with the largest value is selected as the first ray parameter.

[0135] Substituting the first ray parameter and the plane wave propagation velocity v into the ray parameter calculation formula, we obtain the angle formula θ = arcsin(pv), thus obtaining the incident ray angle θ corresponding to the first ray parameter. S Different ray parameters p represent different plane wave components, representing wave fields propagating in different directions. Therefore, the incident angle θ at the underground imaging point is calculated from the maximum first ray parameter corresponding to the first plane wave component with the highest energy. S This determined a unique angle of incidence, such as Figure 2 As shown.

[0136] In this embodiment, since only one incident ray angle is needed on the source side of the shot point, the amount of calculation is greatly reduced, the calculation speed is accelerated, and the calculation cost is reduced.

[0137] S120, respectively performing frequency Radon transform on the second seismic wave field of all effective frequencies at the detection point to obtain each second plane wave component corresponding to all effective frequencies, obtaining the value of each second ray parameter under all effective frequencies by averaging sampling in the preset scattering ray angle range, and obtaining each scattering ray angle corresponding to each second ray parameter based on each second ray parameter and the ray parameter calculation formula.

[0138] Specifically, in the embodiment, by performing linear frequency Radon transform (i.e., Radon transform) on the second seismic wave field of all effective frequencies at the detection point side of the preset depth slice on the GPU device, the method similar to S110 obtains each second plane wave component related to the detection point, each second plane wave component represents a wave field of a specific propagation angle, and therefore each second plane wave component corresponds to a scattering angle of a scattered wave from the imaging point to the ground surface. The second plane wave component is a function containing frequency and second ray parameter. The use of the GPU device for calculation improves the calculation efficiency.

[0139] The embodiment does not perform screening of the second ray parameter corresponding to the maximum plane wave for the detection point side. For all second ray parameters p of the detection point side, such as n ray parameters, p∈[p1, p2, …, pn], the corresponding scattering angle θ is calculated. n R .

[0140] The step of obtaining the second ray parameter p is:

[0141] 1) Set the upper limit value and the lower limit value of the preset incident ray angle, for example, set to [-30 degrees, +30 degrees];

[0142] 2) Based on the upper limit value, the lower limit value of the preset incident ray angle, and the ray parameter calculation formula, obtain the maximum value and the minimum value of the ray parameter corresponding to the upper limit and the lower limit of the preset incident ray angle;

[0143] 3) Average sampling based on the maximum value and the minimum value in the value of each ray parameter, and obtain the sampling interval of the ray parameter by dividing the difference between the maximum value and the minimum value of the ray parameter by the sampling number;

[0144] The sampling interval calculation formula is as follows: Δp=(p max-pmin) / n, Δp represents the sampling interval of the ray parameter.

[0145] 4) Based on the sampling interval and the minimum value of the ray parameter, obtain each ray parameter corresponding to the propagation direction of each first plane wave component by multiplying the sampling interval by the product of the sampling number and adding the minimum value of the ray parameter;

[0146] Then the kth ray parameter p k ​may be expressed as, p k = p min+ Δp*k.

[0147] After obtaining each second ray parameter, each scattering angle θ R corresponding to each second ray parameter is obtained according to the ray parameter calculation formula.

[0148] In this embodiment, since all the second plane wave components on the receiver side are reserved, multiple scattering wave fields are reserved, and the multi-value travel time effect of the complex seismic wave field is contained.

[0149] S130, based on the incident ray angle and each scattering ray angle, the half of the difference between the incident ray angle and each scattering ray angle is calculated to obtain each dip angle;

[0150] Specifically, in this embodiment, the unique incident ray angle θ S and the scattering ray angle θ R corresponding to any second plane wave component on the receiver side are obtained, and the dip angle calculation formula is used to calculate the difference between the incident ray angle and a certain scattering ray angle, and then divided by 2 to obtain a dip angle θ D . Repeating the calculation of all scattering ray angles can obtain multiple dip angles corresponding to the unique incident ray angle and multiple scattering ray angles, as shown in Figure 2 .

[0151] The dip angle calculation formula is:

[0152] S140, based on each dip angle, incident ray angle, and each scattering ray angle, the single shot imaging function of each dip angle is obtained by using the single shot imaging calculation formula.

[0153] Specifically, in this embodiment, the incident wave field on the shot side and a certain plane wave component on the receiver side are cross-correlated by the single shot imaging calculation formula to obtain a single shot imaging function, wherein the incident angles of all frequencies of the incident wave field are set to be uniform θ S , S(x, z, θ S ; ω) is the seismic wave field of all effective frequencies and the unique incident ray angle corresponding to each effective frequency, and R(x, z, θ R ; ω) is the seismic wave field of all effective frequencies and each scattering ray angle corresponding to each effective frequency. Selecting only the unique incident ray angle corresponding to each effective frequency on the shot side can greatly reduce the calculation amount, but does not reduce the imaging effect of the single shot imaging.

[0154] Applying the imaging condition to all plane wave components on the receiver side and the incident wave field can obtain the single shot imaging results on multiple dip angles.

[0155] In this embodiment, the single-shot imaging calculation formula is:

[0156]

[0157] in,

[0158] ω is the frequency.

[0159] x is the horizontal coordinate of the single-shot imaging point.

[0160] z represents the depth coordinates of the point of view formed by a single shot.

[0161] θ S Let θ be the angle of the incident ray. R Let θ be the angle of the scattered ray. D The tilt angle is...

[0162] I(x,z,θ D ) represents the single-shot imaging function for the pre-defined spatial location imaging points underground.

[0163] S(x,z,θ S ;ω) represents the seismic wave field with all effective frequencies and the corresponding unique incident ray angle.

[0164] R(x,z,θ R ;ω) represents the seismic wave field with all effective frequencies and the corresponding angles of each scattered ray.

[0165] S150 is a process of iteratively executing the single-shot imaging function on all shot points to obtain the imaging function of all shot points, thereby obtaining the dip gather. The dip gathers are then superimposed to obtain the seismic profile image.

[0166] By adopting the method of this embodiment, accurate dip gathers can be obtained efficiently, thereby further obtaining high-resolution, high-signal-to-noise ratio seismic profile imaging.

[0167] In this embodiment, after obtaining the imaging results of all shot points by cyclically performing single-shot imaging on all shot points to obtain the dip gather, and then superimposing the dip gather to obtain the seismic profile imaging, the following steps are also included:

[0168] Subsurface structural analysis based on seismic imaging profiles.

[0169] like Figure 4 As shown, the international standard Sigsbee2A model contains a large number of reflection horizons, faults, discontinuities, salt domes, etc., as indicated by the black arrows. These represent geological bodies of varying sizes, such as faults, salt dome boundaries, and isolated diffraction points. The corresponding shot gather seismic data generated by finite difference forward modeling also contains a large amount of reflected wave energy and diffracted wave energy.

[0170] This embodiment uses this famous model and corresponding model data to test the effect, forming the depth domain common imaging point dip angle gather, as shown in Figure 5 .

[0171] Figure 5 The spatial position of the dip angle gather in Figure 6 corresponds to the spatial position indicated by the white line in Figure 5 , and the corresponding position has two isolated diffraction points at different depths, as indicated by the two white arrows in Figure 5 . In , the positive and negative 70-degree dip angles are taken, and the arrangement of the corresponding energy of the dip angle gather from negative 70 degrees to positive 70 degrees can be seen that at the position of the reflection point, the specific reflection stratum at the specific imaging point only corresponds to a specific geological dip angle, so all the reflection energy is relatively concentrated around this dip angle; the isolated diffraction point has no specific geological dip angle, so the diffraction energy is distributed in a very wide range of geological dip angles from the gather.

[0172] Figure 6 is the migration profile obtained by stacking the dip angle gather, from which it can be seen that many faults, salt dome boundaries and isolated diffraction points are imaged very clearly.

[0173] In summary, the embodiment provides a method for seismic profile imaging, which efficiently obtains accurate dip angle gathers by calculating the plane wave component of the preset frequency with the maximum value on the side of the shot point source, obtaining the ray parameter corresponding to the maximum plane wave component, obtaining the unique incident ray angle, and retaining the multiple scattering angles on the side of the receiver, thereby obtaining high-resolution and high-signal-to-noise ratio seismic profile imaging, improving the accuracy of underground structure analysis and facilitating the development of geological exploration work.

[0174] Second embodiment

[0175] Figure 7 is a logical framework schematic diagram of a device for seismic profile imaging according to the embodiment.

[0176] The embodiment provides a device for seismic profile imaging, which comprises:

[0177] The seismic wave field acquisition module is configured to obtain a first seismic wave field on a depth slice of the shot point and a second seismic wave field on the side of the receiver by using one-way wave phase shift and interpolation algorithms based on common shot point seismic data.

[0178] The incidence angle calculation module is used for performing linear frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component of a first preset frequency, performing norm square operation on each first plane wave component in a preset incidence ray angle range to obtain a numerical value of each first plane wave component and sorting, and obtaining a first ray parameter, the first ray parameter being the ray parameter corresponding to the first plane wave component with the maximum numerical value, and obtaining a unique incidence ray angle corresponding to the first ray parameter based on the first ray parameter and a ray parameter calculation formula;

[0179] The scattering angle calculation module is used for performing the frequency Radon transform on the second seismic wave field of all effective frequencies of the receiver point to obtain each second plane wave component corresponding to the all effective frequencies, obtaining the numerical value of each second ray parameter under the all effective frequencies by average sampling in a preset scattering ray angle range, and obtaining each scattering ray angle corresponding to each second ray parameter based on each second ray parameter and the ray parameter calculation formula.

[0180] The single-shot dip angle calculation module is used for performing the operation of half of the difference between the incidence ray angle and each scattering ray angle based on the incidence ray angle and each scattering ray angle to obtain each dip angle.

[0181] The single-shot imaging generation module is used for obtaining the single-shot imaging function of each dip angle by bringing each dip angle, the incidence ray angle and each scattering ray angle into a single-shot imaging calculation formula.

[0182] The seismic profile imaging generation module is used for performing the process of the single-shot imaging result for all shot points to obtain the imaging result of all shot points and thus obtain a dip angle gather, and superimposing the dip angle gather to obtain seismic profile imaging.

[0183] The embodiment provides a device for seismic profile imaging, which obtains the plane wave component of a preset frequency with the maximum value on the shot point side, obtains the ray parameter corresponding to the plane wave component with the maximum value, obtains a unique incidence ray angle, retains multiple scattering angles on the receiver point side, efficiently obtains accurate dip angle gathers, and thus obtains high-resolution and high-signal-to-noise-ratio seismic profile imaging, improves the accuracy of underground structure analysis, and facilitates the development of geological exploration work.

[0184] Third embodiment

[0185] The application provides a storage medium having a computer program stored thereon,

[0186] The program is executed by a processor to implement the steps of the method for seismic profile imaging according to any one of the above embodiments.

[0187] Those skilled in the art will appreciate that embodiments of the application can be devised for a method, a system, or a computer program product. Accordingly, the present application can be embodied in the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects that can all generally be referred to herein as a "service". Furthermore, the present application can take the form of a computer program product on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROMs, optical storage devices, etc.) embodying a computer program of instructions executable by a computer.

[0188] Those skilled in the art will appreciate that all or portions of the methods described above can be embodied in the form of a computer program product stored in a non-transitory computer readable storage medium, which includes steps for execution by a computer. The steps can be embodied in the form of a computer program product stored in a non-transitory computer readable storage medium, which includes steps for execution by a computer. The computer program product can include a computer readable storage medium having stored thereon the computer program of instructions. Accordingly, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects that can all generally be referred to herein as a "service". Furthermore, the present application can take the form of a computer program product on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROMs, optical storage devices, etc.) embodying a computer program of instructions executable by a computer.

[0189] Fourth embodiment

[0190] The present application provides an electronic device, comprising:

[0191] a memory having stored thereon a computer program; and

[0192] a processor configured to execute the computer program stored in the memory to implement the steps of the method of seismic profile imaging according to any one of the above.

[0193] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks. Figure 1 The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks. Figure 1 The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks.

[0194] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks. Figure 1 The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks. Figure 1 The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks.

[0195] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks. Figure 1 The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks. Figure 1 The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart or flows and / or block diagram block or blocks.

[0196] Although the present disclosure has been described with reference to embodiments thereof, it is to be understood that the present disclosure is not limited to the embodiments disclosed but is susceptible to various modifications and alterations without departing from the scope and spirit of the present disclosure.

Claims

1. A method of seismic profile imaging, characterized by, The method comprises the following steps: obtaining a first seismic wave field on a preset depth slice of a downward continuation of a shot point and a second seismic wave field on a downward continuation of a side of a receiver point based on common shot point seismic data by using a one-way wave phase shift and interpolation algorithm; performing frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component of a first preset frequency, performing square norm operation on each first plane wave component within a preset incident ray angle range to obtain a value of each first plane wave component and sort the value, and obtaining a first ray parameter, which is a ray parameter corresponding to the first plane wave component with the maximum value, and obtaining a unique incident ray angle corresponding to the first ray parameter based on the first ray parameter and a ray parameter calculation formula; performing the frequency Radon transform on the second seismic wave field of all effective frequencies of the receiver point to obtain each second plane wave component corresponding to the all effective frequencies, and obtaining a value of each second ray parameter under the all effective frequencies by average sampling within a preset scattering ray angle range, and obtaining each scattering ray angle corresponding to each second ray parameter based on each second ray parameter and the ray parameter calculation formula; performing half value operation of a difference between the incident ray angle and each scattering ray angle based on the incident ray angle and each scattering ray angle to obtain each dip angle; obtaining a single shot imaging function of each dip angle by inputting each dip angle, the incident ray angle and each scattering ray angle into a single shot imaging calculation formula; performing the process of the single shot imaging function for all shot points to obtain imaging functions of all shot points, thereby obtaining a dip angle gather, and stacking the dip angle gather to obtain a seismic profile imaging.

2. The method of claim 1, wherein, The step of performing frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component of a first preset frequency comprises the following steps: obtaining a third seismic wave field by Fourier transform of the first seismic wave field; obtaining each first plane wave component of the first preset frequency by frequency Radon transform of the third seismic wave field.

3. The method according to claim 2, wherein the calculation formula of the Fourier transform is: wherein f(x, t) is a time-space domain seismic wave field, F(ω, x) is a frequency-space domain seismic wave field after Fourier transform, ω is frequency, t is time, x is horizontal coordinate of a single shot imaging point, and i is an imaginary index. The calculation formula of the frequency Radon transform is: wherein G(ω, p) is a plane wave component after frequency Radon transform, F(ω, x) is a frequency-space domain seismic wave field after Fourier transform, ω is frequency, x is horizontal coordinate of a single shot imaging point, i is an imaginary index, and pk is a ray parameter representing a propagation direction of the kth plane wave component. ​ ​ ​ ​ ​ ​ ​ ​ ​ the range of the local spatial window centered at the imaging point x, ​ 4. The method of claim 2, wherein, The step of obtaining the numerical value of each first plane wave component by performing the square operation of the norm of each first plane wave component in the preset incident ray angle range and sorting the numerical value of each first plane wave component in descending order to obtain the first ray parameter, wherein the first ray parameter is the ray parameter corresponding to the first plane wave component with the maximum numerical value, comprises the following steps: setting the upper limit value and the lower limit value of the preset incident ray angle; obtaining the maximum value and the minimum value of the ray parameter corresponding to the upper limit value and the lower limit value of the preset incident ray angle based on the upper limit value, the lower limit value and the ray parameter calculation formula of the preset incident ray angle; performing average sampling based on the maximum value and the minimum value of the ray parameter, and obtaining the sampling interval of the ray parameter by dividing the difference between the maximum value and the minimum value of the ray parameter by the sampling number; obtaining the ray parameter corresponding to the propagation direction of each first plane wave component based on the sampling interval and the minimum value of the ray parameter by multiplying the sampling interval by the product of the sampling number and adding the minimum value of the ray parameter; obtaining each first plane wave component in the preset incident ray angle range based on the frequency Radon transform calculation formula, the third seismic wave field, the first preset frequency and the horizontal coordinate of the single-shot imaging point; performing the square operation of the norm of each first plane wave component to obtain the numerical value of each first plane wave component and sorting the numerical value of each first plane wave component in descending order; selecting the ray parameter corresponding to the first plane wave component with the maximum numerical value as the first ray parameter.

5. The method of claim 4, wherein the ray parameter calculation formula is: where p is the ray parameter, is the angle, and v is the plane wave propagation velocity, which is generally measured.

6. The method of claim 1, wherein, after the step of obtaining the imaging results of all shot points in the process of performing the single-shot imaging result on all shot points to obtain the dip angle gather, and stacking the dip angle gather to obtain the seismic profile imaging, the method further comprises: performing subsurface structure analysis based on the seismic imaging profile.

7. The method of claim 1, wherein the single-shot imaging calculation formula is: wherein, ω is the frequency, x is the horizontal coordinate of the single-shot imaging point, z is the depth coordinate of the single-shot imaging point, for the angle of scattered radiation, for the angle of inclination a single shot imaging function for imaging a point in a predetermined spatial location underground, for all valid frequencies and corresponding unique incident ray angles, is the seismic wavefield for all valid frequencies and corresponding scattering ray angles.

8. An apparatus for seismic profile imaging, characterized by the method further comprises: an acquisition of seismic wave field module, configured to obtain a first seismic wave field on a preset depth slice along the downward direction of a shot point and a second seismic wave field on a preset depth slice along the downward direction of a receiver based on common receiver seismic data by using one-way wave phase shift and interpolation algorithms; a calculation of incident angle module, configured to perform linear frequency Radon transform on the first seismic wave field of all effective frequencies of the shot point to obtain each first plane wave component at a first preset frequency, perform the square operation of the norm of each first plane wave component in the preset incident ray angle range to obtain the numerical value of each first plane wave component and sort the numerical value of each first plane wave component in descending order, and obtain a first ray parameter, wherein the first ray parameter is the ray parameter corresponding to the first plane wave component with the maximum numerical value, and obtain the unique incident ray angle corresponding to the first ray parameter based on the first ray parameter and a ray parameter calculation formula. a computing scattering angle module, configured to perform the frequency Radon transform on the second seismic wave field of all effective frequencies of the receiver point respectively to obtain each second plane wave component corresponding to the all effective frequencies, to obtain the value of each second ray parameter under the all effective frequencies by averaging sampling in a preset scattering ray angle range, and to obtain each scattering ray angle corresponding to each second ray parameter based on each second ray parameter and the ray parameter calculation formula; a computing single-shot dip angle module, configured to perform the operation of half of the difference between the incident ray angle and each scattering ray angle based on the incident ray angle and each scattering ray angle to obtain each dip angle; a generating single-shot imaging module, configured to obtain the single-shot imaging function of each dip angle by bringing each dip angle, the incident ray angle and each scattering ray angle into a single-shot imaging calculation formula; a generating seismic profile imaging module, configured to perform the process of the single-shot imaging function on all shot points, to obtain the imaging function of all shot points, and to obtain the dip angle gather by stacking the dip angle gather, and to obtain the seismic profile imaging.

9. A storage medium having a computer program stored thereon, and the computer program is executed by a processor to implement the steps of the method for seismic profile imaging according to any one of claims 1 to 7. including:

10. An electronic device, comprising: a memory having a computer program stored thereon; and a processor configured to execute the computer program in the memory to implement the steps of the method for seismic profile imaging according to any one of claims 1 to 7. ​

Images