A pre-stack depth migration imaging method based on generalized diffraction stack

Abstract

The present application belongs to the technical field of oil and gas seismic exploration, and particularly relates to a pre-stack depth migration imaging method based on generalized diffraction stacking. The present application uses one-way wave equation to calculate the Green's function required for imaging, thereby avoiding the defects of two-way wave equation, making the imaging result not containing the wave field component useless for imaging, ensuring good imaging effect, and reducing the calculation amount, accelerating the operation speed, and improving the imaging efficiency. Moreover, the Green's function is solved in the frequency domain, and there is no sampling rate dispersion problem in the time domain, so that the calculation mode of one-way wave equation can adapt to larger spatial grid interval. Moreover, the relative positions of the imaging point, the source point and the receiver point are judged through angle calculation, so as to perform regional decomposition on the velocity model, and only the necessary region needs to be recursively calculated when solving the Green's function, without recursively calculating in the whole velocity model space. On the basis of ensuring the recursive effect, the useless calculation region is removed, and the imaging efficiency is further improved.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas seismic exploration technology, specifically relating to a pre-stack depth migration imaging method based on generalized diffraction stacking. Background Technology

[0002] In recent years, with the increasing demands for precision in oil exploration, academia and industry have proposed a series of imaging methods, leading to significant advancements in pre-stack depth migration technology. Generally speaking, pre-stack depth migration imaging methods primarily rely on two types of computational methods: geometric ray methods and wave equation methods. Among these, the Generalized Diffraction Stack (GDS) method is a newly emerging migration imaging method that combines these two types. Essentially, it utilizes the basic idea of ​​geometric ray algorithms (Kirchhoff's method) while employing the solution method of wave equation methods during the migration process to solve for the numerical Green's function. In other words, it uses wave equation methods to achieve geometric ray migration. Therefore, this method simultaneously possesses the limited bandwidth advantage of wave equation methods and the target-oriented imaging advantage of geometric ray methods.

[0003] For example, the first issue of *Seg Technical Program Expanded Abstracts* in 2016 published Schuster's "Reverse-time migration = generalized diffraction stack migration," which proposed a generalized diffraction stacking method based on the two-way wave equation. It proved that this method, similar to reverse-time migration (RTM), has the advantage of imaging complex structures with steep dip angles, but also suffers from the disadvantages of RTM, such as high computational cost and excessive shallow noise. The *Chinese Journal of Engineering Geophysics*, Volume 14, Issue 6, 2017 published Yan Hongqun and Sun Jianguo's "A Depth Domain Migration Imaging Method Based on Generalized Diffraction Stacking," which implemented generalized diffraction stacking migration in China using the finite difference method of the two-way wave equation. Yan Hongqun's "Implementation Technology of Generalized Diffraction Stacking Migration," published in a 2018 master's thesis from Jilin University, optimized generalized diffraction stacking migration by adding a two-way illumination compensation term for the source and receiver points, building upon his previous research. In 2021, the doctoral dissertation of Jilin University published Meng Xiangyu's "Reflection Seismic Response and Related Processing Methods under Complex Marine Acoustic Environments". This work, targeting the characteristics of marine acoustic environments, derives the time-domain imaging formula of Kirchhoff migration into a reverse-time migration imaging formula, which is essentially a reverse-time migration. Through the implementation method of reverse-time migration, Green's function is used to solve for forward and backward propagation wavefield snapshots, and the imaging process is carried out in the time domain, that is, the form of the cross-correlation of the forward propagation term and the backward propagation term of the data along time and the integral.

[0004] The aforementioned generalized diffraction stacking method mainly utilizes the two-way wave equation to solve Green's function in the time domain, that is, to transform the frequency domain GDS imaging formula into the time domain imaging formula. However, the existing generalized diffraction stacking method using the two-way wave equation has a large computational load when performing pre-stack depth migration imaging, resulting in low computational efficiency; and the two-way wave equation contains wave field components such as reflected waves that are not beneficial to imaging, resulting in poor imaging effect. Summary of the Invention

[0005] The purpose of this invention is to provide a pre-stack depth migration imaging method based on generalized diffraction superposition, which solves the problems of low computational efficiency and poor imaging effect of the generalized diffraction superposition method using two-way wave equation in the prior art when performing pre-stack depth migration imaging.

[0006] To achieve the above objectives, the present invention provides a pre-stack depth migration imaging method based on generalized diffraction superposition, which reads relevant data of the target area required for imaging, sets grid points as imaging points at a set interval in the velocity model space, and performs cyclic imaging on each imaging point to obtain the pre-stack depth migration imaging result.

[0007] The relevant data for the target area required for imaging includes the velocity model corresponding to the offset velocity field, the location of the detector point, and the location of the seismic source point.

[0008] The step of performing cyclic imaging of each imaging point to obtain pre-stack depth migration imaging results includes the following steps:

[0009] 1) Select a grid point as the current imaging point, and select the recursive region and recursive direction of the one-way wave field from the source point to the imaging point, and the recursive region and recursive direction of the one-way wave field from the receiver point to the imaging point.

[0010] 2) Based on the recursive region and recursive direction of the single-pass wave field from the source point to the imaging point, and the recursive region and recursive direction of the single-pass wave field from the receiver point to the imaging point, Green's function from the source point to the imaging point and from the receiver point to the imaging point are solved recursively using the single-pass wave.

[0011] 3) Substitute the reflection seismic data of the target area, the solution of Green's function from the source point to the imaging point, and the solution of Green's function from the receiver point to the imaging point into the generalized diffraction superposition formula, and integrate along the receiver points and frequencies at different locations to obtain the imaging value of the current imaging point.

[0012] 4) Determine whether imaging of all grid points in the velocity model has been completed; if yes, proceed to step 5); if no, return to step 1) and continue imaging the next grid point.

[0013] 5) Output the imaging values ​​of all grid points in the velocity model as the final pre-stack depth migration imaging result.

[0014] The beneficial effects of the above technical solution are as follows: Green's function required for imaging is calculated using the one-way wave equation, rather than the two-way wave equation used in the traditional generalized diffraction superposition; thus avoiding the inherent defects of the two-way wave equation, ensuring that the imaging results do not contain wave field components that are not beneficial to imaging, guaranteeing good imaging results, and reducing the computational load of Green's function, speeding up the operation, and improving imaging efficiency.

[0015] Furthermore, the method for selecting the recursive region and recursive direction of the one-way wavefield from the source point to the imaging point is as follows:

[0016] In the velocity model, the positional relationship between the current imaging point and the seismic source point is determined. The formula for the positional relationship between the current imaging point and the seismic source point is:

[0017]

[0018] In the formula, A s This indicates the current imaging point position (x, y) relative to the source point position (x, y). s ,y s The angle of As; if the value of As is within the range of the first set source point angle, then it is determined that the imaging point is mainly located to the right of the source point. The velocity model is laterally decomposed into left and right sides along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region where the source point +x direction is located. S ;

[0019] If the value of As falls within the range of the second or third preset source point angle, the imaging point is determined to be mainly located to the left of the source point. The velocity model is then horizontally decomposed into left and right sides along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region in the -x direction of the source point. S ;

[0020] If the value of As falls within the range of the fourth set source point angle, then the imaging point is determined to be mainly located below the source point. The velocity model is longitudinally decomposed into upper and lower parts along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region containing the source point in the +y direction. S ;

[0021] If the value of As falls within the range of the fifth set source point angle, then the imaging point is determined to be mainly located above the source point. The velocity model is longitudinally decomposed into upper and lower parts along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region in the -y direction of the source point. S .

[0022] The beneficial effects of the above technical solution are as follows: by calculating the angle to determine the relative position between the imaging point and the source point, the recursive region and recursive direction of the one-way wave field from the source point to the imaging point are obtained by performing region decomposition. That is, the necessary recursive regions are selected according to the relative position between the imaging point and the source point, so that when recursively solving the Green's function from the source point to the imaging point, it is only necessary to perform the recursion in the necessary regions, without recursing in the full velocity model space. On the basis of ensuring the recursion effect, the useless calculation regions can be removed to a certain extent after the angle region decomposition, and the calculation efficiency of the Green's function from the source point to the imaging point can be further improved, thereby improving the overall imaging efficiency.

[0023] Furthermore, the method for selecting the recursive region and recursive direction of the one-way wavefield from the detector point to the imaging point is as follows:

[0024] In the velocity model, the positional relationship between the current imaging point and the detector point is determined. The formula for the positional relationship between the current imaging point and the detector point is:

[0025]

[0026] In the formula, A r This indicates the current imaging point position (x, y) relative to the detector point position (x, y). r ,y r Angle of A; if A r If the value is within the range of the first set detector angle, then it is determined that the imaging point is mainly located to the right of the detector. The velocity model is decomposed laterally into left and right sides along the detector. The Green's function Gr from the detector to the imaging point is solved only along the region where the detector +x direction is located.

[0027] If A r If the value is within the range of the second set detector angle or the third set detector angle, then the imaging point is determined to be mainly located to the left of the detector; the velocity model is decomposed laterally into left and right sides along the detector, and Green's function Gr from the detector to the imaging point is solved recursively only along the region where the detector is located in the -x direction;

[0028] If A r If the value is within the range of the fourth set detector angle, then it is determined that the imaging point is mainly located below the detector. The velocity model is longitudinally decomposed into upper and lower parts along the detector. The Green's function Gr from the detector to the imaging point is solved only along the region where the detector +y direction is located.

[0029] If A r If the value is within the range of the fifth set detector angle, then it is determined that the imaging point is mainly located above the detector. The velocity model is longitudinally decomposed into upper and lower parts along the detector, and Green's function Gr from the detector to the imaging point is solved only along the region where the detector is located in the -y direction.

[0030] The beneficial effects of the above technical solution are as follows: by calculating the angle to determine the relative position of the imaging point and the detector point, the recursive region and recursive direction of the one-way wavefield from the detector point to the imaging point are obtained by performing region decomposition. That is, the necessary recursive regions are selected according to the relative position of the imaging point and the detector point, so that when recursively solving the Green's function from the detector point to the imaging point, it is only necessary to perform the recursion in the necessary regions, without recursing in the full velocity model space. On the basis of ensuring the recursion effect, the angle region decomposition can remove useless calculation regions to a certain extent and further improve the calculation efficiency of the Green's function from the detector point to the imaging point, thereby improving the overall imaging efficiency.

[0031] Furthermore, the relevant data for the target region required for imaging also includes the source function in the frequency domain of the target region, seismic data, and seismic wavelet; the generalized diffraction superposition formula is:

[0032]

[0033] Where I(x) represents the imaging result, Gs and Gr are Green's functions from the source point and receiver point to the imaging point, respectively, d() represents the seismic data of the target area, xs and xr are the locations of the source point and receiver point, x is the grid position, and ω is the frequency. Indicates complex conjugation;

[0034] In the generalized diffraction superposition formula:

[0035]

[0036] Where s(ω) is the source function in the frequency domain of the target region, and f1(ω) and f2(ω) are the first and second decomposition terms obtained by the seismic wavelet decomposition corresponding to the source, respectively.

[0037] Furthermore, based on the recursive region and direction of the one-way wavefield from the source point to the imaging point, and the recursive region and direction of the one-way wavefield from the receiver point to the imaging point, the Green's function from the source point to the imaging point and from the receiver point to the imaging point is solved recursively using the one-way wavefield, respectively:

[0038] Using the first decomposition term f1(ω) of the seismic wavelet corresponding to the earthquake source as the boundary condition, the formula for recursively solving the Green's function Gs from the source point to the imaging point along the x-axis is as follows:

[0039]

[0040] In this formula, a "+" indicates recursively solving Green's function along the positive x-axis, and a "-" indicates recursively solving Green's function along the negative x-axis; i 2=-1, Gs is Green's function from the source point to the imaging point, and Ex is the one-way recursive operator for the wave field along the x-direction;

[0041] The formula for recursively solving Green's function Gs from the seismic source point to the imaging point along the y-axis is as follows:

[0042]

[0043] In this formula, a + sign indicates that the Green's function is solved recursively along the positive y-axis, and a - sign indicates that the Green's function is solved recursively along the negative y-axis; i 2 =-1, Gs is Green's function from the source point to the imaging point, and Ey is the one-way recursive operator for the wave field along the y direction;

[0044] Using the second decomposition term f2(ω) of the seismic wavelet corresponding to the earthquake source as the boundary condition, the formula for recursively solving Green's function Gr from the receiver point to the imaging point along the x-axis is as follows:

[0045]

[0046] In this formula, a "+" indicates recursively solving Green's function along the positive x-axis, and a "-" indicates recursively solving Green's function along the negative x-axis; i 2 =-1, Gr is Green's function from the detector point to the imaging point, and Ex is the one-way recursive operator for the wave field along the x-direction;

[0047] The formula for recursively solving Green's function Gr from the receiver point to the imaging point along the y-axis is as follows:

[0048]

[0049] In this formula, a + sign indicates that the Green's function is solved recursively along the positive y-axis, and a - sign indicates that the Green's function is solved recursively along the negative y-axis; i 2 =-1, Gr is Green's function from the detector point to the imaging point, and Ey is the one-way recursive operator for the wave field along the y-direction;

[0050] The approximate one-way recursive operators Ex and Ey, which expand the wave field along the x and y directions, are expressed as follows:

[0051]

[0052]

[0053] In the above formula, c represents the background velocity; v(x,y) represents the velocity field value; x and y are the coordinates of the imaging point in the horizontal and depth directions, respectively; and ω is the frequency.

[0054] The beneficial effects of the above technical solution are as follows: Green's function is solved in the frequency domain without having to convert the generalized diffraction superposition formula used for imaging to the time domain. Therefore, there is no sampling rate dispersion problem inherent in the time domain, which allows the calculation method of Green's function based on the one-way wave equation to adapt to a larger spatial grid interval. Attached Figure Description

[0055] Figure 1 This is a flowchart of the pre-stack depth migration imaging method based on generalized diffraction superposition in an embodiment of the present invention.

[0056] Figure 2 This is a schematic diagram of the velocity model corresponding to the migration velocity field of the target region in an embodiment of the pre-stack depth migration imaging method based on generalized diffraction superposition of the present invention.

[0057] Figure 3 This is a schematic diagram illustrating the region decomposition of the velocity model based on the positional relationship of the imaging point relative to the source point and the receiver point in an embodiment of the pre-stack depth migration imaging method based on generalized diffraction superposition of the present invention.

[0058] Figure 4 This is a schematic diagram of Green's function recursively derived from the source point to the imaging point in an embodiment of the pre-stack depth migration imaging method based on generalized diffraction superposition of the present invention.

[0059] Figure 5 This is a schematic diagram of Green's function recursively derived from the detector point to the imaging point in an embodiment of the pre-stack depth migration imaging method based on generalized diffraction superposition of the present invention.

[0060] Figure 6 This is a schematic diagram of the final pre-stack depth migration imaging result obtained by using the pre-stack depth migration imaging method based on generalized diffraction superposition in this embodiment of the present invention.

[0061] Figure 7 This is a schematic diagram of the imaging result obtained by Laplacian filtering and denoising using conventional generalized diffraction superposition imaging in an embodiment of the pre-stack depth migration imaging method based on generalized diffraction superposition of the present invention.

[0062] Figure 8 This is a comparison chart showing the computational efficiency of Green's function in the pre-stack depth migration imaging method based on generalized diffraction stacking of the present invention, using the conventional generalized diffraction stacking imaging method, the generalized diffraction stacking imaging method based on single-path wavefield recursion of this embodiment, and the generalized diffraction stacking imaging method based on single-path wavefield recursion after region decomposition of this embodiment. Detailed Implementation

[0063] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0064] Example of a pre-stack depth migration imaging method based on generalized diffraction superposition

[0065] This embodiment presents a technical solution for a pre-stack depth migration imaging method based on generalized diffraction superposition, referring to... Figure 1 The imaging method first reads the relevant data of the target area required for imaging, then sets grid points as imaging points in the velocity model space at a set interval, and finally performs cyclic imaging on each imaging point to obtain the pre-stack depth migration imaging result. In this embodiment, grid points are the term used during calculation, and imaging points are the term used during the final imaging. That is, grid points and imaging points are actually the same.

[0066] The relevant data for the target area required for imaging include the velocity model corresponding to the migration velocity field, the location of the receiver point, and the location of the seismic source point. In this embodiment, since imaging needs to be performed in the software system, the relevant data for the target area required for imaging that needs to be read in before imaging includes the reflection seismic data of the common shot domain (i.e., surface observation data obtained from the same shot point), the migration velocity field v(x,y), the seismic wavelet, and the observation system. In addition, based on the above information, the wavefield recursion parameters in the program also need to be given, including the spatial sampling interval (dx = dy = 10m, x represents horizontal distance, y represents depth) and the number of sampling points (nx = 401, ny = 301, x represents horizontal distance, y represents depth), the time sampling interval of the data (dt = 1.0ms) and the number of sampling points (nt = 1500), and the frequency sampling interval (dω = 1Hz) and the number of sampling points (nω = 155). These constitute the known conditions of this invention. Figure 2 The figure shown is a schematic diagram of the velocity model corresponding to the offset velocity field obtained in this embodiment. It can be seen that the minimum velocity of the model is about 1200m / s, the maximum velocity is about 5000m, and the model size is 4000m×3000m.

[0067] In this embodiment, the imaging points cycle from left to right and from top to bottom in the model, specifically as follows: Figure 2 In the velocity model shown, the horizontal direction is cycled from 0 meters to 4000 meters at 10-meter intervals, and the depth direction is cycled from 0 meters to 3000 meters, until all imaging points are cycled through. The specific steps for performing cyclic imaging on each imaging point to obtain the pre-stack depth migration imaging results are as follows:

[0068] 1) Select a grid point as the current imaging point, and select the recursive region and recursive direction of the one-way wave field from the source point to the imaging point, and the recursive region and recursive direction of the one-way wave field from the receiver point to the imaging point.

[0069] In this embodiment, the method for selecting the recursive region and recursive direction of the one-way wavefield from the source point to the imaging point is as follows:

[0070] In the velocity model, the positional relationship between the current imaging point and the seismic source point is determined. The formula for the positional relationship between the current imaging point and the seismic source point is:

[0071]

[0072] In the formula, A s This indicates the current imaging point position (x, y) relative to the source point position (x, y). s ,y s The angle of As; if the value of As is within the range of the first set source point angle, then it is determined that the imaging point is mainly located to the right of the source point. The velocity model is laterally decomposed into left and right sides along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region where the source point +x direction is located. S In this embodiment, the value of As falling within the first predetermined angle range of the seismic source point is represented as 45°≥A. S >-45°;

[0073] If the value of As falls within the range of the second or third preset source point angle, the imaging point is determined to be mainly located to the left of the source point. The velocity model is then horizontally decomposed into left and right sides along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region in the -x direction of the source point. S In this embodiment, the value of As falling within the range of the second predetermined seismic source point angle is represented as 180 ≥ A. S If the angle is greater than 135°, and the value of As falls within the range of the third preset seismic source angle, it is expressed as -135 ≥ A. S >-180°;

[0074] If the value of As falls within the range of the fourth set source point angle, then the imaging point is determined to be mainly located below the source point. The velocity model is longitudinally decomposed into upper and lower parts along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region containing the source point in the +y direction. S In this embodiment, the value of As falling within the fourth predetermined source point angle range is represented as 135° ≥ A. S >45°;

[0075] If the value of As falls within the range of the fifth set source point angle, then the imaging point is determined to be mainly located above the source point. The velocity model is longitudinally decomposed into upper and lower parts along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region in the -y direction of the source point. S In this embodiment, the value of As within the fifth predetermined source point angle range is represented as -45° ≥ A. S >-135°.

[0076] Similarly, the method for selecting the recursive region and recursive direction of the one-way wavefield from the detector point to the imaging point is as follows:

[0077] In the velocity model, the positional relationship between the current imaging point and the detector point is determined. The formula for the positional relationship between the current imaging point and the detector point is:

[0078]

[0079] In the formula, A r This indicates the current imaging point position (x, y) relative to the detector point position (x, y). r ,y r Angle of A; if A r If the value is within the range of the first set detector angle, then the imaging point is determined to be mainly located to the right of the detector. The velocity model is laterally decomposed into left and right sides along the detector, and Green's function Gr from the detector to the imaging point is recursively solved only along the region where the detector +x direction is located; in this embodiment, A r The value within the first set detector angle range is expressed as 45°≥A. r >-45°;

[0080] If A r If the value is within the range of the second or third preset detector angle interval, then the imaging point is determined to be mainly located to the left of the detector point; the velocity model is laterally decomposed into left and right sides along the detector point, and Green's function Gr from the detector point to the imaging point is recursively solved only along the region where the detector point is located in the -x direction; in this embodiment, A r The value within the second set detector angle range is expressed as 180 ≥ A. r >135°, A r The value within the range of the third set detector angle is expressed as -135 ≥ A. r >-180°;

[0081] If A rIf the value is within the range of the fourth set detector angle, then the imaging point is determined to be mainly located below the detector. The velocity model is longitudinally decomposed into upper and lower parts along the detector, and Green's function Gr from the detector to the imaging point is recursively solved only along the region where the detector +y direction is located; In this embodiment, A r The value within the fourth set detector angle range is expressed as 135°≥A. r >45°;

[0082] If A r If the value is within the range of the fifth set detector angle, then the imaging point is determined to be mainly located above the detector. The velocity model is longitudinally decomposed into upper and lower parts along the detector, and Green's function Gr from the detector to the imaging point is recursively solved only along the region where the detector is located in the -y direction. In this embodiment, A r The value within the fifth set detector angle range is expressed as -45° ≥ A. r >-135°.

[0083] like Figure 3 The diagram shows the positional relationship between imaging point I (1000m, 1000m), source point S (1000m, 50m), and receiver point r (3000m, 50m). Through angle calculation in this step, it is known that the angle of imaging point I relative to source point S is approximately 90°, indicating that it is mainly located below the source point. Therefore, the velocity model needs to be decomposed along the depth direction (i.e., Green's function Gs is recursively applied only below the source point). Similarly, the angle of imaging point I relative to receiver point r is approximately 154°, indicating that it is mainly located to the left of the receiver. Therefore, the velocity model needs to be decomposed along the horizontal direction (i.e., Green's function Gr is recursively applied only to the left of the receiver). After angle calculation and region decomposition, the following results are obtained: Figure 3 The recursive region in the model is selected, which means that the necessary recursive regions are selected. This allows the recursion of Green's function to be performed only in the necessary regions, without recursing in the entire model space. While ensuring the recursion effect, the corresponding amount of computation is saved and the overall imaging efficiency is improved.

[0084] 2) Based on the recursive region and recursive direction of the single-pass wave field from the source point to the imaging point, and the recursive region and recursive direction of the single-pass wave field from the receiver point to the imaging point, Green's function from the source point to the imaging point and from the receiver point to the imaging point are solved recursively using the single-pass wave.

[0085] In this embodiment, the generalized diffraction superposition formula used to calculate the image value is as follows:

[0086]

[0087] Where I(x) represents the imaging result, Gs and Gr are Green's functions from the source point and receiver point to the imaging point, respectively, d() represents the seismic data of the target area, xs and xr are the locations of the source point and receiver point, x is the grid position, and ω is the frequency. Indicates complex conjugation;

[0088] Furthermore, in this generalized superposition formula for diffraction:

[0089]

[0090] Wherein, s(ω) is the source function in the frequency domain of the target region, and f1(ω) and f2(ω) are the first and second decomposition terms obtained by the seismic wavelet decomposition corresponding to the source, respectively; in this embodiment, s(ω) is the result of the Ricker wavelet with a dominant frequency of 30Hz after Fourier transform from the time domain to the frequency domain; therefore, the relevant data of the target region required for imaging in this embodiment also includes the source function in the frequency domain of the target region, seismic data, and seismic wavelet, rather than the relevant data in the time domain of the current generalized diffraction superposition;

[0091] In this embodiment, the generalized diffraction superposition formulas Gs(xs|x,ω)f1(ω) and Gr(xr|x,ω)f2(ω) used to calculate the imaging values, obtained by solving Green's functions from the source point to the imaging point and from the receiver point to the imaging point, utilize the recursive method of the one-way wave equation in the frequency domain, instead of the conventional two-way wave equation recursive method for generalized diffraction superposition. Since the one-way wave equation is an approximate equation that only contains the transmitted wave component, it can fully leverage the advantages of the one-way wave equation and overcome the shortcomings of the two-way wave equation, such as the inclusion of components in the entire wave field that are not beneficial to imaging and the large amount of computation. This improves the generalized diffraction superposition method, resulting in a clearer final imaging result and higher computational efficiency in the imaging process.

[0092] The specific method for recursively solving Green's functions from the source point to the imaging point and from the receiver point to the imaging point using single-way waves, based on the recursive region and direction of the single-way wavefield from the source point to the imaging point, and the recursive region and direction of the single-way wavefield from the receiver point to the imaging point, is as follows:

[0093] Using the first decomposition term f1(ω) of the seismic wavelet corresponding to the earthquake source as the boundary condition, the formula for recursively solving the Green's function Gs from the source point to the imaging point along the x-axis is as follows:

[0094]

[0095] In this formula, a "+" indicates recursively solving Green's function along the positive x-axis, and a "-" indicates recursively solving Green's function along the negative x-axis; i 2=-1, Gs is Green's function from the source point to the imaging point, and Ex is the one-way recursive operator for the wave field along the x-direction;

[0096] The formula for recursively solving Green's function Gs from the seismic source point to the imaging point along the y-axis is as follows:

[0097]

[0098] In this formula, a + sign indicates that the Green's function is solved recursively along the positive y-axis, and a - sign indicates that the Green's function is solved recursively along the negative y-axis; i 2 =-1, Gs is Green's function from the source point to the imaging point, and Ey is the one-way recursive operator for the wave field along the y direction;

[0099] Using the second decomposition term f2(ω) of the seismic wavelet corresponding to the earthquake source as the boundary condition, the formula for recursively solving Green's function Gr from the receiver point to the imaging point along the x-axis is as follows:

[0100]

[0101] In this formula, a "+" indicates recursively solving Green's function along the positive x-axis, and a "-" indicates recursively solving Green's function along the negative x-axis; i 2 =-1, Gr is Green's function from the detector point to the imaging point, and Ex is the one-way recursive operator for the wave field along the x-direction;

[0102] The formula for recursively solving Green's function Gr from the receiver point to the imaging point along the y-axis is as follows:

[0103]

[0104] In this formula, a + sign indicates that the Green's function is solved recursively along the positive y-axis, and a - sign indicates that the Green's function is solved recursively along the negative y-axis; i 2 =-1, Gr is Green's function from the detector point to the imaging point, and Ey is the one-way recursive operator for the wave field along the y-direction;

[0105] The approximate one-way recursive operators Ex and Ey, which expand the wave field along the x and y directions, are expressed as follows:

[0106]

[0107]

[0108] In the above equation, c represents the background velocity; v(x,y) represents the velocity field value; x and y are the coordinates of the imaging point in the horizontal and depth directions, respectively; and ω is the frequency. By performing any Taylor expansion on the operators Ex and Ey in the above equation, Green's function can be solved in the frequency domain, yielding Gs(xs|x,ω)f1(ω) and Gr(xr|x,ω)f2(ω). Since the specific method of Taylor expansion is existing technology, it will not be elaborated here.

[0109] In this embodiment, as Figure 4 The diagram shows the Green's function Gs(1000m|x,50Hz)f1(50Hz) derived from the source point S(1000m,50m) to the imaging point, obtained by calculating the positional relationship between the imaging point and the source point as determined in step 1), and by using the recursive region and direction of the one-way wavefield from the source point to the imaging point determined through region decomposition. Here, ω = 50Hz represents the Green's function Gs with a frequency of 50Hz. (Refer to...) Figure 4 The imaging point is mainly located below the seismic source point, and the downward regional decomposition method and recursive direction were indeed selected according to the judgment method in step 1).

[0110] like Figure 5 The diagram shows the Green's function Gr(3000m|x,50Hz)f1(50Hz) derived from the detector point r(3000m,50m) to the imaging point, obtained by calculating the positional relationship between the imaging point and the detector point in step 1), the recursive region of the one-way wavefield from the detector point to the imaging point through region decomposition, and the recursive direction; where ω=50Hz, representing the Green's function Gr with a frequency of 50Hz. (Refer to...) Figure 5 The imaging point is mainly located to the left of the detector point, and the leftward region decomposition method and recursive direction were indeed selected according to the judgment method in step 1).

[0111] As can be seen, compared with the traditional calculation of generalized diffraction superposition in the time domain, this embodiment also implements Green's function solution in the frequency domain. There is no need to convert the generalized diffraction superposition formula used for imaging to the time domain solution. Therefore, there is no sampling rate dispersion problem inherent in the time domain. That is, the frequency domain solution equation of the one-way wave field is unconditionally stable, which makes the above calculation method adaptable to larger spatial grid intervals.

[0112] 3) Substitute the reflection seismic data of the target area, the solution of Green's function from the source point to the imaging point, and the solution of Green's function from the receiver point to the imaging point into the generalized diffraction superposition formula, and integrate along the receiver points and frequencies at different locations to obtain the imaging value of the current imaging point.

[0113] In this embodiment, using Gs(xs|x,ω)f1(ω) and Gr(xr|x,ω)f2(ω) calculated in step 2), and substituting them into the generalized diffraction superposition formula used to calculate the imaging value, along the detector points x at different locations... r By integrating with the frequency ω, the final imaging value I(x) of the front imaging point can be obtained.

[0114] 4) Determine whether imaging of all grid points in the velocity model has been completed; if yes, proceed to step 5); if no, return to step 1) and continue imaging the next grid point.

[0115] 5) Output the imaging values ​​of all grid points in the velocity model as the final pre-stack depth migration imaging result.

[0116] In this embodiment, as Figure 6 The image shown is the final imaging result obtained using the pre-stack depth migration imaging method based on generalized diffraction superposition in this embodiment. It can be seen that, due to the use of Green's function for one-way wavefield solution, which excludes wavefield components such as reflected waves that are not beneficial to imaging, shallow imaging noise is low, and the imaging results can clearly distinguish the location and morphology of the strata (see...). Figure 6 (The area indicated by the middle arrow A); in addition, deep imaging results have been improved, and geological structural energy has been enhanced (see...). Figure 6 (The area indicated by the middle arrow B).

[0117] Figure 7 This is the imaging result obtained by using conventional generalized diffraction stacking imaging followed by Laplacian filtering for noise reduction. It can be seen that, due to the use of the two-way wave equation, wavefield components such as reflected waves, which are not beneficial to imaging, also participate in the imaging process, leading to increased shallow source noise and consequently, a deterioration in the stratigraphic continuity of the shallow imaging effect (see...). Figure 7 (The white noise area indicated by the middle arrow A); Furthermore, the results described above are consistent with those of this embodiment. Figure 6 compared to, Figure 7 Traditional generalized diffraction superposition deep structure imaging suffers from weakened energy and insufficient resolution (see...). Figure 7 (The area indicated by the middle arrow B).

[0118] like Figure 8 The image shows the hardware. Under the computational conditions of a Core™ i7-10700 CPU @ 2.90GHz, with a spatial sampling interval dx = dy = 10m and sampling points nx = 401 and ny = 301, the computational efficiency of the Green's function in the conventional generalized diffraction stacking imaging method, the Green's function in the generalized diffraction stacking imaging method using the single-way wavefield recursion of this embodiment, and the Green's function in the generalized diffraction stacking imaging method using the single-way wavefield recursion after region decomposition of this embodiment is compared. It can be seen that, under the same conditions... Under the given computing hardware conditions, the computational efficiency of the generalized diffraction stacking imaging method based on single-path wavefield recursion in this embodiment is significantly higher than that of the conventional generalized diffraction stacking imaging method. After angular region decomposition, useless computational regions can be removed to a certain extent, and the computational efficiency of Green's function can be further improved, thereby improving the overall imaging efficiency. Overall, the computational efficiency of the pre-stack depth migration imaging method based on generalized diffraction stacking in this embodiment is much higher than that of the conventional generalized diffraction stacking imaging method based on two-path wavefield, with a computational efficiency of approximately 4.23 times that of the conventional method.

[0119] Therefore, the pre-stack depth migration imaging method based on angle domain decomposition and generalized diffraction stacking of one-way wavefields in this embodiment first reads in the reflection seismic data, migration velocity field, seismic wavelet, and wavefield recursion parameters in the program from the common shot domain, and begins the looping and imaging process of imaging points. Then, it calculates the angle of the imaging point relative to the source and receiver points, determines the recursion direction of Green's function, and performs domain decomposition. Afterward, along the determined direction, it uses the one-way wavefield to recursively derive the Green's function from the source point to the imaging point and from the receiver point to the imaging point, respectively. Then, it integrates the calculated Green's function to obtain the imaging value of the imaging point. Finally, the imaging points are looped until imaging of all imaging points is completed. Compared with the two-way wavefield used in conventional generalized diffraction stacking migration, the technical solution of this embodiment uses the one-way wavefield to achieve generalized diffraction stacking migration. While realizing imaging of complex structures, it overcomes the defects of two-way wavefield generalized diffraction stacking and solves the problems of low computational efficiency and high noise in conventional generalized diffraction stacking methods.

[0120] In summary, the pre-stack depth migration imaging method based on generalized diffraction stacking in this embodiment fundamentally changes the two-way wave equation method in generalized diffraction stacking to a one-way wave equation method. By rapidly solving the one-way wave field in the frequency domain, it achieves, for the first time, a generalized diffraction stacking method based on a one-way wave field, overcoming the limitations of the traditional two-way method in generalized diffraction stacking. Furthermore, by decomposing the imaging angle (direction) into regions and determining the recursive region where the imaging point is located, useless computational regions are reduced, further improving computational efficiency. Therefore, it can achieve generalized diffraction stacking migration with both high efficiency and high quality, improving imaging effects at shallow, medium, and deep layers, overcoming the shortcomings of conventional generalized diffraction stacking methods, and further shortening the computation time of Green's function in the generalized diffraction stacking imaging formula through region decomposition, thus resulting in higher overall imaging efficiency.

[0121] This invention has the following characteristics:

[0122] 1) The Green's function required for imaging is calculated using the one-way wave equation, instead of the two-way wave equation used in traditional generalized diffraction superposition. This avoids the inherent defects of the two-way wave equation, ensuring that the imaging results do not contain wavefield components that are not beneficial to imaging, reducing noise in shallow imaging and enhancing the geological structural energy in deep imaging, thus guaranteeing good imaging results. It also reduces the computational cost of the Green's function, speeds up the operation, and improves imaging efficiency. Furthermore, the Green's function is solved in the frequency domain, eliminating the need to convert the generalized diffraction superposition formula used for imaging to the time domain. Therefore, the sampling rate dispersion problem inherent in the time domain does not exist, allowing the calculation method of the Green's function using the one-way wave equation to adapt to larger spatial grid intervals.

[0123] 2) After calculating the angles to determine the relative positions of the imaging point with respect to the source point and the receiver point, the region decomposition is performed to obtain the recursive regions and directions of the one-way wavefield from the source point to the imaging point and from the receiver point to the imaging point. That is, the necessary recursive regions are selected according to the relative positions of the imaging point with respect to the source point and the receiver point, so that when recursively solving the Green's function, it is only necessary to perform the recursion in the necessary regions, without having to perform the recursion in the full velocity model space. While ensuring the recursion effect, the angle region decomposition can remove useless calculation regions to a certain extent and further improve the calculation efficiency of the Green's function, thereby improving the overall imaging efficiency.

[0124] It should be understood that the above-described specific embodiments of the present invention are merely illustrative or explanatory of the principles of the present invention, and do not constitute a limitation thereof.

Claims

1. A method for pre-stack depth migration imaging based on generalized diffraction stack, characterized in that, Read the relevant data of the target area required for imaging, including the velocity model corresponding to the offset velocity field, the location of the receiver point, and the location of the source point, and set grid points as imaging points in the velocity model space at a set interval; Cyclic imaging is performed on each imaging point to obtain pre-stack depth migration imaging results, including: 1) Select a grid point as the current imaging point, and determine the positional relationship between the current imaging point and the source point in the velocity model. The formula for the positional relationship between the current imaging point and the source point is: ； In the formula, A s This indicates the current imaging point position (x, y) relative to the source point position (x, y). s , y s The angle of As; if the value of As is within the range of the first set source point angle, then it is determined that the imaging point is mainly located to the right of the source point. The velocity model is laterally decomposed into left and right sides along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region where the source point +x direction is located. S ; If the value of As falls within the range of the second or third preset source point angle, the imaging point is determined to be mainly located to the left of the source point. The velocity model is then horizontally decomposed into left and right sides along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region in the -x direction of the source point. S ; 2) Based on the recursive region and recursive direction of the single-pass wave field from the source point to the imaging point, and the recursive region and recursive direction of the single-pass wave field from the receiver point to the imaging point, Green's function from the source point to the imaging point and from the receiver point to the imaging point are solved recursively using the single-pass wave. 3) Substitute the reflection seismic data of the target area, the solution of Green's function from the source point to the imaging point, and the solution of Green's function from the receiver point to the imaging point into the generalized diffraction superposition formula, and integrate along the receiver points and frequencies at different locations to obtain the imaging value of the current imaging point. 4) Determine if imaging of all grid points in the velocity model has been completed; if yes, proceed to step 5); if no, return to step 1) and continue imaging the next grid point. 5) Output the imaging values ​​of all grid points in the velocity model as the final pre-stack depth migration imaging result.

2. The generalized diffraction stack-based pre-stack depth migration imaging method of claim 1, wherein, Step 1) also includes: If the value of As falls within the range of the fourth set source point angle, then the imaging point is determined to be mainly located below the source point. The velocity model is longitudinally decomposed into upper and lower parts along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region containing the source point in the +y direction. S ; If the value of As falls within the range of the fifth set source point angle, then the imaging point is determined to be mainly located above the source point. The velocity model is longitudinally decomposed into upper and lower parts along the source point, and Green's function G from the source point to the imaging point is recursively solved only along the region in the -y direction of the source point. S .

3. The generalized diffraction stack-based pre-stack depth migration imaging method according to any one of claims 1 or 2, characterized in that, The method for selecting the recursive region and recursive direction of the one-way wavefield from the detector point to the imaging point is as follows: In the velocity model, the positional relationship between the current imaging point and the detector point is determined. The formula for the positional relationship between the current imaging point and the detector point is: ； In the formula, A r This indicates the current imaging point position (x, y) relative to the detector point position (x, y). r , y r Angle of A; if A r If the value is within the range of the first set detector angle, then it is determined that the imaging point is mainly located to the right of the detector. The velocity model is decomposed laterally into left and right sides along the detector. The Green's function Gr from the detector to the imaging point is solved only along the region where the detector +x direction is located. If A r If the value is within the range of the second set detector angle or the third set detector angle, then the imaging point is determined to be mainly located to the left of the detector; the velocity model is decomposed laterally into left and right sides along the detector, and Green's function Gr from the detector to the imaging point is solved recursively only along the region where the detector is located in the -x direction; If A r If the value is within the range of the fourth set detector angle, then it is determined that the imaging point is mainly located below the detector. The velocity model is longitudinally decomposed into upper and lower parts along the detector. The Green's function Gr from the detector to the imaging point is solved only along the region where the detector +y direction is located. If A r If the value is within the range of the fifth set detector angle, then it is determined that the imaging point is mainly located above the detector. The velocity model is longitudinally decomposed into upper and lower parts along the detector, and Green's function Gr from the detector to the imaging point is solved only along the region where the detector is located in the -y direction.

4. The generalized diffraction stack-based pre-stack depth migration imaging method according to any one of claims 1-2, wherein, The relevant data for the target region required for imaging also include the source function in the frequency domain of the target region, seismic data, and seismic wavelets; the generalized diffraction superposition formula is: ； Where I(x) represents the imaging result, Gs and Gr are Green's functions from the source point and receiver point to the imaging point, respectively, d() represents the seismic data of the target area, xs and xr are the locations of the source point and receiver point, respectively, and x is the grid position. For frequency, Indicates complex conjugation; In the generalized diffraction superposition formula: ； wherein, is a source function of the frequency domain of the target region, and are respectively a first decomposition item and a second decomposition item obtained by decomposing the seismic wavelet corresponding to the source.

5. The generalized diffraction stack-based prestack depth migration imaging method of claim 4, wherein, Based on the recursive region and direction of the one-way wavefield from the source point to the imaging point, and the recursive region and direction of the one-way wavefield from the receiver point to the imaging point, the Green's function from the source point to the imaging point and from the receiver point to the imaging point is solved recursively using the one-way wavefield, respectively: the first decomposition term of the seismic wavelet corresponding to the source For the boundary condition, the formula of Green's function Gsfrom the source point to the imaging point is recursively solved along the x-axis direction as follows: ； In this formula, a "+" indicates recursively solving Green's function along the positive x-axis, and a "-" indicates recursively solving Green's function along the negative x-axis; i 2 =-1, Gs is Green's function from the source point to the imaging point. It is a one-way recursive operator for approximating the wave field along the x-direction; The formula for recursively solving Green's function Gs from the seismic source point to the imaging point along the y-axis is as follows: ； In this formula, a + sign indicates that the Green's function is solved recursively along the positive y-axis, and a - sign indicates that the Green's function is solved recursively along the negative y-axis; i 2 =-1, Gs is Green's function from the source point to the imaging point. This is a one-way recursive operator for approximating the wave field along the y-direction; The second decomposition term of the seismic wavelet corresponding to the source For the boundary condition, the formula for recursively solving the Green's function Gr from the receiver point to the imaging point along the x-axis direction is as follows: ； In this formula, a "+" indicates recursively solving Green's function along the positive x-axis, and a "-" indicates recursively solving Green's function along the negative x-axis; i 2 =-1, Gr is Green's function from the detector point to the imaging point. It is a one-way recursive operator for approximating the wave field along the x-direction; The formula for recursively solving Green's function Gr from the receiver point to the imaging point along the y-axis is as follows: ； In this formula, a + sign indicates that the Green's function is solved recursively along the positive y-axis, and a - sign indicates that the Green's function is solved recursively along the negative y-axis; i 2 =-1, Gr is Green's function from the detector point to the imaging point. This is a one-way recursive operator for approximating the wave field along the y-direction; Approximate one-way propagator operators for wavefield expansion along x, y directions are represented by: ； ； In the above formula, c represents the background velocity; v(x, y) represents the velocity field value; x and y are respectively the coordinates of the imaging point in the horizontal direction and the depth direction; is the frequency.

Images