阅读说明：本技术 弹性波场数值模拟方法及系统 (Elastic wave field numerical simulation method and system ) 是由 李文杰 刘来详 陈世军 于 2018-06-15 设计创作，主要内容包括：公开了一种弹性波场数值模拟方法及系统。该方法可以包括：步骤1：建立速度模型,确定空间和时间的采样间隔,进行有限差分网格离散化；步骤2：根据弹性波动方程计算内部波场数据；步骤3：根据角点波场数值计算公式计算角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算边界波场数据；步骤5：获得速度模型的波场数据,判断波场边界的反射情况是否满足预定标准,若是,执行步骤6,若否,重新确定波场吸收系数,重复步骤4-5；步骤6：输出最终波场数据。本发明的模型波场值能准确反映出弹性波在二维各向同性介质中传播的真实性,达到提高弹性波数值模拟结果精度的目的。(An elastic wave field numerical simulation method and system are disclosed. The method can comprise the following steps: step 1: establishing a speed model, determining sampling intervals of space and time, and performing finite difference grid discretization; step 2: calculating internal wave field data according to an elastic wave equation; and step 3: calculating angular point wave field data according to an angular point wave field numerical calculation formula; and 4, step 4: calculating boundary wave field data according to a wave field absorption boundary equation and a wave field absorption coefficient; and 5: acquiring wave field data of the velocity model, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing the step 6, otherwise, re-determining the wave field absorption coefficient, and repeating the steps 4-5; step 6: and outputting final wave field data. The model wave field value of the invention can accurately reflect the reality of the elastic wave propagating in the two-dimensional isotropic medium, thereby achieving the purpose of improving the accuracy of the numerical simulation result of the elastic wave.)

1. An elastic wave field numerical simulation method comprising:

step 1: establishing a speed model, determining a space sampling interval and a time sampling interval of the speed model, and performing finite difference grid discretization on the speed model;

step 2: calculating wave field values of all moments of all grid points in the velocity model according to an elastic wave equation to obtain internal wave field data;

and step 3: calculating wave field values of all moments corresponding to four angular points of the velocity model according to an angular point wave field numerical calculation formula to obtain angular point wave field data;

and 4, step 4: calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to a wave field absorption boundary equation and a wave field absorption coefficient to obtain boundary wave field data;

and 5: obtaining wave field data of the velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing a step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5;

step 6: outputting final wavefield data of the velocity model.

2. The elastic wave field numerical simulation method of claim 1, wherein the expression of the wave field value at each time for each grid point inside the velocity model is:

wherein the content of the first and second substances,

u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

α, β are the longitudinal and transverse wave velocities, respectively, t is the time, and X and Z are the coordinates along the X and Z axes in the case of a two-dimensional medium.

3. The elastic wavefield numerical simulation method of claim 1, wherein the wavefield absorption boundary equation is obtained by:

performing Fourier transform on the elastic wave equation to obtain an elastic wave equation after Fourier transform;

factorizing and simplifying the elastic wave equation after Fourier transformation to obtain a simplified equation;

and converting the simplified equation into a time domain to obtain a wave field absorption boundary equation.

4. The elastic wave field numerical simulation method of claim 3, wherein the elastic wave equation after Fourier transform is:

5. the elastic wavefield numerical simulation method of claim 3, wherein the simplified equation is:

6. the elastic wave field numerical simulation method according to claim 3, wherein the wave field values at respective time instants corresponding to grid points on the boundary of the velocity model are calculated by equations (4a) - (4d),

the absorption boundary equation corresponding to the bottom edge of the velocity model is as follows:

the absorption boundary equation corresponding to the top edge of the velocity model is:

the absorption boundary equation corresponding to the left side of the velocity model is:

the absorption boundary equation corresponding to the right side of the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

r₁，r₂、r₃，r₄is wave field absorption coefficient, 1 is more than or equal to r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄More than or equal to 0, and alpha and beta are longitudinal wave speed and transverse wave speed respectively.

7. The elastic wave field numerical simulation method of claim 1, wherein the wave field values at each time corresponding to four corner points of the velocity model are calculated by equations (5a) - (5 d):

the wavefield calculation equation for the lower right corner point is:

the wave field calculation equation of the corner point at the upper left corner is as follows:

the wave field calculation equation of the corner point at the upper right corner is as follows:

the wavefield calculation equation for the lower left corner point is:

wherein:u and w are respectively a horizontal displacement field and a vertical displacement field corresponding to the elastic wave field;

alpha and beta are respectively the longitudinal wave speed and the transverse wave speed, and X and Z are coordinates along an X axis and a Z axis under the condition of a two-dimensional medium;

wherein:

θ is the angle of incidence of the wave.

8. An elastic wave field numerical simulation system having a computer program stored thereon, wherein the program when executed by a processor performs the steps of:

step 1: establishing a speed model, determining a space sampling interval and a time sampling interval of the speed model, and performing finite difference grid discretization on the speed model;

step 2: calculating wave field values of all moments of all grid points in the velocity model according to an elastic wave equation to obtain internal wave field data;

and step 3: calculating wave field values of all moments corresponding to four angular points of the velocity model according to an angular point wave field numerical calculation formula to obtain angular point wave field data;

and 4, step 4: calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to a wave field absorption boundary equation and a wave field absorption coefficient to obtain boundary wave field data;

and 5: obtaining wave field data of the velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing a step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5;

step 6: outputting final wavefield data of the velocity model.

9. The elastic wave field numerical simulation system of claim 8, wherein obtaining a wave field absorption boundary equation from the elastic wave equation, calculating wave field values at respective time instants corresponding to grid points on a boundary of the velocity model from the wave field absorption boundary equation, and obtaining boundary wave field data comprises:

performing Fourier transform on the elastic wave equation to obtain an elastic wave equation after Fourier transform;

factorizing and simplifying the elastic wave equation after Fourier transformation to obtain a simplified equation;

converting the simplified equation into a time domain to obtain a wave field absorption boundary equation;

and calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to the wave field absorption boundary equation to obtain boundary wave field data.

10. The elastic wave field numerical simulation system of claim 9, wherein the wave field values at respective time instants corresponding to grid points on the boundary of the velocity model are calculated by equations (4a) - (4d),

the absorption boundary equation corresponding to the bottom edge of the velocity model is as follows:

the absorption boundary equation corresponding to the top edge of the velocity model is:

the absorption boundary equation corresponding to the left side of the velocity model is:

the absorption boundary equation corresponding to the right side of the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

r₁，r₂、r₃，r₄is wave field absorption coefficient, 1 is more than or equal to r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄More than or equal to 0, and alpha and beta are longitudinal wave speed and transverse wave speed respectively.

Technical Field

The invention relates to the technical field of elastic wave field numerical simulation, in particular to an elastic wave field numerical simulation method and system.

Background

The elastic wave field numerical simulation technology is an important technology in seismic exploration, and by utilizing a forward modeling technology in the elastic wave field numerical simulation technology, whether underground physical parameters finally obtained through seismic data processing and geological interpretation and an inverted speed model can truly reflect an underground oil storage structure can be checked, if the elastic wave field numerical simulation technology is left, the processing and interpretation work of a plurality of seismic data at present is difficult, the elastic wave field numerical simulation technology widely applied to actual work at present is mainly an elastic wave field finite difference wave field numerical simulation technology, in the elastic wave finite difference wave field numerical simulation, due to the limitation of a calculation model, an inevitable problem encountered by people is how to solve artificial boundary reflection energy caused by calculating grid boundaries, and the artificial reflection from the boundaries distorts the authenticity of waves propagating in an infinite medium, therefore, it is necessary to provide an algorithm for propagating waves on the boundary to enable the energy of the outward propagating waves to penetrate through the boundary, otherwise, we have to enlarge the boundary of the model, and the calculation workload is greatly increased. In order to solve the boundary reflection problem, a great deal of research has been carried out for decades, and many different methods for calculating boundary wavefields have been proposed, but until now, there is no perfect solution, and the boundary problem in the elastic wave finite difference wavefield numerical simulation is more complicated and more difficult to solve than the acoustic wave field numerical simulation, and at present, the boundary problem research in the elastic wave finite difference wavefield numerical simulation is still an important subject in the field, and algorithms for different boundary problems are also successively produced, wherein the method for eliminating the boundary reflection by the wavefield absorption is more general and has more obvious effect, and many research results and documents in this respect appear in domestic and foreign countries, and at present, two types of absorption boundary conditions are widely used in exploration: sponge absorption boundary conditions and paraxial approximation absorption boundary conditions. The sponge absorption boundary condition attenuates incident waves in a band range of a viscous boundary or a band range close to the boundary; the paraxial approximate absorption boundary condition is based on the principle of a one-way wave propagation method, and one-way wave equations under different approximate conditions are used as the absorption boundary conditions at different boundaries. In China, Clayton populates the concept of borderless reflection in one-dimensional situations proposed by predecessors into numerical simulation of three-dimensional anisotropic medium elastic waves using a feature analysis method in literature (Doryngo. absorption boundary conditions in elastic wave numerical simulation, 1999,34(1): 45-56) to obtain absorption boundary conditions in TI media, and in abroad, Clayton and Engquist et al in literature (Robert Clayton, B jorn Engquist. absorbing boundary conditions for both acidic and elastic wave propagation. Shell. Seis. Soeq. Am.,1977,67(6): 1529-1540) obtain absorption boundary conditions for elastic wave numerical simulation in isotropic media using the elastic wave propagation equation axis approximation theory, and so far, the absorption boundary conditions proposed by predecessors still perform elastic wave differential numerical simulation in isotropic media to solve the problem of borderline difference and the problem in Clayton (later), bjorn Engquist, absorbing boundary conditions for wave equalizations, geophysics,1980,45(3): 895-904) have also made more intensive studies on the absorption boundary conditions in the elastic wave finite difference wave field numerical simulation. In China, a plurality of researchers do more detailed research on the absorption boundary conditions under the condition of performing the numerical simulation of the elastic wave field by using the finite difference method to obtain a plurality of valuable research results, and the algorithm of the researchers obtains some better results on certain models. Through existing elastic wave field numerical simulation algorithms, each method has own characteristics, but all the absorption boundary conditions cannot completely absorb the reflected wave field energy from the boundary, and only partially absorb the reflected energy from the boundary. Therefore, for the elastic wave number value simulation, it is really needed to invent an absorption boundary condition that can completely absorb the reflected energy from the boundary, so as to improve the accuracy of the elastic wave numerical simulation result and make the numerical simulation result more conform to the propagation rule of the actual wave field. Therefore, there is a need to develop an elastic wave field numerical 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 an elastic wave field numerical simulation method and system, which can effectively eliminate the reflection energy generated from the limited speed model boundary by continuously adjusting the wave field absorption coefficient, and can completely absorb the reflection energy from the model boundary when the absorption coefficient is adjusted to the optimal state, so that the calculated model wave field value can accurately reflect the reality of the elastic wave propagating in a two-dimensional isotropic medium, and the aim of improving the accuracy of the elastic wave numerical simulation result is fulfilled.

According to an aspect of the invention, an elastic wave field numerical simulation method is provided. The method may include: step 1: establishing a speed model, determining a space sampling interval and a time sampling interval of the speed model, and performing finite difference grid discretization on the speed model; step 2: calculating wave field values of all moments of all grid points in the velocity model according to an elastic wave equation to obtain internal wave field data; and step 3: calculating wave field values of all moments corresponding to four angular points of the velocity model according to an angular point wave field numerical calculation formula to obtain angular point wave field data; and 4, step 4: calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to a wave field absorption boundary equation and a wave field absorption coefficient to obtain boundary wave field data; and 5: obtaining wave field data of the velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing a step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5; step 6: outputting final wavefield data of the velocity model.

Preferably, the expression of the wave field value at each time of each grid point inside the velocity model is:

wherein the content of the first and second substances,

u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

α, β are the longitudinal and transverse wave velocities, respectively, t is the time, and X and Z are the coordinates along the X and Z axes in the case of a two-dimensional medium.

Preferably, the wavefield absorption boundary equation is obtained by: performing Fourier transform on the elastic wave equation to obtain an elastic wave equation after Fourier transform; factorizing and simplifying the elastic wave equation after Fourier transformation to obtain a simplified equation; and converting the simplified equation into a time domain to obtain a wave field absorption boundary equation.

Preferably, the elastic wave equation after fourier transform is:

preferably, the simplified equation is:

preferably, wave field values at respective time instants corresponding to grid points on the boundary of the velocity model are calculated by equations (4a) - (4d),

the absorption boundary equation corresponding to the bottom edge of the velocity model is as follows:

the absorption boundary equation corresponding to the top edge of the velocity model is:

the absorption boundary equation corresponding to the left side of the velocity model is:

the absorption boundary equation corresponding to the right side of the velocity model is:

wherein the content of the first and second substances,

u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

r₁，r₂、r₃，r₄is wave field absorption coefficient, 1 is more than or equal to r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄Not less than 0, alpha and beta are longitudinal and transverse waves respectivelySpeed.

Preferably, wave field values of each time corresponding to four corner points of the velocity model are calculated through equations (5a) - (5d),

the wavefield calculation equation for the lower right corner point is:

the wave field calculation equation of the corner point at the upper left corner is as follows:

the wave field calculation equation of the corner point at the upper right corner is as follows:

the wavefield calculation equation for the lower left corner point is:

wherein:

u and w are respectively a horizontal displacement field and a vertical displacement field corresponding to the elastic wave field;

alpha and beta are respectively the longitudinal wave speed and the transverse wave speed, and X and Z are coordinates along an X axis and a Z axis under the condition of a two-dimensional medium;

wherein:

θ is the angle of incidence of the wave.

According to another aspect of the invention, an elastic wave field numerical simulation system is proposed, on which a computer program is stored, wherein the program, when being executed by a processor, realizes the following steps: step 1: establishing a speed model, determining a space sampling interval and a time sampling interval of the speed model, and performing finite difference grid discretization on the speed model; step 2: calculating wave field values of all moments of all grid points in the velocity model according to an elastic wave equation to obtain internal wave field data; and step 3: calculating wave field values of all moments corresponding to four angular points of the velocity model according to an angular point wave field numerical calculation formula to obtain angular point wave field data; and 4, step 4: calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to a wave field absorption boundary equation and a wave field absorption coefficient to obtain boundary wave field data; and 5: obtaining wave field data of the velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing a step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5; step 6: outputting final wavefield data of the velocity model.

Preferably, the expression of the wave field value at each time of each grid point inside the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

α, β are the longitudinal and transverse wave velocities, respectively, t is the time, and X and Z are the coordinates along the X and Z axes in the case of a two-dimensional medium.

Preferably, the wavefield absorption boundary equation is obtained by: performing Fourier transform on the elastic wave equation to obtain an elastic wave equation after Fourier transform; factorizing and simplifying the elastic wave equation after Fourier transformation to obtain a simplified equation; and converting the simplified equation into a time domain to obtain a wave field absorption boundary equation.

Preferably, the elastic wave equation after fourier transform is:

preferably, the simplified equation is:

preferably, wave field values at respective time instants corresponding to grid points on the boundary of the velocity model are calculated by equations (4a) - (4d),

the absorption boundary equation corresponding to the bottom edge of the velocity model is as follows:

the absorption boundary equation corresponding to the top edge of the velocity model is:

the absorption boundary equation corresponding to the left side of the velocity model is:

the absorption boundary equation corresponding to the right side of the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

r₁，r₂、r₃，r₄is wave field absorption coefficient, 1 is more than or equal to r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄More than or equal to 0, and alpha and beta are longitudinal wave speed and transverse wave speed respectively.

Preferably, wave field values of each time corresponding to four corner points of the velocity model are calculated through equations (5a) - (5d),

the wavefield calculation equation for the lower right corner point is:

the wave field calculation equation of the corner point at the upper left corner is as follows:

the wave field calculation equation of the corner point at the upper right corner is as follows:

the wavefield calculation equation for the lower left corner point is:

wherein:u and w are levels corresponding to the elastic wave field respectivelyA displacement field and a vertical displacement field;

alpha and beta are respectively the longitudinal wave speed and the transverse wave speed, and X and Z are coordinates along an X axis and a Z axis under the condition of a two-dimensional medium;

wherein:

θ is the angle of incidence of the wave.

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 an elastic wave field numerical simulation method according to the invention.

FIG. 2 shows a schematic diagram of a two-dimensional isotropic medium constant velocity model according to one embodiment of the invention.

FIGS. 3a and 3b show diagrams of the X and Z components of the wavefield propagating to 120 ms under rigid boundary conditions, respectively, in accordance with one embodiment of the present invention.

FIGS. 4a and 4b show diagrams of the X and Z components of the wavefield propagating to 120 ms under absorption boundary conditions, respectively, in accordance with one embodiment of the present invention.

FIG. 5 shows a schematic diagram of a two-layer media velocity model according to one embodiment of the invention.

Fig. 6a and 6b show schematic diagrams of horizontal displacement and vertical displacement, respectively, when propagating to 260 ms under the condition of elastic wave double absorption boundary according to an embodiment of the present invention.

FIGS. 7a and 7b show diagrams of wavefield horizontal displacement and vertical displacement, respectively, as propagated to 260 milliseconds under the Clayton absorption boundary condition, in accordance with one 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 an elastic wave field numerical simulation method according to the invention.

In this embodiment, the elastic wave field numerical simulation method according to the present invention may include: step 1: establishing a speed model, determining a space sampling interval and a time sampling interval of the speed model, and performing finite difference grid discretization on the speed model; step 2: calculating wave field values of all moments of all grid points in the velocity model according to an elastic wave equation to obtain internal wave field data; and step 3: calculating wave field values of each moment corresponding to four angular points of the velocity model according to an angular point wave field numerical calculation formula to obtain angular point wave field data; and 4, step 4: calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to the wave field absorption boundary equation and the wave field absorption coefficient to obtain boundary wave field data; and 5: obtaining wave field data of a velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing a step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5; step 6: and outputting final wave field data of the velocity model.

In one example, the expression of the wave field values at various times for various grid points within the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

alpha and beta are the longitudinal and transverse wave speeds, respectively, and t is time.

In one example, the wavefield absorption boundary equation is obtained by: performing Fourier transform on the elastic wave equation to obtain the elastic wave equation after Fourier transform; factorizing and simplifying the elastic wave equation after Fourier transformation to obtain a simplified equation; and converting the simplified equation into a time domain to obtain a wave field absorption boundary equation.

In one example, the elastic wave equation after fourier transform is:

in one example, the simplified equation is:

in one example, the wave field values at respective time instants corresponding to the grid points on the boundary of the velocity model are calculated by equations (4a) - (4d),

the absorption boundary equation corresponding to the bottom edge of the velocity model is as follows:

the absorption boundary equation corresponding to the top edge of the velocity model is:

the absorption boundary equation corresponding to the left side of the velocity model is:

the absorption boundary equation corresponding to the right side of the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

r₁，r₂、r₃，r₄is wave field absorption coefficient, 1 is more than or equal to r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄More than or equal to 0, and alpha and beta are longitudinal wave speed and transverse wave speed respectively.

In one example, wave field values at respective time instants corresponding to four corner points of the velocity model are calculated by equations (5a) - (5d),

the wavefield calculation equation for the lower right corner point is:

the wave field calculation equation of the corner point at the upper left corner is as follows:

the wave field calculation equation of the corner point at the upper right corner is as follows:

the wavefield calculation equation for the lower left corner point is:

wherein:u and w are respectively a horizontal displacement field and a vertical displacement field corresponding to the elastic wave field;

alpha and beta are respectively the longitudinal wave speed and the transverse wave speed, and X and Z are coordinates along an X axis and a Z axis under the condition of a two-dimensional medium;

wherein:

θ is the angle of incidence of the wave.

Specifically, the elastic wave field numerical simulation method according to the present invention may include:

step 1: establishing a velocity model, determining the position of a shot point (namely an excitation source), selecting seismic wavelets, determining the space sampling interval and the time sampling interval of the velocity model, and performing finite difference grid discretization on the velocity model.

Step 2: under the condition of an isotropic medium, assuming that X and Z are coordinates along an X axis and a Z axis in the case of a two-dimensional medium, the X axis is directed to the right, and the Z axis is directed downward, we can describe the motion of a P wave in the medium and the motion of a SV wave polarized vertically by using two conjugated second-order partial differential equations, where, without considering the SH wave polarized along the horizontal direction, u and w are displacements in the horizontal direction and the vertical direction, ρ is the density of the medium, t is time, and λ and μ are the ramet coefficients in the specific medium, then the elastic wave equation about the isotropic inhomogeneous medium can be obtained as formula (5):

assuming that the density ρ is constant, equation (6) can be viewed as a function of P-wave and SV-wave velocities as a function of spatial position. Wherein λ and μ are related to longitudinal and transverse wave velocities α (x, z) and β (x, z) by equation (6):

and then, the formula (6) is rewritten into a matrix form, namely the formula (1), the wave field value of each time of each grid point in the velocity model is calculated into the formula (1), and the internal wave field data is obtained.

And step 3: calculating wave field values of four corner points of the velocity model according to a corner point wave field numerical calculation formula to be respectively formulas (5a) - (5d), obtaining corner point wave field data, and calculating sin theta and cos theta according to the position coordinates of the shot point and the position coordinates corresponding to the four corner points, thereby obtaining each item value in the J matrix vector.

And 4, step 4: obtaining a wave field absorption boundary equation according to the elastic wave equation and the wave field absorption coefficient, and calculating wave field values at each moment corresponding to grid points on a boundary of a velocity model corresponding to the wave field absorption coefficient according to the wave field absorption boundary equation to obtain boundary wave field data, wherein the wave field absorption boundary equation comprises the following steps:

performing Fourier transform on the elastic wave equation to obtain an elastic wave equation after Fourier transform as a formula (2); to be able to decompose the left term of equation (2), consider changing the left term of equation (2) to equation (8):

wherein:

and (3) finishing the formula (8) to obtain a formula (9):

obtaining formula (10) from formula (9):

order:then, formula (11):

further, equation (12) is obtained:

wherein: 1 is more than or equal to r₁≥0；1≥r₂≥0。

Thus, the expression (10) can be written as the following derivation:

formula (13) may be rewritten into the following form:

I+C(k_z/ω)-A(k_x/ω)-B(k_z/ω)＝0 (14)

wherein:

thus, by factoring, equation (2) can be rewritten into the following form:

let G be I + C (k)_z/ω)+A(k_x/ω)+B(k_zω) on both sides of formula (15) by G^-1Obtaining a simplified equation as a formula (3); converting the simplified equation into a time domain to obtain a wave field absorption boundary equation; and calculating wave field values at each moment corresponding to grid points on the boundary of the velocity model corresponding to the wave field absorption coefficient according to a wave field absorption boundary equation to obtain boundary wave field data, wherein the wave field values are in formulas (4a) - (4 d).

And 5: and (3) acquiring wave field data of a velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing the step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5, wherein the preset standard can be determined by a person skilled in the art according to different velocity models.

Step 6: and outputting final wave field data of the velocity model.

The method effectively eliminates the reflection energy generated by the limited speed model boundary by continuously adjusting the wave field absorption coefficient, and can completely absorb the reflection energy from the model boundary when the absorption coefficient is adjusted to the optimal state, so that the calculated model wave field value can accurately reflect the reality of the elastic wave propagating in the two-dimensional isotropic medium, and the aim of improving the accuracy of the numerical simulation result of the elastic wave is fulfilled.

To facilitate understanding of the aspects of the embodiments of the present invention and their effects, two specific application examples are 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.

Application example 1

FIG. 2 shows a schematic diagram of a two-dimensional isotropic medium constant velocity model according to one embodiment of the invention.

Step 1: and (3) establishing a two-dimensional isotropic medium constant velocity model, as shown in FIG. 2, wherein the longitudinal wave velocity is 3000 m/s, the transverse wave velocity is 2000 m/s, and the excitation source is a longitudinal wave seismic source and is positioned in the center of the velocity model. Selecting a Rake wavelet with a main frequency of 50Hz as a seismic wavelet, wherein the space sampling interval is 1 meter, the time sampling interval is 0.25ms, and then carrying out grid discretization on the selected velocity model according to the selected space sampling.

Step 2: and (3) calculating wave field values of all moments of all grid points in the velocity model according to the formula (1) to obtain internal wave field data.

And step 3: and calculating wave field values of the angular points of the velocity model into formulas (5a) - (5d) according to the angular point wave field numerical calculation formula to obtain angular point wave field data.

And 4, step 4: and obtaining a wave field absorption boundary equation according to the elastic wave equation and the wave field absorption coefficient, calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model corresponding to the wave field absorption coefficient according to the wave field absorption boundary equation to be formulas (4a) - (4d), and obtaining boundary wave field data.

And 5: and (3) acquiring wave field data of the velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing the step 6, otherwise, re-determining the wave field absorption coefficient, and repeating the steps 4-5.

Step 6: and outputting final wave field data of the velocity model.

FIGS. 3a and 3b show diagrams of the X and Z components of a wavefield propagating to 120 milliseconds under rigid boundary conditions (i.e., total reflection occurs as the wavefield propagates to the boundary), respectively, according to one embodiment of the present invention. FIGS. 4a and 4b show diagrams of the X and Z components of the wavefield propagating to 120 ms under absorption boundary conditions, respectively, in accordance with one embodiment of the present invention. Comparing fig. 3 and 4, it can be easily seen that the absorption boundary condition used in the present invention can perform a good absorption function for the wave field propagating to the boundary.

Application example 2

FIG. 5 shows a schematic diagram of a two-layer media velocity model according to one embodiment of the invention.

Step 1: a two-dimensional isotropic laminar medium velocity model is established, as shown in fig. 5, the longitudinal wave velocity of the first layer of medium is 2000 m/s, the transverse wave velocity is 1000 m/s, the longitudinal wave velocity of the second layer of medium is 4000 m/s, the transverse wave velocity is 2000 m/s, the excitation source is located 300 m below the top of the model and 300 m away from the left side and the right side respectively. Selecting a Rake wavelet with a main frequency of 50Hz as a seismic wavelet, wherein the space sampling interval is 1 meter, the time sampling interval is 0.25ms, and then carrying out grid discretization on the selected velocity model according to the selected space sampling.

Step 2: and (3) calculating wave field values of all moments of all grid points in the velocity model according to the formula (1) to obtain internal wave field data.

And step 3: and calculating wave field values of the angular points of the velocity model according to the angular point wave field value calculation formula into formulas (5a) - (5d) respectively to obtain angular point wave field data.

And 4, step 4: and obtaining a wave field absorption boundary equation according to the elastic wave equation and the wave field absorption coefficient, calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model corresponding to the wave field absorption coefficient according to the wave field absorption boundary equation to be formulas (4a) - (4d), and obtaining boundary wave field data.

And 5: and (3) acquiring wave field data of the velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing the step 6, otherwise, re-determining the wave field absorption coefficient, and repeating the steps 4-5.

Step 6: and outputting final wave field data of the velocity model.

Fig. 6a and 6b show schematic diagrams of horizontal displacement and vertical displacement, respectively, when propagating to 260 ms under the condition of elastic wave double absorption boundary according to an embodiment of the present invention. FIGS. 7a and 7b show diagrams of wavefield horizontal displacement and vertical displacement, respectively, as propagated to 260 milliseconds under the Clayton absorption boundary condition, in accordance with one embodiment of the present invention. In the upper medium, the direct longitudinal wave is transmitted out of the model, and the reflected wave from the interface of the two layers of media also quickly reaches the top boundary; in the lower medium, the transmitted wave has also passed to the bottom, and only the transmitted wave and the reflected wave exist on the wave field distribution diagram of the model, and the wave field distribution conditions under the Clayton absorption boundary condition and the elastic wave double absorption boundary condition are compared on the diagram, so that the following can be seen: reflection from the boundary appears on a wave field distribution diagram obtained by utilizing the Clayton absorption boundary condition, but the elastic wave double absorption boundary condition provided by the invention can receive good absorption effect.

In conclusion, the invention effectively eliminates the reflection energy generated from the limited speed model boundary by continuously adjusting the wave field absorption coefficient, and when the absorption coefficient is adjusted to the optimal state, the reflection energy from the model boundary can be completely absorbed, so that the calculated model wave field value can accurately reflect the reality of the elastic wave propagating in the two-dimensional isotropic medium, and the aim of improving the accuracy of the elastic wave numerical simulation result is fulfilled.

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.

According to another aspect of the invention, an elastic wave field numerical simulation system is proposed, on which a computer program is stored, wherein the program, when executed by a processor, performs the steps of: step 1: establishing a speed model, determining a space sampling interval and a time sampling interval of the speed model, and performing finite difference grid discretization on the speed model; step 2: calculating wave field values of all moments of all grid points in the velocity model according to an elastic wave equation to obtain internal wave field data; and step 3: calculating wave field values of each moment corresponding to four angular points of the velocity model according to an angular point wave field numerical calculation formula to obtain angular point wave field data; and 4, step 4: calculating wave field values of all moments corresponding to grid points on the boundary of the velocity model according to the wave field absorption boundary equation and the wave field absorption coefficient to obtain boundary wave field data; and 5: obtaining wave field data of a velocity model according to the boundary wave field data, the internal wave field data and the angular point wave field data, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing a step 6, otherwise, re-determining a wave field absorption coefficient, and repeating the steps 4-5; step 6: and outputting final wave field data of the velocity model.

In one example, the expression of the wave field values at various times for various grid points within the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field;

alpha and beta are the longitudinal and transverse wave speeds, respectively, and t is time.

In one example, the wavefield absorption boundary equation is obtained by: performing Fourier transform on the elastic wave equation to obtain the elastic wave equation after Fourier transform; factorizing and simplifying the elastic wave equation after Fourier transformation to obtain a simplified equation; and converting the simplified equation into a time domain to obtain a wave field absorption boundary equation.

In one example, the elastic wave equation after fourier transform is:

in one example, the simplified equation is:

in one example, the wave field values at respective time instants corresponding to the grid points on the boundary of the velocity model are calculated by equations (4a) - (4d),

the absorption boundary equation corresponding to the bottom edge of the velocity model is as follows:

the absorption boundary equation corresponding to the top edge of the velocity model is:

the absorption boundary equation corresponding to the left side of the velocity model is:

the absorption boundary equation corresponding to the right side of the velocity model is:

wherein the content of the first and second substances,u and w are respectively horizontal displacement and vertical displacement corresponding to the elastic wave field; r₁，r₂、r₃，r₄is wave field absorption coefficient, 1 is more than or equal to r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄More than or equal to 0, and alpha and beta are longitudinal wave speed and transverse wave speed respectively.

In one example, wave field values at respective time instants corresponding to four corner points of the velocity model are calculated by equations (5a) - (5d),

the wavefield calculation equation for the lower right corner point is:

the wave field calculation equation of the corner point at the upper left corner is as follows:

the wave field calculation equation of the corner point at the upper right corner is as follows:

the wavefield calculation equation for the lower left corner point is:

wherein:u and w are respectively a horizontal displacement field and a vertical displacement field corresponding to the elastic wave field;

alpha and beta are respectively the longitudinal wave speed and the transverse wave speed, and X and Z are coordinates along an X axis and a Z axis under the condition of a two-dimensional medium;

wherein:

θ is the angle of incidence of the wave.

The system effectively eliminates the reflection energy generated by the limited speed model boundary by continuously adjusting the wave field absorption coefficient, and can completely absorb the reflection energy from the model boundary when the absorption coefficient is adjusted to the optimal state, so that the facticity of the elastic wave propagating in the two-dimensional isotropic medium can be accurately reflected by the model wave field value obtained by calculation, and the aim of improving the accuracy of the elastic wave numerical simulation result is fulfilled.

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.