阅读说明：本技术 地震波模拟方法及系统 (Seismic wave simulation method and system ) 是由 刘炯 刘喜武 刘振峰 张金强 钱恪然 刘志远 霍志周 刘宇巍 张远银 于 2018-06-27 设计创作，主要内容包括：公开了一种地震波模拟方法及系统。该方法可以包括：步骤1：建立地震波波动方程,以g表示其中任意一个变量；步骤2：对三维空间进行网格剖分；步骤3：x方向采用傅里叶方程求解下一时刻点的内部网格点的变量g,z方向采用采用切比雪夫方程求解下一时刻点的内部网格点的变量g；步骤4：针对某一时刻点进行边界方程求解,获得下一时刻点的边界上的变量g；步骤5：随时间步进,重复步骤3-4,直至获得每个时刻点的内部网格点和边界上的变量g。本发明通过傅里叶伪谱法计算水平导数,用切比雪夫伪谱法计算垂直方向导数,在解决自由面模拟的同时,保证了模拟的高精度。(A method and system for simulating seismic waves are disclosed. The method can comprise the following steps: step 1: establishing a seismic wave equation, and expressing any variable in g; step 2: mesh generation is carried out on the three-dimensional space; and step 3: solving a variable g of an internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction; and 4, step 4: solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point; and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained. According to the method, the horizontal derivative is calculated by a Fourier pseudo-spectrum method, the vertical derivative is calculated by a Chebyshev pseudo-spectrum method, and high simulation precision is ensured while free surface simulation is solved.)

1. A seismic wave simulation method, comprising:

step 1: establishing a seismic wave equation in a stratum medium containing a free surface, and expressing any variable in the equation by g;

step 2: mesh subdivision is carried out on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto points are adopted in the z direction to obtain space discrete points;

and step 3: aiming at the internal grid point of a certain time point, solving a variable g of the internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction;

and 4, step 4: applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point;

and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained.

2. The seismic wave simulation method of claim 1, wherein solving a first order differential of a variable g in the x-direction using a fourier equation comprises:

where k denotes a wave number, Δ k is a wave number interval, Δ k is 2 pi/(N Δ x), i is an imaginary unit, G (l Δ k) represents a fourier transform of a function G (N Δ x), and G (N Δ x) is a discrete expression of G in the x direction.

3. A seismic wave simulation method according to claim 2, wherein g (n Δ x) is expressed as:

in the formula, N and l are node indexes, and N is the highest order of the Fourier polynomial.

4. The seismic simulation method of claim 1, wherein the z-direction gaussian-chebyshev-lowbarton discrete points are expressed by:

wherein the calculation interval z is limited to [ -1,1 [ ]]Between the end point z₀＝1，z_MM is the order of the chebyshev polynomial taken in the operation-1.

5. The seismic wave simulation method of claim 4, wherein solving the first order differential of the variable g in the z direction using the chebyshev equation comprises:

wherein the content of the first and second substances,

Wherein D is_cIs a matrix of (M +1) × (M +1),

6. a seismic wave simulation system having a computer program stored thereon, the program when executed by a processor implementing the steps of:

step 1: establishing a seismic wave equation in a stratum medium containing a free surface, and expressing any variable in the equation by g;

step 2: mesh subdivision is carried out on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto points are adopted in the z direction to obtain space discrete points;

and step 3: aiming at the internal grid point of a certain time point, solving a variable g of the internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction;

and 4, step 4: applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point;

and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained.

7. The seismic wave simulation system of claim 6, wherein solving a first order differential of a variable g in the x-direction using a Fourier equation comprises:

where k denotes a wave number, Δ k is a wave number interval, Δ k is 2 pi/(N Δ x), i is an imaginary unit, G (l Δ k) represents a fourier transform of a function G (N Δ x), and G (N Δ x) is a discrete expression of G in the x direction.

8. The seismic wave simulation system of claim 7, wherein the expression of g (n Δ x) is:

in the formula, N and l are node indexes, and N is the highest order of the Fourier polynomial.

9. The seismic simulation system of claim 6, wherein the z-direction gaussian-chebyshev-lowbarton discrete points are expressed as:

wherein the calculation interval z is limited to [ -1,1 [ ]]Between the end point z₀＝1，z_MM is the order of the chebyshev polynomial taken in the operation-1.

10. The seismic wave simulation system of claim 9, wherein solving the first derivative of the variable g in the z-direction using the chebyshev equation comprises:

wherein the content of the first and second substances,is a discrete expression of

Wherein D is_cIs a matrix of (M +1) × (M +1),

Technical Field

The invention relates to the technical field of oil-gas geophysical, in particular to a seismic wave simulation method and system.

Background

The pseudo-spectrum method, the finite difference method and the finite element method are combined as three main methods for seismic wave numerical simulation. The method converts the differential of physical variable to space into algebraic operation in conversion domain by mathematical conversion method, and then converts the result into physical space by inverse conversion, thus obtaining the space differential of corresponding quantity. Compared with the finite difference method and the finite element method, the pseudo-spectrum method has the characteristics of high precision and high operation speed. Therefore, the pseudo-spectrum method has wider application in seismic exploration.

At present, the pseudo spectrum method in the seismic field mainly refers to a Fourier pseudo spectrum method. The basis function in the fourier pseudo-spectrum method is a trigonometric function and has a periodic characteristic, and therefore, the spatial derivative value obtained by the fourier pseudo-spectrum method is also periodic, and the pseudo-spectrum method is generally used in the case of a periodic boundary condition and is not suitable for the case of a fixed condition. The earth medium is bounded, e.g. by the earth surface, by a fixed subsurface, and therefore the boundary conditions are often considered in the wavefield simulation. Such boundary-containing problems cannot be handled well with conventional fourier pseudospectroscopy. Therefore, there is a need to develop a seismic wave simulation method and system.

The information disclosed in this background section is only for enhancement of understanding of the general background of the invention and should not be taken as an acknowledgement or any form of suggestion that this information forms the prior art already known to a person skilled in the art.

Disclosure of Invention

The invention provides a seismic wave simulation method and system, which can calculate a horizontal derivative by a Fourier pseudo-spectrum method and a vertical derivative by a Chebyshev pseudo-spectrum method, and ensure high simulation precision while solving the problem of free surface simulation.

According to an aspect of the present invention, a seismic wave simulation method is provided. The method may include: step 1: establishing a seismic wave equation in a stratum medium containing a free surface, and expressing any variable in the equation by g; step 2: mesh subdivision is carried out on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto points are adopted in the z direction to obtain space discrete points; and step 3: aiming at the internal grid point of a certain time point, solving a variable g of the internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction; and 4, step 4: applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point; and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained.

Preferably, solving the first derivative of the variable g in the x-direction using a fourier equation comprises:

where k denotes a wave number, Δ k is a wave number interval, Δ k is 2 pi/(N Δ x), i is an imaginary unit, G (l Δ k) represents a fourier transform of a function G (N Δ x), and G (N Δ x) is a discrete expression of G in the x direction.

Preferably, the expression of g (n Δ x) is:

in the formula, N and l are node indexes, and N is the highest order of the Fourier polynomial.

Preferably, the expression for the z-direction gaussian-chebyshev-lobott discrete points is:

wherein the calculation interval z is limited to [ -1,1 [ ]]Between the end point z₀＝1，z_MM is the order of the chebyshev polynomial taken in the operation-1.

Preferably, solving the first derivative of the variable g in the z direction using the chebyshev equation comprises:

wherein the content of the first and second substances,

is a discrete expression of

Wherein D is_cIs a matrix of (M +1) × (M +1),

according to another aspect of the invention, a seismic wave simulation system is proposed, on which a computer program is stored, characterized in that the program, when executed by a processor, carries out the steps of: step 1: establishing a seismic wave equation in a stratum medium containing a free surface, and expressing any variable in the equation by g; step 2: mesh subdivision is carried out on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto points are adopted in the z direction to obtain space discrete points; and step 3: aiming at the internal grid point of a certain time point, solving a variable g of the internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction; and 4, step 4: applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point; and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained.

Preferably, solving the first derivative of the variable g in the x-direction using a fourier equation comprises:

where k denotes a wave number, Δ k is a wave number interval, Δ k is 2 pi/(N Δ x), i is an imaginary unit, G (l Δ k) represents a fourier transform of a function G (N Δ x), and G (N Δ x) is a discrete expression of G in the x direction.

Preferably, the expression of g (n Δ x) is:

in the formula, N and l are node indexes, and N is the highest order of the Fourier polynomial.

Preferably, the expression for the z-direction gaussian-chebyshev-lobott discrete points is:

wherein the calculation interval z is limited to [ -1,1 [ ]]Between the end point z₀＝1，z_MM is the order of the chebyshev polynomial taken in the operation-1.

Preferably, solving the first derivative of the variable g in the z direction using the chebyshev equation comprises:

wherein the content of the first and second substances,

is a discrete expression of

Wherein D is_cIs a matrix of (M +1) × (M +1),

the present invention has other features and advantages which will be apparent from or are set forth in detail in the accompanying drawings and the following detailed description, which are incorporated herein, and which together serve to explain certain principles of the invention.

Drawings

The above and other objects, features and advantages of the present invention will become more apparent by describing in more detail exemplary embodiments thereof with reference to the attached drawings, in which like reference numerals generally represent like parts.

FIG. 1 shows a flow chart of the steps of a seismic wave simulation method according to the invention.

Fig. 2 shows a schematic diagram of mesh generation according to an embodiment of the invention.

FIG. 3 shows a schematic diagram of a two-dimensional x-z space model according to one embodiment of the invention.

FIG. 4 shows a schematic of a wavefield snapshot of a source wavelet at a time of 0.35 seconds, in accordance with one embodiment of the present invention.

Fig. 5a, 5b, 5c show graphs comparing displacement recordings and theoretical results of detectors R1, R2, R3, respectively, obtained by simulations according to an embodiment of the present invention.

Detailed Description

The invention will be described in more detail below with reference to the accompanying drawings. While the preferred embodiments of the present invention are shown in the drawings, it should be understood that the present invention may be embodied in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.

FIG. 1 shows a flow chart of the steps of a seismic wave simulation method according to the invention.

In this embodiment, the seismic wave simulation method according to the present invention may include: step 1: establishing a seismic wave equation in a stratum medium containing a free surface, and expressing any variable in the equation by g; step 2: mesh subdivision is carried out on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto points are adopted in the z direction to obtain space discrete points; and step 3: aiming at the internal grid point of a certain time point, solving a variable g of the internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction; and 4, step 4: applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point; and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained.

In one example, solving the first derivative of the variable g in the x-direction using a fourier equation includes:

where k denotes a wave number, Δ k is a wave number interval, Δ k is 2 pi/(N Δ x), i is an imaginary unit, G (l Δ k) represents a fourier transform of a function G (N Δ x), and G (N Δ x) is a discrete expression of G in the x direction.

In one example, the expression for g (n Δ x) is:

in the formula, N and l are node indexes, and N is the highest order of the Fourier polynomial.

In one example, the expression for the gaussian-chebyshev-lowbarton discrete points in the z-direction is:

wherein the calculation interval z is limited to [ -1,1 [ ]]Between the end point z₀＝1，z_MM is the order of the chebyshev polynomial taken in the operation-1.

In one example, solving the first differential of the variable g in the z-direction using the chebyshev equation includes:

wherein the content of the first and second substances,

is a discrete expression of

Wherein D is_cIs a matrix of (M +1) × (M +1),

specifically, taking a first-order SH wave equation as an example, the seismic wave simulation method according to the present invention may include:

and (3) carrying out mesh subdivision on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto point dispersion is adopted in the z direction. According to the principle of elasticity mechanics, a first-order SH wave equation is established for a two-dimensional x-z space and is shown as a formula (6):

where V represents the velocity of the medium in the y-direction perpendicular to the x-z plane, τ_xy、τ_yzIs a stress component, f_yIs the external force per unit mass in the y-direction, ρ represents the density of the medium and μ is the shear modulus of the medium.

And aiming at a certain time point, carrying out first-order SH wave equation solution on the internal grid point of the two-dimensional x-z space to obtain the speed and the stress of the internal grid point of the next time point, wherein the X direction is solved by adopting a Fourier equation, and the Z direction is solved by adopting a Chebyshev equation.

For the x direction, equidistant grids are adopted to obtain corresponding function discrete value tau_xyV, then transforming the discretized function into the wavenumber domain as equations (7), (8), respectively:

wherein the content of the first and second substances,

representing a function τ_xyAnd Fourier transform of V. According to the relation between the Fourier transform of the function differential and the Fourier transform of the functionAnd il Δ k, and inverse fourier transforming the result back to the spatial domain to obtain a 1 st order spatial differential result of the function, i.e. the first differential of the stress component in the x direction is solved by equation (9):

wherein, tau_xyIs the stress component, k represents the wavenumber, the wavenumber spacing Δ k is 2 pi/(N Δ x), i is the imaginary unit, and the first differential of velocity is solved by equation (10):

where V denotes the velocity of the medium in the y direction perpendicular to the x-z plane, k denotes the wavenumber, and the wavenumber spacing Δ k is 2 pi/(N Δ x).

For the z-direction, the z-direction takes the Gaussian-Chebyshev-Lobarton point dispersion, expressed in the form of equation (3), where the calculation interval z is limited to [ -1,1]Between the end point z₀＝1，z_M＝-1，[a,b]Is a physical interval in the z direction (a and b can be any bounded real number) and is converted by coordinatesMapping to a calculated coordinate z ', where z' belongs to [ -1,1 [ ]]Then, the first differential of the stress component in the z direction is solved by equation (11):

wherein

Is a discrete expression of

Wherein the content of the first and second substances,

in formula (11)

[a,b]Is the physical interval in the z direction.

The first differential of the velocity in the z direction is solved by equation (13):

wherein

Is a discrete expression of

Wherein the content of the first and second substances,

in formula (14)[a,b]Is the physical interval in the z direction.

To obtainAnd

then, the value of the right term of equation (6) can be obtained as (f)_yρ and μ are known), then equation (6) is changed from partial differential equation set to ordinary differential equation set, and then numerical integration technique (such as eulerian method, 4-step Rungestota method, etc.) is used to find the next timeV, τ of_xyAnd τ_yz。

Applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation (15) of the boundary of the two-dimensional x-z space aiming at a certain time point:

τ _yz0, when z is 0 (15)

The displacement and stress on the boundary at the next point in time are obtained.

And (4) obtaining the speed, stress and displacement and stress on the boundary of the inner grid point of each time point along with time stepping. And obtaining the wave field record of the SH wave according to the speed and the stress of the internal grid point of each time point and the displacement and the stress on the boundary, and obtaining the change relation of the displacement of each space field along with time.

According to the method, the horizontal derivative is calculated by a Fourier pseudo-spectrum method, the vertical derivative is calculated by a Chebyshev pseudo-spectrum method, and high simulation precision is guaranteed while free surface simulation is achieved.

Application example

To facilitate understanding of the solution of the embodiments of the present invention and the effects thereof, a specific application example is given below. It will be understood by those skilled in the art that this example is merely for the purpose of facilitating an understanding of the present invention and that any specific details thereof are not intended to limit the invention in any way.

Fig. 2 shows a schematic diagram of mesh generation according to an embodiment of the invention.

FIG. 3 shows a schematic diagram of a two-dimensional x-z space model according to one embodiment of the invention.

FIG. 4 shows a schematic of a wavefield snapshot of a source wavelet at a time of 0.35 seconds, in accordance with one embodiment of the present invention.

The seismic wave simulation method comprises the following steps:

and (2) carrying out mesh subdivision on a three-dimensional space, wherein equidistant meshes are adopted in the x direction, Gaussian-Chebyshev-Lobarton point dispersion is adopted in the z direction, as shown in figure 2, SH wave number value simulation is carried out on the two-dimensional x-z space shown in figure 2, 3 virtual detectors R1, R2 and R3 are arranged to record the displacement of a space point, as shown in figure 3, and a seismic source in the simulation adopts a 50Hz Ricker wavelet, as shown in figure 4.

According to the principle of elasticity mechanics, a first-order SH wave equation is established for a two-dimensional x-z space as a formula (6).

And aiming at a certain time point, carrying out first-order SH wave equation solution on the internal grid point of the two-dimensional x-z space to obtain the speed and the stress of the internal grid point of the next time point, wherein the X direction is solved by adopting a Fourier equation, and the Z direction is solved by adopting a Chebyshev equation.

For the x direction, equidistant grids are adopted to obtain corresponding function discrete value tau_xyV, and then transforming the discretized function into a wavenumber domain as formulas (7) and (8). According to the relation between the Fourier transform of the function differential and the Fourier transform of the function

And il deltak, and inverse fourier transforming the result back to the spatial domain to obtain a 1 st order spatial differential result of the function, i.e. solving the first order differential of the stress component in the x direction by equation (9) and the first order differential of the velocity by equation (10).

For the z-direction, the z-direction employs Gaussian-Chebyshev-Lobarton point dispersion, where the computation interval z is limited to [ -1,1]Between the end point z₀＝1，z_M＝-1，[a,b]Is a physical interval in the z direction (a and b can be any bounded real number) and is converted by coordinates

Mapping to a calculated coordinate z ', where z' belongs to [ -1,1 [ ]]Then, the first order differential of the stress component in the z direction is solved by the formula (11), and the first order differential of the velocity in the z direction is solved by the formula (13), wherein the formulas (12) and (14) represent

A discrete expression of (c).

To obtain

And

then, the value of the right term of equation (6) can be obtained as (f)_yρ and μ are known), then equation (6) is changed from partial differential equation set to ordinary differential equation set, and then the numerical integration technique (e.g. eulerian, 4-step Rungestota, etc.) is used to find V, τ at the next time_xyAnd τ_yz。

And applying a physical boundary condition and a numerical boundary condition, and solving a boundary equation (15) of the boundary of the two-dimensional x-z space aiming at a certain time point to obtain the displacement and the stress on the boundary of the next time point.

And (4) obtaining the speed, stress and displacement and stress on the boundary of the inner grid point of each time point along with time stepping. And obtaining the wave field record of the SH wave according to the speed and the stress of the internal grid point of each time point and the displacement and the stress on the boundary, and obtaining the change relation of the displacement of each space field along with time.

5a, 5b and 5c show graphs comparing displacement records of detectors R1, R2 and R3 obtained by simulation according to an embodiment of the invention and theoretical results, and it can be seen that the numerical solution and the theoretical solution are well matched, so that the method is proved to have high simulation precision.

In conclusion, the horizontal derivative is calculated by the Fourier pseudo-spectrum method, the vertical derivative is calculated by the Chebyshev pseudo-spectrum method, and the high precision of the simulation is ensured while the free surface simulation is solved.

It will be appreciated by persons skilled in the art that the above description of embodiments of the invention is intended only to illustrate the benefits of embodiments of the invention and is not intended to limit embodiments of the invention to any examples given.

A seismic wave simulation system according to the invention, having stored thereon a computer program which, when executed by a processor, performs the steps of: step 1: establishing a seismic wave equation in a stratum medium containing a free surface, and expressing any variable in the equation by g; step 2: mesh subdivision is carried out on the three-dimensional space, wherein equidistant meshes are adopted in the x direction, and Gaussian-Chebyshev-Luobatto points are adopted in the z direction to obtain space discrete points; and step 3: aiming at the internal grid point of a certain time point, solving a variable g of the internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction; and 4, step 4: applying physical boundary conditions and numerical boundary conditions, and solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point; and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained.

In one example, solving the first derivative of the variable g in the x-direction using a fourier equation includes:

where k denotes a wave number, Δ k is a wave number interval, Δ k is 2 pi/(N Δ x), i is an imaginary unit, G (l Δ k) represents a fourier transform of a function G (N Δ x), and G (N Δ x) is a discrete expression of G in the x direction.

In one example, the expression for g (n Δ x) is:

in the formula, N and l are node indexes, and N is the highest order of the Fourier polynomial.

In one example, the expression for the gaussian-chebyshev-lowbarton discrete points in the z-direction is:

wherein the calculation interval z is limited to [ -1,1 [ ]]Between the end point z₀＝1，z_MM is the order of the chebyshev polynomial taken in the operation-1.

In one example, solving the first differential of the variable g in the z-direction using the chebyshev equation includes:

wherein the content of the first and second substances,

is a discrete expression of

Wherein D is_cIs a matrix of (M +1) × (M +1),

the system calculates the horizontal derivative by a Fourier pseudo-spectrum method and calculates the vertical derivative by a Chebyshev pseudo-spectrum method, thereby solving the problem of free surface simulation and ensuring the high precision of the simulation.

Having described embodiments of the present invention, the foregoing description is intended to be exemplary, not exhaustive, and not limited to the embodiments disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments.