CN105093265A - Method of simulating the transmission rules of seismic waves in a TI medium - Google Patents

Abstract

The invention provides a method of simulating the transmission rules of seismic waves in a TI medium, and belongs to the seismic prospecting field. The method includes the steps: firstly, determining a six-order matrix coefficient on the basis of a to-be-solved problem; secondly, solving the first-order speed, i.e., a stress equation; thirdly, solving a convolution differential operator in the second step; fourthly, calculating a window function coefficient, and conducting windowing cross cutting on the convolution differential operator; fifthly, setting the parameter, the epicenter parameter and the epicenter frequency of a model, and carrying out discretization on the model according to geological bodies; and sixthly, conducting time iteration, obtaining the speed and the stress value within a whole area, and ending a simulation process and outputting a wave field snapshot and an earthquake record if the setting time is up.

Description

A kind of seismic wave is in analog in the method for TI Propagation rule

Technical field

The invention belongs to field of seismic exploration, be specifically related to the method for a kind of seismic wave in analog in TI Propagation rule, in transverse isotropy (TI) Propagation rule, for the forward simulation of wave equation, seismic wave propagation law in media as well in analog can be carried out by seismic wave in analog.

Background technology

The propagation law of Study of Seismic ripple in complex dielectrics is an important content in geophysics, is the effective simulation means helping people to be familiar with ball medium.

Seismic wave numerical modeling carrys out the propagation law of modeling effort seismic event in the various medium in underground according to known underground structure and physical parameter, be widely used in seismic prospecting and earthquake field, and the anisotropy of actual underground medium is ubiquitous. describe underground medium wave propagation process and can use wave equation, therefore based on the numerical solution of wave equation as method of finite difference, finite element method, pseudo-spectrometry, spectral element methods etc. are all widely used in anisotropic medium simulation. and these methods have respective relative merits: method of finite difference theoretical foundation is for launch based on Taylor, computing velocity is fast, efficiency is high, but it is more difficult for irregular codes problem, finite element method is based on Theory of Variational Principles, and simulation precision is high, can process the geologic body of complicated shape, but computation process is complicated, and calculated amount is large, pseudo-spectrometry is also one of numerical solution solving partial differential equation, quick Fourier transformation is introduced in computation process, precision can reach the precision of spectrum, but needs global information because Fourier converts, and pseudo-spectrometry seems unable to do what one wishes when processing wave field propagation problem in complex dielectrics, the spectral element method that finite element and pseudo-spectrometry advantages are got up also is also existed for example border and new method and be difficult to the problems such as coupling.

Summary of the invention

The object of the invention is to solve the difficult problem existed in above-mentioned prior art, a kind of convolution differential method that can be used for TI medium forward simulation is provided, describes ripple propagation characteristic in media as well better, and reduce numerical solidification phenomenon in Seismic wave numerical modeling.

The present invention is achieved by the following technical solutions:

Seismic wave is in a method for TI Propagation rule in analog, comprising:

Step 1: determine six rank matrix coefficients according to problem to be asked for:

Step 2: solve one-order velocity---stress equation;

Step 3: ask for the Convolution Differential Operator in step 2;

Step 4: ask for window function coefficient, carries out windowing to Convolution Differential Operator and blocks;

Step 5: base area plastid, setting models parameter, focal shock parameter, focus frequency, carry out discretize to model;

Step 6: carry out time iteration, obtains the speed in whole region and stress value, if setting-up time arrives, then simulation process terminates, and exports wave field snapshot and seismologic record.

Described step 1 is achieved in that

Be six rank matrix forms by TI medium anisotropy Parametric Representation:

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right],$

Wherein C ₁₂=C ₁₁-2C ₆₆. in formula: C is Anisotropic Parameters. in isotropy situation, this six rank matrix is weakened as to describe medium with Lame's constant λ and modulus of shearing μ, now meets following relational expression between each parameter: C ₁₁=C ₃₃=λ+2 μ, C ₁₃=λ, C ₅₅=μ.

Described step 2 is achieved in that

Consider the two-dimensional problems of x-z plane, speed-stress equation is

$\left(\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}\right)=C_{11}{D}_x{v}_x+C_{13}{D}_z{v}_z---(1-a)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}{D}_x{v}_x+C_{33}{D}_z{v}_z---(1-b)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{55}(D_z{v}_x+D_x{v}_z)---(1-c)$

$\frac{{\∂v}_x}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xx}}+D_z{\mathrm{\σ}}_{\mathrm{xz}})---(1-d)$

$\frac{{\∂v}_z}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xz}}+D_z{\mathrm{\σ}}_{\mathrm{zz}})---(1-e)$

(6)

In formula, σ _xx, σ _xz, σ _zzfor the components of stress; v _x, v _zbe respectively velocity level's component and vertical component, ρ is density, D _x, D _zbe respectively the Convolution Differential Operator of space along horizontal and vertical direction.

Described step 3 is achieved in that

Function u (x) is expressed as the convolution form of singular kernel function and meta-function:

$u\left(x\right)={\∫}_{-\∞}^{\∞}\mathrm{\δ}(x-t)u\left(t\right)\mathrm{dt}---\left(2\right)$

The δ function in (2) formula is approached with Shannon singular kernel function:

${\mathrm{\δ}}_{\mathrm{\Δx}}\left(x\right)=\frac{\sin (\mathrm{\πx}/\mathrm{\Δx})}{\mathrm{\πx}/\mathrm{\Δx}}---\left(7\right)$

Wherein, Δ x is spatial mesh size. when Δ x → 0, gained function is Delta function.

For the discrete series u (x of a function _m), by the differentiate of singular kernel function, obtain this function in the interpolation of the total space and differential form:

Wherein represent the lattice point closest to x, m represents the number of discrete point, q representation space differential order, W is operator half width. according to formula (2), (3), (4), the single order Convolution Differential Operator obtained based on Shannon singular kernel is:

$d_1\left(\mathrm{m\Δx}\right)=\left\{\begin{array}{cc}{\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left(\mathrm{m\Δx}\right),& m=\±1,\±2,...\\ 0,& m=0\end{array}\right.---\left(9\right)$

Its staggered-mesh form is:

${\hat{d}}_s\left(\mathrm{m\Δx}\right)={\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left[(m+1/2)\mathrm{\Δx}\right],m=-W,-W+1,...,W-1---\left(10\right)$

The Convolution Differential Operator D of what formula (5) and (6) obtained the be discrete form in formula (1) _xand D _z, wherein formula (5) is applicable to the Seismic wave numerical modeling of common grid, and formula (6) is applicable to the Seismic wave numerical modeling in staggered-mesh situation.

Described step 6 is achieved in that the discrete value of formulae discovery stress and the speed provided by step 2: wherein in formula (1), (1-a), (1-b), (1-c) are by speed component v _x, v _zask for components of stress σ _xx, σ _xz, σ _zzformula, (1-d), (1-e) are by components of stress σ _xx, σ _xz, σ _zzask for speed component v _x, v _zformula, it is 0 that the initial value of these variablees is all composed, and the focal shock parameter in step 5 is relevant with the time, along with time iteration, by the iteration of these five formula, obtain the speed in whole region and stress value, if setting-up time arrives, then simulation process terminates, and exports wave field snapshot and seismologic record.

Compared with prior art, the invention has the beneficial effects as follows:

(1) be more conducive to describing ripple propagation characteristic in media as well;

(2) numerical solidification phenomenon in Seismic wave numerical modeling can be reduced.

Accompanying drawing explanation

Fig. 1 is TI dielectric space staggered-mesh schematic diagram.

Fig. 2-1 optimizes convolution differential method horizontal component.

Fig. 2-2 optimizes convolution differential method vertical component.

Fig. 2-3 is method of finite difference horizontal components.

Fig. 2-4 is method of finite difference vertical components.

Fig. 3 is that finite difference method and convolution differential method are positioned at the seismic record comparison figure (x is horizontal component, and z is vertical component) at (500m, 400m) place at wave detector.

Fig. 4 is the step block diagram of this method.

Embodiment

Below in conjunction with accompanying drawing, the present invention is described in further detail:

The space differentiation of wave equation can adopt various ways, and except above-described finite difference, finite element, pseudo-spectrometry, spectral element method etc., the present invention proposes a kind of Convolution Differential Operator and carries out differential to the space of wave equation.Early stage at 20 century 70s, just there is scholar to propose the concept of Convolution Differential Operator, but apply it in Simulation of Seismic Wave more late.Holberg uses Fourier conversion, designs, and achieve good result to Convolution Differential Operator in three dimensional elasticity ripple numerical simulation.Zhou etc., Zhang Zhongjie etc. also conduct in-depth research this method with Dai Zhiyang etc. and develop.Long Guihua etc., on the working foundation of forefathers, according to function distribution theory, introduce singular kernel Convolution Differential Operator, carry out convolution, ask for wave field variable about space differentiation by one group of differentiating operator through windowing process and wave field variable.Convolution Differential Operator is compared with pseudo-spectrometry, it is advantageous that Convolution Differential Operator is the short operator of a kind of optimization, more can give prominence to the local attribute of space differentiation itself; And pseudo-spectrometry must utilize global information when computing differential, with regard to the propagation problem of complex heterogeneous medium medium wave, run counter to physics law of causality.Convolution Differential Operator is compared with finite difference operator, and precision is higher, can make up the deficiency that finite difference operator describes complex dielectrics.Therefore Convolution Differential Operator be applied to the spatial spreading of wave equation or space differentiation approximate be very effective.

As shown in Figure 4, the present invention solves by the following technical programs, to simulate the propagation law of TI medium medium wave, specifically comprises following steps:

Step 1: determine six rank matrix coefficients according to problem to be asked for:

TI medium anisotropy parameter can be expressed as six rank matrix forms:

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right],$

Wherein C ₁₂=C ₁₁-2C ₆₆. in formula: C is Anisotropic Parameters. in isotropy situation, this six rank matrix can weaken as describing medium with Lame's constant λ (elastic modulus) and μ (modulus of shearing), now meets following relational expression between each parameter: C ₁₁=C ₃₃=λ+2 μ, C ₁₃=λ, C ₅₅=μ.

Step 2: solve one-order velocity---stress equation:

Consider the two-dimensional problems of x-z plane, speed-stress equation is

$\left(\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}\right)=C_{11}{D}_x{v}_x+C_{13}{D}_z{v}_z---(1-a)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}{D}_x{v}_x+C_{33}{D}_z{v}_z---(1-b)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{55}(D_z{v}_x+D_x{v}_z)---(1-c)$

$\frac{{\∂v}_x}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xx}}+D_z{\mathrm{\σ}}_{\mathrm{xz}})---(1-d)$

$\frac{{\∂v}_z}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xz}}+D_z{\mathrm{\σ}}_{\mathrm{zz}})---(1-e)$

(11)

In formula, σ _xx, σ _xz, σ _zzfor the components of stress; v _x, v _zbe respectively velocity level's component and vertical component, ρ is density, D _x, D _zbeing respectively the Convolution Differential Operator of space along horizontal and vertical direction. speed-stress equation that the present invention proposes and conventional speeds-stress equation difference are: utilize convolution operation to instead of derivative operation in conventional speeds-stress equation.(1-a) is example (similar in following formula) with the formula, D _xv _xrepresent v _xspatially carry out convolution operation in the horizontal direction, in conventional speeds-stress equation, need v _xspatially differentiate in the horizontal direction, is then designated as with mathematic(al) representation the method that the present invention proposes avoids the derivative operation that sky is asked.

Formula (1) is the relation between stress and speed, if when service routine carrys out physical simulation calculating, in above-mentioned formula, all variablees all get discrete scheme.The left side is component to time differentiate, and under discrete scheme, available common time difference calculates.

Numerical simulation can be carried out on common grid or staggered-mesh.According to common net case form, these variablees all get the numerical value on net point; According to staggered-mesh form, these variablees some need adopt numerical value on net point, some need adopt the numerical value on half net point (namely getting the mean value of net point around), in this example, each parameter distribution situation of TI medium as shown in Figure 1, wherein the intersection point of solid line and solid line is net point, the intersection point of solid line and dotted line is half net point (value is the mean value of former and later two points of this point), and the intersection point of dotted line and dotted line is also half net point (value is the mean value of four points around this point).

Step 3: ask for the Convolution Differential Operator in step 2.

Circular is:

By function u (x) according to function distribution theory, the convolution form of singular kernel function (desirable Delta (δ) function) and meta-function (desirable u (x) itself) can be expressed as:

$u\left(x\right)={\∫}_{-\∞}^{\∞}\mathrm{\δ}(x-t)u\left(t\right)\mathrm{dt}---\left(2\right)$

And the δ function in (2) formula can approach with Shannon singular kernel function:

${\mathrm{\δ}}_{\mathrm{\Δx}}\left(x\right)=\frac{\sin (\mathrm{\πx}/\mathrm{\Δx})}{\mathrm{\πx}/\mathrm{\Δx}}---\left(12\right)$

Wherein, Δ x is spatial mesh size. when Δ x → 0, gained function is Delta function.

For the discrete series u (x of a function _m), by the differentiate of singular kernel function, can this function in the interpolation of the total space and differential form:

Wherein represent the lattice point closest to x, m represents the number of discrete point, q representation space differential order, W is operator half width., according to formula (2), (3), (4), can be based on the single order Convolution Differential Operator of Shannon singular kernel:

$d_1\left(\mathrm{m\Δx}\right)=\left\{\begin{array}{cc}{\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left(\mathrm{m\Δx}\right),& m=\±1,\±2,...\\ 0,& m=0\end{array}\right.---\left(14\right)$

Its staggered-mesh form is:

${\hat{d}}_s\left(\mathrm{m\Δx}\right)={\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left[(m+1/2)\mathrm{\Δx}\right],m=-W,-W+1,...,W-1---\left(15\right)$

The Convolution Differential Operator D of what formula (5) and (6) obtained the be discrete form in formula (1) _x(horizontal direction) and D _z(vertical direction), wherein formula (5) is applicable to the Seismic wave numerical modeling of common grid, and formula (6) is applicable to the Seismic wave numerical modeling in staggered-mesh situation.

Step 4: ask for window function coefficient, windowing is carried out to differentiating operator and blocks:

In order to eliminate the Gibbs phenomenon because interruption operator causes, windowing process can be carried out to operator.For staggered-mesh situation, windowing is carried out to formula (6):

$\stackrel{\‾}{d_s}\left(\mathrm{m\Δx}\right)={\hat{d}}_s\left(\mathrm{m\Δx}\right)w(m+1/2),m=-W,-W+1,...,W-1---\left(16\right)$

Institute's windowed function is Hanning window:

$w\left(m\right)={[2\mathrm{\α}-1+2(1-\mathrm{\α}){\cos }^2\frac{\mathrm{m\π}}{2(W+2)}]}^{\frac{\mathrm{\β}}{2}}---\left(17\right)$

Wherein, α, β are the parameter of portraying window function character. when operator half width W mono-timing, the precision of numerical simulation can be regulated further by regulating the factor alpha of window function and β. about α, β to ask for mathematical derivation process comparatively complicated, less with relevance degree of the present invention, do not repeating herein.In this example, operator length W is 9, is α=0.4015 with the window function weight coefficient that optimization method (as simulated annealing method) is tried to achieve, β=2.7069.Can certainly according to target problem need select suitable operator length, ask for new window function weight coefficient by optimization method.

Step 5: base area plastid, carries out discretize to model.Setting models parameter, focal shock parameter, focus frequency.

In this example, setting model parameter is C ₁₁=25.5 × 10 ⁹n.m, C ₁₃=1.0 × 10 ⁹n.m, C ₃₃=18.4 × 10 ⁹n.m, C ₅₅=5.6 × 10 ⁹n.m, density is 2440kg/m ³. it is 301 × 301 that model meshes is counted, sizing grid is Δ x=Δ z=10m, time sampling interval is 1ms, focus is positioned at model centre, coordinate is (1500m, 1500m). focus adopts vertical concentrated force source, is the Ricker wavelet of dominant frequency 25Hz, the reception arrangement of seismologic record is positioned at underground 500m depth level arrangement. and above parameter only need meet the physical property of plastid routinely, and the span of the present invention to above parameter does not limit.

This step and routine fluctuations equation the Forward Modeling similar.

Step 6: the formulae discovery stress provided by step 2 and the discrete value of speed.Wherein in formula (1), (1-a), (1-b), (1-c) are by speed component v _x, v _zask for components of stress σ _xx, σ _xz, σ _zzformula, (1-d), (1-e) are by components of stress σ _xx, σ _xz, σ _zzask for speed component v _x, v _zformula.It is 0 that the initial value of these variablees is all composed, and in step 5, focal shock parameter is relevant with the time, along with time iteration, by the iteration of these five formula, obtains the speed in whole region and stress value.If setting-up time arrives, then simulation process terminates, and exports wave field snapshot and seismologic record.

The wave field snapshot of elastic wave propagation when Fig. 2-1 to Fig. 2-4 is this medium 400ms, wherein Fig. 2-1, Fig. 2-2 is 9 analog results optimizing Convolution Differential Operator method, the former is horizontal component, the latter is vertical component, Fig. 2-3, 2-4 is the analog result of identical operator length (8 rank) method of finite difference, the former is horizontal component, the analog result that the latter optimizes the analog result of Convolution Differential Operator method and 8 rank staggering mesh finite-difference for vertical component .9 point is basically identical, difference is that comparatively significantly dispersion phenomenon has appearred in staggering mesh finite-difference method when qSV wave simulation, and frequency dispersion obtains good compacting when simulating with optimization Convolution Differential Operator. this shows, feasible by the simulation that this optimization Convolution Differential Operator is used for TI medium Elastic Wave, and precision is high, the method of finite difference of identical operator length is better than in compacting numerical solidification effect. therefore can strengthen mesh spacing as required in computation process thus reduce calculated amount, and then raising counting yield.

For the precision of further analysis optimization Convolution Differential Operator method, the position of receiver is have chosen for (500m in example, 400m) the single-channel seismic record at place, be illustrated in figure 3 finite difference method and convolution differential method is positioned at (500m at wave detector, (x is horizontal component to the seismic record comparison figure at 400m) place, z is vertical component), in figure, solid line is the analog result of 8 rank method of finite difference, dotted line is 9 analog results optimizing Convolution Differential Operator method, x is horizontal component, z is vertical component. can find out, except except the crest of amplitude and two kinds, trough place method slightly deviation, other place all coincide better. in addition, after 0.9s staggered-mesh method of finite difference seismogram on there is obvious dispersion phenomenon, and the record figure optimizing Convolution Differential Operator method does not observe frequency dispersion effect. the wave field snapshot in this and prior figures 2-1 to Fig. 2-4 matches, demonstrate and optimize the finite difference operator method of Convolution Differential Operator method simulation precision higher than same operator length. and in realistic simulation, the coefficient regulating window function neatly can also be needed according to precision, thus regulate the precision of Convolution Differential Operator. this explanation, this Convolution Differential Operator method except precision high, also have and use feature flexibly.

This method is by being optimized the optimizing process achieved the staggered-mesh Convolution Differential Operator based on Shannon singular kernel theory to window function parameter. and apply this optimization Convolution Differential Operator method and numerical simulation has been carried out to anisotropic medium, analyze the propagation characteristic of elastic wave in this type of medium, and contrast with high-order staggering mesh finite-difference method. numerical experiment results shows, the method is applicable to anisotropic medium Elastic Wave field stimulation, precision is high, good stability is a kind of Effective Numerical method studying seismic wave propagation in complex dielectrics.

Technique scheme is one embodiment of the present invention, for those skilled in the art, on the basis that the invention discloses application process and principle, be easy to make various types of improvement or distortion, and the method be not limited only to described by the above-mentioned embodiment of the present invention, therefore previously described mode is just preferred, and does not have restrictive meaning.

Claims (5)

1. seismic wave, in a method for TI Propagation rule, is characterized in that in analog: described method comprises:

Step 1: determine six rank matrix coefficients according to problem to be asked for:

Step 2: solve one-order velocity---stress equation;

Step 3: ask for the Convolution Differential Operator in step 2;

Step 4: ask for window function coefficient, carries out windowing to Convolution Differential Operator and blocks;

Step 5: base area plastid, setting models parameter, focal shock parameter, focus frequency, carry out discretize to model;

Step 6: carry out time iteration, obtains the speed in whole region and stress value, if setting-up time arrives, then simulation process terminates, and exports wave field snapshot and seismologic record.

2. seismic wave in analog according to claim 1 is in the method for TI Propagation rule, it is characterized in that: described step 1 is achieved in that

Be six rank matrix forms by TI medium anisotropy Parametric Representation:

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right],$

Wherein C ₁₂=C ₁₁-2C ₆₆. in formula: C is Anisotropic Parameters. in isotropy situation, this six rank matrix is weakened as to describe medium with Lame's constant λ and modulus of shearing μ, now meets following relational expression between each parameter: C ₁₁=C ₃₃=λ+2 μ, C ₁₃=λ, C ₅₅=μ.

3. seismic wave in analog according to claim 2 is in the method for TI Propagation rule, it is characterized in that: described step 2 is achieved in that

Consider the two-dimensional problems of x-z plane, speed-stress equation is

$\left(\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}\right)=C_{11}{D}_x{v}_x+C_{13}{D}_z{v}_z---(1-a)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}{D}_x{v}_x+C_{33}{D}_z{v}_z---(1-b)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{55}(D_z{v}_x+D_x{v}_z)---(1-c)$

$\frac{{\∂v}_x}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xx}}+D_z{\mathrm{\σ}}_{\mathrm{xz}})---(1-d)$

$\frac{{\∂v}_z}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xz}}+D_z{\mathrm{\σ}}_{\mathrm{zz}})---(1-e)$

(1)

In formula, σ _xx, σ _xz, σ _zzfor the components of stress; v _x, v _zbe respectively velocity level's component and vertical component, ρ is density, D _x, D _zbe respectively the Convolution Differential Operator of space along horizontal and vertical direction.

4. seismic wave in analog according to claim 3 is in the method for TI Propagation rule, it is characterized in that: described step 3 is achieved in that

Function u (x) is expressed as the convolution form of singular kernel function and meta-function:

$u\left(x\right)={\∫}_{-\∞}^{\∞}\mathrm{\δ}(x-t)u\left(t\right)\mathrm{dt}---\left(2\right)$

The δ function in (2) formula is approached with Shannon singular kernel function:

${\mathrm{\δ}}_{\mathrm{\Δx}}\left(x\right)=\frac{\sin (\mathrm{\πx}/\mathrm{\Δx})}{\mathrm{\πx}/\mathrm{\Δx}}---\left(2\right)$

Wherein, Δ x is spatial mesh size. when Δ x → 0, gained function is Delta function.

For the discrete series u (x of a function _m), by the differentiate of singular kernel function, obtain this function in the interpolation of the total space and differential form:

Wherein represent the lattice point closest to x, m represents the number of discrete point, q representation space differential order, W is operator half width. according to formula (2), (3), (4), the single order Convolution Differential Operator obtained based on Shannon singular kernel is:

$d_1\left(\mathrm{m\Δx}\right)=\left\{\begin{array}{cc}{\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left(\mathrm{m\Δx}\right),& m=\±1,\±2,...\\ 0,& m=0\end{array}\right.---\left(4\right)$

Its staggered-mesh form is:

${\hat{d}}_s\left(\mathrm{m\Δx}\right)={\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left[(m+1/2)\mathrm{\Δx}\right],m=-W,-W+1,...,W-1---\left(5\right)$

The Convolution Differential Operator D of what formula (5) and (6) obtained the be discrete form in formula (1) _xand D _z, wherein formula (5) is applicable to the Seismic wave numerical modeling of common grid, and formula (6) is applicable to the Seismic wave numerical modeling in staggered-mesh situation.

5. seismic wave in analog according to claim 4 is in the method for TI Propagation rule, it is characterized in that: described step 6 is achieved in that the discrete value of formulae discovery stress and the speed provided by step 2: wherein in formula (1), (1-a), (1-b), (1-c) are by speed component v _x, v _zask for components of stress σ _xx, σ _xz, σ _zzformula, (1-d), (1-e) are by components of stress σ _xx, σ _xz, σ _zzask for speed component v _x, v _zformula, it is 0 that the initial value of these variablees is all composed, and the focal shock parameter in step 5 is relevant with the time, along with time iteration, by the iteration of these five formula, obtain the speed in whole region and stress value, if setting-up time arrives, then simulation process terminates, and exports wave field snapshot and seismologic record.