CN105652320A - Reverse time migration imaging method and apparatus - Google Patents

Abstract

The invention provides a reverse time migration imaging method and device. The method includes the following steps that: a shot of data are obtained; a depth domain velocity field is transformed to pyramid grids from regular grids; the coordinates of anisotropy parameters delta and epsilon, a symmetry axis tile angle parameter theta and a symmetry axis azimuth parameter phi are transformed to the pyramid grids from the regular grids; finite difference wave field simulation is performed on a wave field generated by excitation by using a coupling second-order partial differential equation of the pyramid grids, so that a wave field corresponding to a shot point and a wave field corresponding to a detection point can be obtained; imaging is performed on the wave field corresponding to the shot point and the wave field corresponding to the detection point by using a cross-correlation imaging condition; an imaging result is interpolated to the regular grids, the imaging result which has been interpolated to the regular grids is adopted as a single-shot imaging result of the shot point. With the method and device of the invention adopted, the technique problem that changing grids cannot be utilized to perform imaging in a horizontal direction and a vertical direction simultaneously can be solved, and changing grids can be utilized to perform imaging in a horizontal direction and a vertical direction simultaneously in a horizontal direction and a vertical direction.

Description

Reverse time migration imaging method and device

Technical Field

The invention relates to the technical field of geological exploration, in particular to a reverse time migration imaging method and device.

Background

The Reverse Time Migration (RTM) technology is based on a full acoustic wave equation, completely follows the physical law of fluctuation in an acoustic wave medium, and completely simulates various wave phenomena in the process of seismic wave field propagation with high precision, so that the RTM technology is more suitable for the seismic wave imaging problem of complex-structure areas with violent speed change, and especially when the wave field has the phenomena of rotating waves and prism reflection, the RTM technology can show the advantage over one-way wave equation migration.

Since Reverse Time Migration (RTM) simulates the propagation of seismic waves in the time-space domain by means of a difference method, this approach can be easily and naturally generalized to the problem of imaging non-uniform anisotropic media. These important features of Reverse Time Migration (RTM) make it the first choice for imaging complex regions in the seismic exploration industry, an important leading-edge technique for imaging complex structures. Research on Reverse Time Migration (RTM) focuses on improving the work efficiency of Reverse Time Migration (RTM) while ensuring the imaging accuracy.

Wave field simulation for acoustic wave and elastic wave equations is the basis of wave equation-based migration imaging methods, and commonly used numerical calculation methods include: finite element methods and finite difference methods, where finite element methods have great flexibility in dealing with irregular meshes, but suffer from the lack of capability to be applied to higher order algorithms, particularly to account for frequency dependent dispersion. In contrast, the finite difference can be well applied to high-order algorithms, and algorithm optimization methods for the finite difference are many (for example, variable grids are adopted in the depth Z direction), and the method can adopt larger Z direction grids for places with large deep layer speed without reducing the high-frequency precision of simulation

With the widespread use of Reverse Time Migration (RTM) and other techniques, there is an increasing demand for improving the computational efficiency of wavefield simulation, and therefore, it is considered whether a varying grid can be used in the horizontal direction.

However, no effective solution has been proposed at present how to perform imaging in both the horizontal direction (X \ Y) and the vertical direction (Z) with varying grids.

Disclosure of Invention

The embodiment of the invention provides a reverse time migration imaging method, which aims to achieve the purpose of imaging by simultaneously adopting a changed grid in the horizontal direction and the vertical direction, and comprises the following steps:

s1: acquiring shot data;

s2: reading a depth domain velocity field corresponding to an imaging space of a shot point of the shot data, and converting the depth domain velocity field from a regular grid to a pyramid grid;

s3: reading anisotropic parameter sums corresponding to the imaging space, and converting the coordinates of the anisotropic parameter sums from a regular grid to a pyramid grid;

s4: reading a symmetrical axis inclination angle parameter theta and a symmetrical axis azimuth angle parameter phi corresponding to the imaging space, and converting the coordinates of the symmetrical axis inclination angle parameter theta and the symmetrical axis azimuth angle parameter phi from a regular grid to a pyramid grid;

s5: exciting the shot point by placing a wavelet, and performing finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the shot point;

s6: performing finite difference wave field simulation on the shot data through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the detection point;

s7: applying a cross-correlation imaging condition to image a wave field corresponding to the shot point and a wave field corresponding to the wave detection point;

s8: interpolating the imaging result to a regular grid, and taking the imaging result interpolated to the regular grid as the single shot imaging result of the shot point;

repeating the steps S1 to S8 on a plurality of guns to obtain the imaging results of the single guns, and superposing the imaging results of the single guns to obtain the reverse time migration imaging section.

In one embodiment, the coordinates are transformed from a regular grid to a pyramid grid according to the following formula:

$\left\{\begin{array}{c}{x}_p=\frac{x}{1+\mathrm{\α}z}\\ y_p=\frac{y}{1+\mathrm{\α}z}\\ z_p=\frac{1}{\mathrm{\α}}ln(1+\mathrm{\α}z)\end{array}\right.$

wherein (x, y, z) represents the coordinate system of a regular grid, (x)_p,y_p,z_p) A coordinate system representing a pyramid grid, α coordinate transform coefficients.

In one embodiment, the coupled second order partial differential equation for the pyramid grid is:

$\left\{\begin{array}{c}\frac{1}{v_{pz}^2}\frac{{\∂}^{2}p}{\∂t^2}=(1+2\mathrm{\ϵ})H_{2}p+H_{1}q+\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_1(p-q)\\ \frac{1}{v_{pz}^2}\frac{{\∂}^{2}q}{\∂t^2}=(1+2\mathrm{\δ})H_{2}p+H_{1}q-\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_2(p-q)\end{array}\right.$

wherein,

$\begin{array}{c}{H}_1=\frac{e^{-2{\mathrm{\αz}}_p}}{1+{\left(xdip\right)}^2+{\left(ydip\right)}^2}\left(\begin{array}{c}{(xdip-{\mathrm{\αx}}_p)}^2\frac{{\∂}^2}{\∂x_p^2}+{(ydip-{\mathrm{\αy}}_p)}^2\frac{{\∂}^2}{\∂y_p^2}\\ +\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +2(xdip-{\mathrm{\αx}}_p)(ydip-{\mathrm{\αy}}_p)\frac{{\∂}^2}{\∂x_p\∂y_p}\\ -2(xdip-{\mathrm{\αx}}_p)(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})\\ -2(ydip-{\mathrm{\αy}}_p)(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\end{array}\right)\\ H_2=(1+{\left({\mathrm{\αx}}_p\right)}^2)\frac{{\∂}^2}{\∂x_p^2}+(1+{\left({\mathrm{\αy}}_p\right)}^2)\frac{{\∂}^2}{\∂y_p^2}+\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +{\mathrm{\αx}}_p(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})+2{\mathrm{\αy}}_p(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\\ +2{\mathrm{\α}}^2{x}_p{y}_p\frac{{\∂}^2}{\∂x_p\∂y_p}-H_1\end{array}$

wherein xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ),p denotes the P wave wavefield, q denotes the q wave wavefield, v_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

In one embodiment, the finite difference wave field simulation of the wave field generated by excitation through the coupling second-order partial differential equation of the pyramid grid to obtain the wave field corresponding to the shot point comprises the following steps: performing wave field compression on the wave field corresponding to the shot point according to a preset time interval, and storing a compressed result in a local disk;

the imaging of the wave field corresponding to the shot point and the wave field corresponding to the wave detection point by applying the cross-correlation imaging condition comprises the following steps: decompressing and reading the wave field corresponding to the shot point from the local disk; and imaging the wave field corresponding to the wave detection point and the wave field corresponding to the shot point which is decompressed by applying a cross-correlation imaging condition.

In one embodiment, interpolating the imaging results to a regular grid comprises: interpolating the imaging results into the regular grid in the X-direction, Y-direction and Z-direction.

The embodiment of the invention also provides a reverse time migration imaging device, which is used for achieving the purpose of imaging in the horizontal direction and the vertical direction by simultaneously adopting the changed grids, and comprises the following components:

the acquisition module is used for acquiring shot data;

the first transformation module is used for reading a depth domain velocity field corresponding to an imaging space of a shot point of the shot data and transforming the depth domain velocity field from a regular grid to a pyramid grid;

the second transformation module is used for reading anisotropic parameter sums corresponding to the imaging space and transforming the coordinates of the anisotropic parameter sums from a regular grid to a pyramid grid;

the third transformation module is used for reading a symmetry axis inclination angle parameter theta and a symmetry axis azimuth angle parameter phi corresponding to the imaging space and transforming the coordinates of the symmetry axis inclination angle parameter theta and the symmetry axis azimuth angle parameter phi to a pyramid grid from a regular grid;

the shot point difference simulation module is used for exciting the shot point by placing a wavelet and carrying out finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the shot point;

the wave detection point difference simulation module is used for carrying out finite difference wave field simulation on the shot data through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the wave detection point;

the imaging module is used for imaging the wave field corresponding to the shot point and the wave field corresponding to the wave detection point by applying a cross-correlation imaging condition;

and the interpolation module is used for interpolating the imaging result to the regular grid and taking the imaging result interpolated to the regular grid as the single shot imaging result of the shot point.

In one embodiment, the first transformation module and the second transformation module are specifically configured to transform the coordinates from a regular grid to a pyramid grid according to the following formula:

$\left\{\begin{array}{c}{x}_p=\frac{x}{1+\mathrm{\α}z}\\ y_p=\frac{y}{1+\mathrm{\α}z}\\ z_p=\frac{1}{\mathrm{\α}}ln(1+\mathrm{\α}z)\end{array}\right.$

wherein (x, y, z) represents the coordinate system of a regular grid, (x)_p,y_p,z_p) A coordinate system representing a pyramid grid, α coordinate transform coefficients.

In one embodiment, the coupled second order partial differential equation for the pyramid grid is:

$\left\{\begin{array}{c}\frac{1}{v_{pz}^2}\frac{{\∂}^{2}p}{\∂t^2}=(1+2\mathrm{\ϵ})H_{2}p+H_{1}q+\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_1(p-q)\\ \frac{1}{v_{pz}^2}\frac{{\∂}^{2}q}{\∂t^2}=(1+2\mathrm{\δ})H_{2}p+H_{1}q-\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_2(p-q)\end{array}\right.$

wherein,

$\begin{array}{c}{H}_1=\frac{e^{-2{\mathrm{\αz}}_p}}{1+{\left(xdip\right)}^2+{\left(ydip\right)}^2}\left(\begin{array}{c}{(xdip-{\mathrm{\αx}}_p)}^2\frac{{\∂}^2}{\∂x_p^2}+{(ydip-{\mathrm{\αy}}_p)}^2\frac{{\∂}^2}{\∂y_p^2}\\ +\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +2(xdip-{\mathrm{\αx}}_p)(ydip-{\mathrm{\αy}}_p)\frac{{\∂}^2}{\∂x_p\∂y_p}\\ -2(xdip-{\mathrm{\αx}}_p)(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})\\ -2(ydip-{\mathrm{\αy}}_p)(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\end{array}\right)\\ H_2=(1+{\left({\mathrm{\αx}}_p\right)}^2)\frac{{\∂}^2}{\∂x_p^2}+(1+{\left({\mathrm{\αy}}_p\right)}^2)\frac{{\∂}^2}{\∂y_p^2}+\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +{\mathrm{\αx}}_p(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})+2{\mathrm{\αy}}_p(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\\ +2{\mathrm{\α}}^2{x}_p{y}_p\frac{{\∂}^2}{\∂x_p\∂y_p}-H_1\end{array}$

wherein xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ),p denotes the P wave wavefield, q denotes the q wave wavefield, v_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

In one embodiment, the shot point difference simulation module further performs wave field compression on the wave field corresponding to the shot point according to a preset time interval in the process of performing finite difference wave field simulation on the wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain the wave field corresponding to the shot point, and stores the compressed result in a local disk; the imaging module is specifically used for decompressing and reading the wave field corresponding to the shot point from the local disk; and imaging the wave field corresponding to the wave detection point and the wave field corresponding to the shot point which is decompressed by applying a cross-correlation imaging condition.

In one embodiment, the interpolation module is specifically configured to interpolate the imaging results into the regular grid in an X direction, a Y direction, and a Z direction.

In the embodiment, the seismic data are converted from the regular grid to the pyramid grid, reverse time migration imaging is carried out through a TTI seismic anisotropic medium of the pyramid grid, and differential solution of a wave equation is realized based on a coupling second-order partial differential equation of the pyramid grid, so that the calculation efficiency can be remarkably improved, the memory consumption can be reduced, and the problem of three-dimensional complex structure imaging with rapid speed change can be finally solved on the premise of not influencing the imaging effect. By means of the method, the technical problem that the imaging cannot be carried out by simultaneously adopting the changed grids in the horizontal direction and the vertical direction in the prior art is solved, and the purpose of simultaneously adopting the changed grids in the horizontal direction and the vertical direction to carry out imaging is achieved.

Drawings

The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this application, illustrate embodiment(s) of the invention and together with the description serve to explain the principles of the invention. In the drawings:

FIG. 1 is a flow chart of a reverse time migration imaging method according to an embodiment of the present invention;

FIG. 2 is a diagram illustrating the relationship between an original grid and a pyramid grid according to an embodiment of the present invention;

FIG. 3 is a time-reversed migration shot wavefield snapshot comparison of a pyramid grid and a regular grid in accordance with an embodiment of the present invention;

FIG. 4 is a time-reversed migration of a pyramid grid compared to a regular grid for a snapshot of a wave-field of a receiver in accordance with an embodiment of the present invention;

FIG. 5 is a pyramid grid TTI medium reverse time migration single shot migration result according to an embodiment of the present invention;

FIG. 6 is a diagram illustrating interpolation of single shot migration results from a pyramid grid to a regular grid in accordance with an embodiment of the present invention;

FIG. 7 is a schematic diagram of a pyramid grid and a regular grid comparing the TTI medium anisotropic reverse time migration 16 gun migration superposition result according to an embodiment of the present invention;

fig. 8 is a block diagram of a reverse time shift imaging device according to an embodiment of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail with reference to the following embodiments and accompanying drawings. The exemplary embodiments and descriptions of the present invention are provided to explain the present invention, but not to limit the present invention.

The inventors consider that imaging can be performed using a pyramid Grid (P-Grid), thereby achieving the purpose of imaging using a varied Grid in both the horizontal and vertical directions.

Specifically, a reverse time migration imaging method is provided, as shown in fig. 1, including the following steps:

s1: acquiring a shot of data, for example, a shot of data may be retrieved from a task list in preparation for Reverse Time Migration (RTM) of the shot;

s2: and reading a depth domain velocity field corresponding to the imaging space of the shot point of the shot data, and converting the depth domain velocity field from a regular grid to a pyramid grid, namely, reading the depth domain velocity field corresponding to the imaging space, and converting the coordinates to the pyramid grid.

Let u (x, y, z) be the spatial function in the original grid (i.e. Cartesian coordinate system) and (x) be the coordinate system of the new grid_p,y_p,z_p) And the conversion between the two needs to satisfy:

$\left\{\begin{array}{c}{x}_p=\frac{x}{1+\mathrm{\α}z}\\ y_p=\frac{y}{1+\mathrm{\α}z}\\ z_p=\frac{1}{\mathrm{\α}}ln(1+\mathrm{\α}z)\end{array}\right.$

wherein (x, y, z) represents the coordinate system of a regular grid, (x)_p,y_p,z_p) A coordinate system representing a pyramid grid, α coordinate transform coefficients.

FIG. 2 is a schematic diagram of coordinates of an original grid and a pyramid grid, wherein the dark color is the pyramid grid (x)_p,z_p) The light color is the original grid (x, z), and the velocity field can be converted from the original grid (i.e. regular grid) to the pyramid grid by the coordinate transformation formula.

S3: and reading anisotropic parameter sums corresponding to the imaging space, and transforming the coordinates of the anisotropic parameter sums from a regular grid to a pyramid grid.

S4: and reading a symmetrical axis inclination angle parameter theta and a symmetrical axis azimuth angle parameter phi corresponding to the imaging space, and transforming the coordinates of the symmetrical axis inclination angle parameter theta and the symmetrical axis azimuth angle parameter phi from a regular grid to a pyramid grid.

S5: exciting the shot point by placing a wavelet, and performing finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the shot point;

during specific implementation, in the process of performing finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of a pyramid grid to obtain a wave field corresponding to a shot point, wave field compression can be performed on the wave field corresponding to the shot point according to a preset time interval, and a compressed result is stored in a local disk, namely, a wavelet is placed at the shot point position corresponding to the shot point, finite difference wave field simulation is performed, and meanwhile, the wave field compression is performed on the wave field at a certain time interval and is stored in the local disk.

By the above original grid (x, z) and pyramid grid (x)_p,z_p) The relation between the first and second partial derivatives of the two grids can be obtained by the coordinate transformation formula as follows:

$\frac{\∂u}{\∂x}=\frac{\∂u}{\∂x_p}{e}^{-{\mathrm{\αz}}_p}$

$\frac{\∂u}{\∂y}=\frac{\∂u}{\∂y_p}{e}^{-{\mathrm{\αz}}_p}$

$\frac{\∂u}{\∂z}=(-{\mathrm{\αx}}_p\frac{\∂u}{\∂x_p}-{\mathrm{\αy}}_p\frac{\∂u}{\∂y_p}+\frac{\∂u}{\∂z_p})e^{-{\mathrm{\αz}}_p}$

$\frac{{\∂}^{2}u}{\∂x^2}=\frac{{\∂}^{2}u}{\∂x_p^2}{e}^{-2{\mathrm{\αz}}_p}$

$\frac{{\∂}^{2}u}{\∂y^2}=\frac{{\∂}^{2}u}{\∂y_p^2}{e}^{-2{\mathrm{\αz}}_p}$

$\frac{{\∂}^{2}u}{\∂z^2}=e^{-2{\mathrm{\αz}}_p}\left(\begin{array}{c}{\mathrm{\α}}^2{x}_p^2\frac{{\∂}^{2}u}{\∂x_p^2}+{\mathrm{\α}}^2{y}_p^2\frac{{\∂}^{2}u}{\∂y_p^2}+\frac{{\∂}^{2}u}{\∂z_p^2}\\ +2{\mathrm{\α}}^2{x}_p\frac{\∂u}{\∂x_p}+2{\mathrm{\α}}^2{y}_p\frac{\∂u}{\∂y_p}-\mathrm{\α}\frac{\∂u}{\∂z_p}\\ +2{\mathrm{\α}}^2{x}_p{y}_p\frac{{\∂}^{2}u}{\∂x_p\∂y_p}-2{\mathrm{\αx}}_p\frac{{\∂}^{2}u}{\∂x_p\∂z_p}-2{\mathrm{\αy}}_p\frac{{\∂}^{2}u}{\∂y_p\∂z_p}\end{array}\right)$

$\frac{{\∂}^{2}u}{\∂x\∂y}=e^{-2{\mathrm{\αz}}_p}\frac{{\∂}^{2}u}{\∂x_p\∂y_p}$

$\frac{{\∂}^{2}u}{\∂x\∂z}=e^{-2{\mathrm{\αz}}_p}(\frac{{\∂}^{2}u}{\∂x_p\∂z_p}-{\mathrm{\αx}}_p\frac{{\∂}^{2}u}{\∂x_p^2}-\mathrm{\α}\frac{\∂u}{\∂x_p}-{\mathrm{\αy}}_p\frac{{\∂}^{2}u}{\∂x_p\∂y_p})$

$\frac{{\∂}^{2}u}{\∂y\∂z}=e^{-2{\mathrm{\αz}}_p}(\frac{{\∂}^{2}u}{\∂y_p\∂z_p}-{\mathrm{\αy}}_p\frac{{\∂}^{2}u}{\∂y_p^2}-\mathrm{\α}\frac{\∂u}{\∂y_p}-{\mathrm{\αx}}_p\frac{{\∂}^{2}u}{\∂x_p\∂y_p})$

through the derivation, a first-order partial derivative, a second-order partial derivative and a second-order mixed partial derivative in the X direction, the Y direction and the Z direction can be obtained, and based on the partial derivative equations, an acoustic wave equation of the TTI anisotropic medium under the pyramid Grid (P-Grid) can be derived, that is, a coupling second-order partial differential equation of the pyramid Grid is expressed as:

$\left\{\begin{array}{c}\frac{1}{v_{pz}^2}\frac{{\∂}^{2}p}{\∂t^2}=(1+2\mathrm{\ϵ})H_{2}p+H_{1}q+\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_1(p-q)\\ \frac{1}{v_{pz}^2}\frac{{\∂}^{2}q}{\∂t^2}=(1+2\mathrm{\δ})H_{2}p+H_{1}q-\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_2(p-q)\end{array}\right.$

wherein:

$\left\{\begin{array}{c}{H}_1=\frac{1}{1+{\left(xdip\right)}^2+{\left(ydip\right)}^2}\left(\begin{array}{c}{\left(xdip\right)}^2\frac{{\∂}^2}{\∂x^2}+{\left(ydip\right)}^2\frac{{\∂}^2}{\∂y^2}+\frac{{\∂}^2}{\∂z^2}\\ +2\left(xdip\right)\left(ydip\right)\frac{{\∂}^2}{\∂x\∂y}\\ +2\left(xdip\right)\frac{{\∂}^2}{\∂x\∂z}+2\left(ydip\right)\frac{{\∂}^2}{\∂y\∂z}\end{array}\right)\\ H_2=\frac{{\∂}^2}{\∂x^2}+\frac{{\∂}^2}{\∂y^2}+\frac{{\∂}^2}{\∂y^2}+\frac{{\∂}^2}{\∂z^2}-H_1\end{array}\right.$

where xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ)), P represents a P-wave wavefield, q represents a q-wave wavefield, and v represents a_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

Substituting the derived partial derivative formula under the pyramid Grid (P-Grid), and further obtaining a coupling second-order partial differential equation of the pyramid Grid, which is finally expressed as:

$\begin{array}{c}{H}_1=\frac{e^{-2{\mathrm{\αz}}_p}}{1+{\left(xdip\right)}^2+{\left(ydip\right)}^2}\left(\begin{array}{c}{(xdip-{\mathrm{\αx}}_p)}^2\frac{{\∂}^2}{\∂x_p^2}+{(ydip-{\mathrm{\αy}}_p)}^2\frac{{\∂}^2}{\∂y_p^2}\\ +\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +2(xdip-{\mathrm{\αx}}_p)(ydip-{\mathrm{\αy}}_p)\frac{{\∂}^2}{\∂x_p\∂y_p}\\ -2(xdip-{\mathrm{\αx}}_p)(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})\\ -2(ydip-{\mathrm{\αy}}_p)(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\end{array}\right)\\ H_2=(1+{\left({\mathrm{\αx}}_p\right)}^2)\frac{{\∂}^2}{\∂x_p^2}+(1+{\left({\mathrm{\αy}}_p\right)}^2)\frac{{\∂}^2}{\∂y_p^2}+\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +{\mathrm{\αx}}_p(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})+2{\mathrm{\αy}}_p(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\\ +2{\mathrm{\α}}^2{x}_p{y}_p\frac{{\∂}^2}{\∂x_p\∂y_p}-H_1\end{array}$

wherein xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ),p denotes the P wave wavefield, q denotes the q wave wavefield, v_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

It can be seen that the above formula seems to be complicated, but in practice the number of partial derivatives to be calculated is the same, that is, the calculation amount is not increased significantly because the formula is complicated, but since the number of grid points required for the pyramid grid is significantly reduced compared to the original grid, the calculation efficiency can be significantly improved and the memory consumption can be reduced.

Fig. 3 to 7 show three-dimensional SEG \ EAGE salt dome model TTI medium anisotropy Reverse Time Migration (RTM) tests, where fig. 3 shows a comparison between pyramid Grid (P-Grid) and regular Grid Reverse Time Migration (RTM) shot point (Source) wave field snapshots, fig. 4 shows a comparison between pyramid Grid (P-Grid) and regular Grid Reverse Time Migration (RTM) wave field snapshots, fig. 5 shows a comparison between pyramid Grid (P-Grid) and regular Grid Reverse Time Migration (RTM) single shot migration result (depth slice), fig. 6 shows an interpolation of the above single shot migration result from pyramid Grid (P-Grid) to regular Grid, which will be saved to a local disk, and fig. 7 shows a comparison between TTI medium anisotropy Reverse Time Migration (RTM)16 shot migration superposition result pyramid Grid and regular Grid. It can be seen from fig. 3 to 7 that the imaging effect is basically equivalent, but the time required for computing one shot by pyramid grid reverse time migration is 50 minutes and the memory is 9.7G, while the regular grid reverse time migration is 105 minutes and the memory is 18.8G, so that the computing efficiency is obviously improved, and the memory consumption is reduced.

FIG. 3 is a shot wave field snapshot (comparison of pyramid grid and regular grid) of three-dimensional SEG \ EAGE salt dome model TTI medium anisotropy Reverse Time Migration (RTM) medium anisotropy test, wherein 416, 910 grid points are respectively used in X \ Y \ Z three directions of the regular grid, and 256, 698 grid points are only used in the pyramid grid, and shot wave fields at certain time intervals are compressed and then stored in a local disk.

S6: and performing finite difference wave field simulation on the shot data through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to a wave detection point, namely reading shot gather data to perform wave field simulation of a wave detection point, and performing differential wave field simulation by using the partial differential equation same as the shot point to obtain a wave field corresponding to the wave detection point.

S7: applying a cross-correlation imaging condition to image a wave field corresponding to the shot point and a wave field corresponding to the wave detection point;

s8: interpolating the imaging result to a regular grid, and taking the imaging result interpolated to the regular grid as the single shot imaging result of the shot point;

that is, the shot wavefield may be decompressed and read during the geophone wavefield simulation, the cross-correlation imaging conditions are applied for imaging, the imaging results are interpolated to a regular grid, and the single-shot migration results for that shot are stored in a local disk.

FIG. 4 is a three-dimensional SEG \ EAGE salt dome model TTI dielectric anisotropy Reverse Time Migration (RTM) test demodulator probe wave field snapshot (comparison of pyramid grid and regular grid), wherein 416, 910 grid points are respectively used in X \ Y \ Z three directions of the regular grid, and 256, 698 grid points are used in the pyramid grid.

And (3) imaging by applying a cross-correlation imaging condition to obtain a Reverse Time Migration (RTM) migration result of the single shot, interpolating the migration result of the single shot to a regular grid in the X \ Y direction as shown in figure 5, and storing the result in a local disk as shown in figure 6.

Repeating the above S1-S8 for a plurality of guns to obtain the imaging results of the single guns, and superposing the imaging results of the single guns to obtain a reverse time migration imaging section, namely superposing all the imaging results of the single guns to form a final Reverse Time Migration (RTM) imaging section.

As shown in fig. 7, when the TTI medium anisotropic Reverse Time Migration (RTM) 16-shot migration superposition result is compared with the regular grid, it can be seen that the imaging effect is basically equivalent, but the time required for computing one shot by the pyramid grid reverse time migration is 50 minutes and the memory is 9.7G, while the time required for computing the regular grid reverse time migration is 105 minutes and the memory is 18.8G, which significantly improves the computing efficiency and reduces the memory consumption.

In the above example, a TTI seismic anisotropic medium reverse time migration imaging method based on a pyramid grid is provided, which can solve the imaging problem of a three-dimensional complex structure, because the differential solution of a wave equation is realized based on a coupling second-order partial differential equation of the pyramid grid, the calculation efficiency can be significantly improved, the memory consumption can be reduced, and the imaging problem of the three-dimensional complex structure with rapidly changing speed can be finally solved without affecting the imaging effect, and the method is easy to implement, and is suitable for the development of Reverse Time Migration (RTM) commercial software and the requirement of industrial production.

Based on the same inventive concept, the embodiment of the present invention further provides a reverse time shift imaging device, as described in the following embodiments. Because the principle of solving the problem of the reverse time migration imaging device is similar to that of the reverse time migration imaging method, the implementation of the reverse time migration imaging device can refer to the implementation of the reverse time migration imaging method, and repeated details are not repeated. As used hereinafter, the term "unit" or "module" may be a combination of software and/or hardware that implements a predetermined function. Although the means described in the embodiments below are preferably implemented in software, an implementation in hardware, or a combination of software and hardware is also possible and contemplated. Fig. 8 is a block diagram of a reverse time shift imaging apparatus according to an embodiment of the present invention, as shown in fig. 8, including: the structure of the system is described below by an obtaining module 801, a first transformation module 802, a second transformation module 803, a third transformation module 804, a shot point difference simulation module 805, a demodulator probe difference simulation module 806, an imaging module 807, and an interpolation module 808.

An obtaining module 801, configured to obtain shot data;

a first transformation module 802, configured to read a depth domain velocity field corresponding to an imaging space of a shot of the shot data, and transform the depth domain velocity field from a regular grid to a pyramid grid;

a second transformation module 803, configured to read an anisotropic parameter sum corresponding to the imaging space, and transform coordinates of the anisotropic parameter sum from a regular grid to a pyramid grid;

a third transformation module 804, configured to read a symmetry axis inclination angle parameter θ and a symmetry axis azimuth angle parameter Φ corresponding to the imaging space, and transform coordinates of the symmetry axis inclination angle parameter θ and the symmetry axis azimuth angle parameter Φ from a regular grid to a pyramid grid;

a shot point difference simulation module 805, configured to place a wavelet to excite the shot point, and perform finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of a pyramid grid to obtain a wave field corresponding to the shot point;

a demodulator probe difference simulation module 806, configured to perform finite difference wave field simulation on the shot data through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the demodulator probe;

an imaging module 807 for applying cross-correlation imaging conditions to image the wavefield corresponding to the shot point and the wavefield corresponding to the demodulator probe;

and the interpolation module 808 is configured to interpolate the imaging result to a regular grid, and use the imaging result interpolated to the regular grid as the single shot imaging result of the shot point.

In one embodiment, the first transformation module 802 and the second transformation module 803 may be specifically configured to transform coordinates from a regular grid to a pyramid grid according to the following formula:

$\left\{\begin{array}{c}{x}_p=\frac{x}{1+\mathrm{\α}z}\\ y_p=\frac{y}{1+\mathrm{\α}z}\\ z_p=\frac{1}{\mathrm{\α}}ln(1+\mathrm{\α}z)\end{array}\right.$

wherein (x, y, z) represents the coordinate system of a regular grid, (x)_p,y_p,z_p) A coordinate system representing a pyramid grid, α coordinate transform coefficients.

In one embodiment, the coupled second order partial differential equation for the pyramid grid is:

$\left\{\begin{array}{c}\frac{1}{v_{pz}^2}\frac{{\∂}^{2}p}{\∂t^2}=(1+2\mathrm{\ϵ})H_{2}p+H_{1}q+\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_1(p-q)\\ \frac{1}{v_{pz}^2}\frac{{\∂}^{2}q}{\∂t^2}=(1+2\mathrm{\δ})H_{2}p+H_{1}q-\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_2(p-q)\end{array}\right.$

wherein,

$\begin{array}{c}{H}_1=\frac{e^{-2{\mathrm{\αz}}_p}}{1+{\left(xdip\right)}^2+{\left(ydip\right)}^2}\left(\begin{array}{c}{(xdip-{\mathrm{\αx}}_p)}^2\frac{{\∂}^2}{\∂x_p^2}+{(ydip-{\mathrm{\αy}}_p)}^2\frac{{\∂}^2}{\∂y_p^2}\\ +\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +2(xdip-{\mathrm{\αx}}_p)(ydip-{\mathrm{\αy}}_p)\frac{{\∂}^2}{\∂x_p\∂y_p}\\ -2(xdip-{\mathrm{\αx}}_p)(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})\\ -2(ydip-{\mathrm{\αy}}_p)(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\end{array}\right)\\ H_2=(1+{\left({\mathrm{\αx}}_p\right)}^2)\frac{{\∂}^2}{\∂x_p^2}+(1+{\left({\mathrm{\αy}}_p\right)}^2)\frac{{\∂}^2}{\∂y_p^2}+\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +{\mathrm{\αx}}_p(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p})+2{\mathrm{\αy}}_p(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p})\\ +2{\mathrm{\α}}^2{x}_p{y}_p\frac{{\∂}^2}{\∂x_p\∂y_p}-H_1\end{array}$

wherein xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ),p denotes the P wave wavefield, q denotes the q wave wavefield, v_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

In one embodiment, the shot point difference simulation module 805 further performs wave field compression on the wave field corresponding to the shot point according to a predetermined time interval in the process of performing finite difference wave field simulation on the wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain the wave field corresponding to the shot point, and stores the compressed result in the local disk; the imaging module 807 may specifically be configured to decompress and read the wavefield corresponding to the shot from the local disk; and imaging the wave field corresponding to the wave detection point and the wave field corresponding to the shot point which is decompressed by applying a cross-correlation imaging condition.

In one embodiment, the interpolation module 808 may be specifically configured to interpolate the imaging results into the regular grid in the X direction, the Y direction, and the Z direction.

From the above description, it can be seen that the embodiments of the present invention achieve the following technical effects: in the embodiment, the seismic data are converted from the regular grid to the pyramid grid, reverse time migration imaging is carried out through a TTI seismic anisotropic medium of the pyramid grid, and differential solution of a wave equation is realized based on a coupling second-order partial differential equation of the pyramid grid, so that the calculation efficiency can be remarkably improved, the memory consumption can be reduced, and the problem of three-dimensional complex structure imaging with rapid speed change can be finally solved on the premise of not influencing the imaging effect. By means of the method, the technical problem that the imaging cannot be carried out by simultaneously adopting the changed grids in the horizontal direction and the vertical direction in the prior art is solved, and the purpose of simultaneously adopting the changed grids in the horizontal direction and the vertical direction to carry out imaging is achieved.

It will be apparent to those skilled in the art that the modules or steps of the embodiments of the invention described above may be implemented by a general purpose computing device, they may be centralized on a single computing device or distributed across a network of multiple computing devices, and alternatively, they may be implemented by program code executable by a computing device, such that they may be stored in a storage device and executed by a computing device, and in some cases, the steps shown or described may be performed in an order different than that described herein, or they may be separately fabricated into individual integrated circuit modules, or multiple ones of them may be fabricated into a single integrated circuit module. Thus, embodiments of the invention are not limited to any specific combination of hardware and software.

The above description is only a preferred embodiment of the present invention, and is not intended to limit the present invention, and various modifications and changes may be made to the embodiment of the present invention by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (10)

1. A method of reverse time migration imaging, comprising:

s1: acquiring shot data;

s2: reading a depth domain velocity field corresponding to an imaging space of a shot point of the shot data, and converting the depth domain velocity field from a regular grid to a pyramid grid;

s3: reading anisotropic parameter sums corresponding to the imaging space, and converting the coordinates of the anisotropic parameter sums from a regular grid to a pyramid grid;

s4: reading a symmetrical axis inclination angle parameter theta and a symmetrical axis azimuth angle parameter phi corresponding to the imaging space, and converting the coordinates of the symmetrical axis inclination angle parameter theta and the symmetrical axis azimuth angle parameter phi from a regular grid to a pyramid grid;

s5: exciting the shot point by placing a wavelet, and performing finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the shot point;

s6: performing finite difference wave field simulation on the shot data through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the detection point;

s7: applying a cross-correlation imaging condition to image a wave field corresponding to the shot point and a wave field corresponding to the wave detection point;

s8: interpolating the imaging result to a regular grid, and taking the imaging result interpolated to the regular grid as the single shot imaging result of the shot point;

repeating the steps S1 to S8 on a plurality of guns to obtain the imaging results of the single guns, and superposing the imaging results of the single guns to obtain the reverse time migration imaging section.

2. The method of claim 1, wherein the coordinates are transformed from a regular grid to a pyramid grid according to the following formula:

$\left\{\begin{array}{c}{x}_p=\frac{x}{1+\mathrm{\α}z}\\ y_p=\frac{y}{1+\mathrm{\α}z}\\ z_p=\frac{1}{\mathrm{\α}}ln(1+\mathrm{\α}z)\end{array}\right.$

wherein (x, y, z) represents the coordinate system of a regular grid, (x)_p,y_p,z_p) A coordinate system representing a pyramid grid, α coordinate transform coefficients.

3. The method of claim 2, wherein the coupled second order partial differential equations of the pyramid grid are:

$\left\{\begin{array}{c}\frac{1}{v_{pz}^2}\frac{{\∂}^{2}p}{\∂t^2}=(1+2\mathrm{\ϵ})H_{2}p+H_{1}q+\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_1(p-q)\\ \frac{1}{v_{pz}^2}\frac{{\∂}^{2}q}{\∂t^2}=(1+2\mathrm{\δ})H_{2}p+H_{1}q-\frac{\mathrm{\ϵ}-\mathrm{\δ}}{\mathrm{\σ}}{H}_2(p-q)\end{array}\right.$

wherein,

$\left\{\begin{array}{c}{H}_1=\frac{e^{-2{\mathrm{\αz}}_p}}{1+{\left(xdip\right)}^2+{\left(ydip\right)}^2}\left(\begin{array}{c}{\left(xdip-{\mathrm{\αx}}_p\right)}^2\frac{{\∂}^2}{\∂x_p^2}+{\left(ydip-{\mathrm{\αy}}_p\right)}^2\frac{{\∂}^2}{\∂y_p^2}\\ +\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +2\left(xdip-{\mathrm{\αx}}_p\right)\left(ydip-{\mathrm{\αy}}_p\right)\frac{{\∂}^2}{\∂x_p\∂y_p}\\ -2\left(xdip-{\mathrm{\αx}}_p\right)\left(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p}\right)\\ -2\left(ydip-{\mathrm{\αy}}_p\right)\left(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p}\right)\end{array}\right)\\ H_2=\left(1+{\left({\mathrm{\αx}}_p\right)}^2\right)\frac{{\∂}^2}{\∂x_p^2}+\left(1+{\left({\mathrm{\αy}}_p\right)}^2\right)\frac{{\∂}^2}{\∂y_p^2}+\frac{{\∂}^2}{\∂z_p^2}-\mathrm{\α}\frac{\∂}{\∂z_p}\\ +2{\mathrm{\αx}}_p\left(\mathrm{\α}\frac{\∂}{\∂x_p}-\frac{{\∂}^2}{\∂x_p\∂z_p}\right)+2{\mathrm{\αy}}_p\left(\mathrm{\α}\frac{\∂}{\∂y_p}-\frac{{\∂}^2}{\∂y_p\∂z_p}\right)\\ +2{\mathrm{\α}}^2{x}_p{y}_p\frac{{\∂}^2}{\∂x_p\∂y_p}-H_1\end{array}\right.$

wherein xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ),p denotes the P wave wavefield, q denotes the q wave wavefield, v_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

4. The method of any of claims 1 to 3, wherein:

the method comprises the following steps of performing finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to a shot point: performing wave field compression on the wave field corresponding to the shot point according to a preset time interval, and storing a compressed result in a local disk;

the imaging of the wave field corresponding to the shot point and the wave field corresponding to the wave detection point by applying the cross-correlation imaging condition comprises the following steps: decompressing and reading the wave field corresponding to the shot point from the local disk; and imaging the wave field corresponding to the wave detection point and the wave field corresponding to the shot point which is decompressed by applying a cross-correlation imaging condition.

5. The method of any of claims 1 to 3, wherein interpolating imaging results to a regular grid comprises:

interpolating the imaging results into the regular grid in the X-direction, Y-direction and Z-direction.

6. A reverse time migration imaging apparatus, comprising:

the acquisition module is used for acquiring shot data;

the first transformation module is used for reading a depth domain velocity field corresponding to an imaging space of a shot point of the shot data and transforming the depth domain velocity field from a regular grid to a pyramid grid;

the second transformation module is used for reading anisotropic parameter sums corresponding to the imaging space and transforming the coordinates of the anisotropic parameter sums from a regular grid to a pyramid grid;

the third transformation module is used for reading a symmetry axis inclination angle parameter theta and a symmetry axis azimuth angle parameter phi corresponding to the imaging space and transforming the coordinates of the symmetry axis inclination angle parameter theta and the symmetry axis azimuth angle parameter phi to a pyramid grid from a regular grid;

the shot point difference simulation module is used for exciting the shot point by placing a wavelet and carrying out finite difference wave field simulation on a wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the shot point;

the wave detection point difference simulation module is used for carrying out finite difference wave field simulation on the shot data through a coupling second-order partial differential equation of the pyramid grid to obtain a wave field corresponding to the wave detection point;

the imaging module is used for imaging the wave field corresponding to the shot point and the wave field corresponding to the wave detection point by applying a cross-correlation imaging condition;

and the interpolation module is used for interpolating the imaging result to the regular grid and taking the imaging result interpolated to the regular grid as the single shot imaging result of the shot point.

7. The apparatus of claim 6, wherein the first transformation module and the second transformation module are specifically configured to transform coordinates from a regular grid to a pyramid grid according to the following formula:

$\left\{\begin{array}{c}{x}_p=\frac{x}{1+\mathrm{\α}z}\\ y_p=\frac{y}{1+\mathrm{\α}z}\\ z_p=\frac{1}{\mathrm{\α}}ln(1+\mathrm{\α}z)\end{array}\right.$

wherein (x, y, z) represents the coordinate system of a regular grid, (x)_p,y_p,z_p) A coordinate system representing a pyramid grid, α coordinate transform coefficients.

8. The apparatus of claim 7, wherein the coupled second order partial differential equations of the pyramid grid are:

wherein,

wherein xdip (tan (θ) × cos (Φ)), ydip (tan (θ) × sin (Φ),p denotes the P wave wavefield, q denotes the q wave wavefield, v_pzRepresenting the velocity, v, of the P-wave_szRepresenting the velocity of the S-wave.

9. The apparatus of any one of claims 6 to 8, wherein:

the shot point difference simulation module is used for carrying out wave field compression on the wave field corresponding to the shot point according to a preset time interval in the process of carrying out finite difference wave field simulation on the wave field generated by excitation through a coupling second-order partial differential equation of the pyramid grid to obtain the wave field corresponding to the shot point, and storing the compressed result in a local disk;

the imaging module is specifically used for decompressing and reading the wave field corresponding to the shot point from the local disk; and imaging the wave field corresponding to the wave detection point and the wave field corresponding to the shot point which is decompressed by applying a cross-correlation imaging condition.

10. The apparatus of any of claims 6 to 8, wherein the interpolation module is specifically configured to interpolate the imaging results into the regular grid in an X-direction, a Y-direction, and a Z-direction.