CN104597488A - Optimum design method of finite difference template of non-equiangular long-grid wave equation - Google Patents

Abstract

The invention discloses an optimum design method of a finite difference template of a non-equiangular long-grid wave equation. The method comprises the following steps: according to a time sampling interval and a spatial sampling interval, mesh subdivision is executed to the simulating region of an actual geologic model; according to the maximum appointed permissible error and the wave number range, the corresponding operator lengths of different acoustic velocities are calculated; based on the finite difference method of the least square optimization time space domain and the corresponding operator length, the optimized finite difference coefficient of the grid is acquired; the acquired optimized finite difference coefficient is brought into the wave equation with a difference scheme to perform forward modeling of the wave equation. According to the method, the finite difference method of the time space domain only applied to square and cube grids is expanded into the rectangular or rectangular parallelepiped grids, thus a special simulation precision requirement and a requirement for saving a calculated amount are met in actual production.

Description

Non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method

Technical field

The present invention relates to seismic forward modeling numerical simulation technology field, particularly the non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method of one.

Background technology

Seismic forward modeling numerical simulation technology is when complex underground structure group (comprising isotropic medium, anisotropic medium, the heterogeneous anisotropic medium of Biot, random hole medium etc.) is known, numerical computation method is utilized to make ripple at this Propagation, through the repeatedly transmission of subsurface geological structure, reflection, scattering, by the process that the wave detector of earth's surface or underground receives.Accurate wave equation numerical solution is utilized to carry out the seismic response of simulate formation complex geological structure, for many aspects such as Study of Seismic wave traveling mechanism, the special treatment method of seismic data and the explanation of bad ground provide the mathematical physics foundation of more science.In recent years, Wave Equation Numerical method is widely used in reverse-time migration and full waveform inversion.

Wave equation forward modeling has multiple method, comparatively common are: finite difference method, pseudo-spectrometry, finite element method, boundary element method, spectral element method etc.Wherein finite difference method because its calculated amount is little, counting yield is high, more complicated rate pattern can be adapted to and be widely used.Method of finite difference can be divided into according to different standards: Explicit finite difference and implicit expression finite difference; Regular grid finite difference, staggering mesh finite-difference and rotationally staggered grid finite difference.In method of finite difference, difference coefficient can be tried to achieve by Taylor series expansion or optimization method, respectively corresponding finite difference based on Taylor series expansion and the finite difference based on optimization.In conventional finite method of difference, difference coefficient is obtained by the dispersion relation of minimization spatial domain.In recent years, occurred a kind of time-space domain method of finite difference, the method asks for difference coefficient by the dispersion relation of minimization time domain and spatial domain, has higher simulation precision and better stability.

Current time-space domain finite difference method requires that the spatial sampling interval in all directions is equal, and namely needing model facetization is square or square grid.And in actual production, in order to meet specific accuracy requirement or in order to save calculated amount, we usually need to be rectangle or rectangular parallelepiped grid (normally the mesh spacing of depth direction is different from horizontal direction) by model facetization.Be applicable to the time-space domain finite difference method of square and square grid, specific accuracy requirement can not be met and maybe can not save calculated amount.

Summary of the invention

Embodiments provide a kind of non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method, the time-space domain finite difference method being only applicable to square and square grid is expanded in rectangle or rectangular parallelepiped grid, meet the demand of specific simulation precision requirement and saving calculated amount in actual production, the method comprises:

According to time sampling interval and spatial sampling interval, mesh generation is carried out to the simulated domain of actual geologic model;

According to the given limits of error and wave-number range, corresponding operator length is obtained to different acoustic velocity;

Based on the time-space domain method of finite difference of Least-squares minimization and the operator length of described correspondence, obtain the optimization finite difference coefficient of grid;

The optimization finite difference coefficient obtained is brought into the wave equation of difference scheme, carry out Wave equation forward modeling.

In one embodiment, describedly according to time sampling interval and spatial sampling interval, difference gridding subdivision is carried out to the simulated domain of actual geologic model, comprising: the simulated domain of actual geologic model is split into rectangular parallelepiped grid or rectangular node;

Based on the time-space domain method of finite difference of Least-squares minimization and the operator length of described correspondence, obtain the optimization finite difference coefficient of grid, comprising: obtain the optimization finite difference coefficient of rectangular parallelepiped grid or the optimization finite difference coefficient of rectangular node.

In one embodiment, the optimization finite difference coefficient of described acquisition rectangular parallelepiped grid, according to following formulae discovery:

Wherein,

Wherein, b is wave number, and M is operator length, a _mfor the finite difference coefficient after optimization, θ is the angle of plane wave propagation direction and surface level, θ ∈ [0, π]; φ is the position angle of plane wave propagation, φ ∈ [0,2 π]; v is acoustic velocity, and τ is time sampling interval, and h is x, y direction sampling interval; C=Δ z/h, c are parametric variable, and Δ z is z direction sampling interval; β=kh, k are parametric variable, and β is wave-number range, β ∈ [0, b]; M is parametric variable, and m is integer, m ∈ [1, M]; N is parametric variable, and n is integer, n ∈ [1, M].

In one embodiment, the maximum error that the optimization finite difference coefficient of described rectangular parallelepiped grid is corresponding meets following constraint condition:

ξ _1max<η；

Wherein, ξ _1maxfor the maximum error that the optimization finite difference coefficient of rectangular parallelepiped grid is corresponding; η is that maximum permission mistake 5 is poor.

In one embodiment, the maximum error ξ that the optimization finite difference coefficient of described rectangular parallelepiped grid is corresponding _1maxbe calculated as follows:

${\mathrm{\&xi;}}_{1\max }=\underset{\mathrm{\&beta;}\&Element;[0,b],\mathrm{\&theta;}\&Element;\left[\mathrm{0,2}\mathrm{\&pi;}\right]}{\max }\left|{\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})\right|;$

Wherein,

${\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{\mathrm{r\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\begin{array}{c}\left[\begin{array}{c}\mathrm{si}{n}^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\right)-\mathrm{si}{n}^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\right)\\ -\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\&beta;}}{2}\sin \mathrm{\&theta;}\right)\end{array}\right]\end{array}}-1.$

In one embodiment, the described wave equation optimization finite difference coefficient obtained being brought into difference scheme, carry out Wave equation forward modeling, comprise: the optimization finite difference coefficient of the rectangular parallelepiped grid obtained is substituted in the three-dimensional acoustic wave wave equation of difference scheme, carries out three-dimensional acoustic wave Wave equation forward modeling; The three-dimensional acoustic wave wave equation of described difference scheme is:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,y,z}^{t-1}-2P_{x,y,z}^t+P_{x,y,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,y,z-m}^t+P_{x,y,z+m}^t)]\\ +\frac{1}{h^2}[2a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,y,z}^t+P_{x+m,y,z}^t+P_{x,y-m,z}^t+P_{x,y+m,z}^t)]\end{array};$

Wherein, P is acoustic pressure.The optimization finite difference coefficient of described acquisition rectangular node, according to following formulae discovery:

Wherein,

Wherein, b is wave number, and M is operator length, a _mfor the finite difference coefficient after optimization, θ is the angle of plane wave propagation direction and surface level, θ ∈ [0, π]; v is acoustic velocity, and τ is time sampling interval, and h is x direction sampling interval; C=Δ z/h, c are parametric variable, and Δ z is z direction sampling interval; β=kh, k are parametric variable, and β is wave-number range, β ∈ [0, b]; M is parametric variable, and m is integer, m ∈ [1, M]; N is parametric variable, and n is integer, n ∈ [1, M].

In one embodiment, the maximum error that the optimization finite difference coefficient of described rectangular node is corresponding meets following constraint condition:

ξ _2max<η；

Wherein, ξ _2maxfor the maximum error that the optimization finite difference coefficient of rectangular node is corresponding; η is the limits of error.

In one embodiment, the maximum error ξ that the optimization finite difference coefficient of described rectangular node is corresponding _2maxbe calculated as follows:

${\mathrm{\&xi;}}_{2\max }=\underset{\mathrm{\&beta;}\&Element;[0,b],\mathrm{\&theta;}\&Element;\left[\mathrm{0,2}\mathrm{\&pi;}\right]}{\max }\left|{\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})\right|;$

Wherein, ${\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{\mathrm{r\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m[\mathrm{si}{n}^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\right)+\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\&beta;}}{2}\sin \mathrm{\&theta;}\right)]}-1.$

In one embodiment, described the optimization finite difference coefficient obtained is brought in the wave equation of difference scheme, carry out Wave equation forward modeling, also comprise: the optimization finite difference coefficient of the rectangular node obtained is substituted in the two-dimentional Acoustic Wave-equation of difference scheme, carries out two-dimentional Acoustic Wave-equation forward simulation; The two-dimentional Acoustic Wave-equation of described difference scheme is:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,z}^{t-1}-2P_{x,z}^t+P_{x,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,z-m}^t+P_{x,z+m}^t)]\\ +\frac{1}{h^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,z}^t+P_{x+m,z}^t)]\end{array};$

Wherein, P is acoustic pressure.

In embodiments of the present invention, propose a kind of non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method, the time-space domain finite difference method being only applicable to square and square grid extend in rectangle or rectangular parallelepiped grid by the method, meets the demand of specific simulation precision requirement and saving calculated amount in actual production.

Accompanying drawing explanation

Accompanying drawing described herein is used to provide a further understanding of the present invention, forms a application's part, does not form limitation of the invention.In the accompanying drawings:

Fig. 1 is the non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method flow diagram of one that the embodiment of the present invention provides;

Fig. 2 is that the graph of errors obtained by traditional finite difference method that provides of the embodiment of the present invention is with the Changing Pattern figure of the direction of propagation;

Fig. 3 is the Changing Pattern figure of the graph of errors obtained by non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method that provides of the embodiment of the present invention with the direction of propagation;

Fig. 4 is that the graph of errors obtained by traditional finite difference method that provides of the embodiment of the present invention is with the Changing Pattern figure of operator length;

Fig. 5 is the Changing Pattern figure of the graph of errors obtained by non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method that provides of the embodiment of the present invention with operator length;

Fig. 6 is the wave field snapshot comparison diagram of the even sound wave medium obtained respectively by traditional finite difference method and non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method that provides of the embodiment of the present invention in the 1.0s moment;

Fig. 7 is the comparison of wave shape figure of dotted line position in Fig. 6;

Fig. 8 is the wave field snapshot comparison diagram of the Marmousi model obtained respectively by traditional finite difference method and non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method that provides of the embodiment of the present invention in the 1.0s moment;

Embodiment

For making the object, technical solutions and advantages of the present invention clearly understand, below in conjunction with embodiment and accompanying drawing, the present invention is described in further details.At this, exemplary embodiment of the present invention and illustrating for explaining the present invention, but not as a limitation of the invention.

Inventor finds, current time-space domain finite difference method requires that the spatial sampling interval in all directions is equal, namely needing model facetization is square or square grid, but this kind of method can not meet specific accuracy requirement maybe can not save calculated amount.If be rectangle or rectangular parallelepiped grid by model facetization, then can meet specific accuracy requirement or save calculated amount.Based on this, the present invention proposes a kind of non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method.

Fig. 1 is the non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method flow diagram of one that the embodiment of the present invention provides, and as shown in Figure 1, the method comprises:

Step 101: mesh generation is carried out to the simulated domain of actual geologic model according to time sampling interval and spatial sampling interval;

Step 102: according to the given limits of error and wave-number range, obtains corresponding operator length to different acoustic velocity;

Step 103: based on the time-space domain method of finite difference of Least-squares minimization and the operator length of described correspondence, obtains the optimization finite difference coefficient of grid;

Step 104: the wave equation optimization finite difference coefficient obtained being brought into difference scheme, carries out Wave equation forward modeling.

During concrete enforcement, need to be rectangular node or rectangular parallelepiped grid (normally the mesh spacing of depth direction is different from horizontal direction) by the simulated domain of actual geologic model (or zoning) subdivision.

Three-dimensional acoustic wave wave equation is as follows:

$\frac{{\&PartialD;}^{2}P}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}P}{\&PartialD;y^2}+\frac{{\&PartialD;}^{2}P}{\&PartialD;z^2}=\frac{1}{V^2}\frac{{\&PartialD;}^{2}P}{\&PartialD;t^2}---\left(1\right)$

In formula, P represents acoustic pressure, and V represents acoustic velocity.

Be different from the situation in x, y direction for z direction mesh spacing, the optimization finite difference scheme of concrete rectangular parallelepiped grid is as follows:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,y,z}^{t-1}-2P_{x,y,z}^t+P_{x,y,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,y,z-m}^t+P_{x,y,z+m}^t)]\\ +\frac{1}{h^2}[2a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,y,z}^t+P_{x+m,y,z}^t+P_{x,y-m,z}^t+P_{x,y+m,z}^t)]\end{array}---\left(2\right)$

Wherein, τ is time sampling interval; H is x, y direction sampling interval; C=Δ z/h, c are parametric variable, and Δ z is z direction sampling interval; a _mfor the finite difference coefficient after optimization, M represents operator length; M is parametric variable, and m is integer, m ∈ [1, M].

In order to contracted notation, we will be abbreviated as according to Plane wave theory, Wo Menling:

$P_{m,l,j}^n=e^{i[k_x(x+\mathrm{mh})+k_y(y+\mathrm{lh})+k_z(z+\mathrm{jch})-\mathrm{\&omega;}(t+\mathrm{n\&tau;})]}---\left(3\right)$

Wherein,

$\left\{\begin{array}{c}{k}_x\text{=kcos}\left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&phi;}\right)\\ k_y=k\cos \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&phi;}\right)\\ k_z=k\sin \left(\mathrm{\&theta;}\right)\end{array}\right.---\left(4\right)$

Wherein, i, l, j, n are parametric variable, and n is integer, n ∈ [1, M]; K is parametric variable; The angle that θ ∈ [0, π] is plane wave propagation direction and surface level, the position angle that φ ∈ [0,2 π] is plane wave propagation.Equation (3), (4) are substituted into equation (2), and suitable abbreviation, can obtain:

$\begin{array}{c}2r^{-2}[\cos \left(\mathrm{\&omega;\&tau;}\right)-1]=\frac{1}{c^2}{a}_0+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{2}{c^2}{a}_m\left[\cos \left(m{k}_{z}h\right)\right]\\ +2a_0+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}2a_m[\cos \left(m{k}_{x}h\right)+\cos \left(m{k}_{y}h\right)]\end{array}---\left(5\right)$

Wherein, again according to constraint condition;

$a_0+2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m=0---\left(6\right)$

Can obtain:

$r^{-2}[1-\cos \left(\mathrm{\&omega;\&tau;}\right)]=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\left[\begin{array}{c}2-\cos \left(\mathrm{m\&beta;}\cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\right)\text{-cos}\left(\mathrm{m\&beta;}\cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\right)\\ +\frac{1}{c^2}[1-\cos \left(\mathrm{mc\&beta;}\sin \mathrm{\&theta;}\right)]\end{array}\right]---\left(7\right)$

Wherein, β=kh, β are wave-number range, β ∈ [0, b].Optimization finite-difference algorithm will be obtained a set of difference coefficient exactly and make equation (7) two ends error minimum.According to Least Square Theory, this problem can be converted into and solve following system of equations:

Wherein,

By solving system of linear equations (8), we can obtain the optimization difference coefficient of rectangular parallelepiped grid.Again difference coefficient is updated in equation (2), just can carries out solving of three-dimensional acoustic wave wave equation.

In order to ensure the validity of the inventive method, need the constraint condition added in the optimization procedure above below, the maximum error that namely the optimization finite difference coefficient of rectangular parallelepiped grid is corresponding meets following constraint condition:

ξ _1max<η (10)

We pass through maximum error corresponding to following formula calculation optimization difference coefficient:

${\mathrm{\&xi;}}_{1\max }=\underset{\mathrm{\&beta;}\&Element;[0,b],\mathrm{\&theta;}\&Element;\left[\mathrm{0,2}\mathrm{\&pi;}\right]}{\max }\left|{\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})\right|---\left(11\right)$

Wherein,

${\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{\mathrm{r\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\begin{array}{c}\left[\begin{array}{c}\mathrm{si}{n}^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\right)-\mathrm{si}{n}^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\right)\\ -\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\&beta;}}{2}\sin \mathrm{\&theta;}\right)\end{array}\right]\end{array}}-1---\left(12\right)$

Wherein, η is the limits of error.

ξ _1maxonly relevant with b and M.When M fixes, ξ _1maxbe decided by b, ξ _1maxincrease along with b and become large.Therefore, if initial b does not meet equation (10), we should by reducing b until ξ gradually _1max< η obtains best b.In addition, when b fixes, ξ _1maxonly relevant with M, ξ _1maxincrease along with M and diminish.So if initial M does not meet equation (10), we should by increasing M until ξ gradually _1max< η obtains best M.

For two-dimensional rectangle grid, be different from the situation in x direction for z direction mesh spacing, concrete optimization finite difference scheme is as follows:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,z}^{t-1}-2P_{x,z}^t+P_{x,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,z-m}^t+P_{x,z+m}^t)]\\ +\frac{1}{h^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,z}^t+P_{x+m,z}^t)]\end{array}---\left(13\right)$

In formula, P represents acoustic pressure, and V represents acoustic velocity; Wherein, τ is time sampling interval; H is x direction sampling interval; C=Δ z/h, c are parametric variable, and Δ z is z direction sampling interval; a _mfor the finite difference coefficient after optimization, M represents operator length; M is parametric variable, and m is integer, m ∈ [1, M].

We can obtain the difference coefficient of the optimization of rectangular node by solving following linear equation:

Wherein,

By solving system of linear equations (15), we can obtain the optimization difference coefficient of rectangular node.Again difference coefficient is updated in equation (13), just can carries out solving of two-dimentional Acoustic Wave-equation.

Same, in order to ensure the validity of the inventive method, need the constraint condition added in the optimization procedure above below, the maximum error that namely the optimization finite difference coefficient of rectangular node is corresponding meets following constraint condition:

ξ _2max<η (16)

We pass through maximum error corresponding to following formula calculation optimization difference coefficient:

${\mathrm{\&xi;}}_{2\max }=\underset{\mathrm{\&beta;}\&Element;[0,b],\mathrm{\&theta;}\&Element;\left[\mathrm{0,2}\mathrm{\&pi;}\right]}{\max }\left|{\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})\right|---\left(17\right)$

Wherein,

${\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{\mathrm{r\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m[\mathrm{si}{n}^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\right)+\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\&beta;}}{2}\sin \mathrm{\&theta;}\right)]}-1---\left(18\right)$

Wherein, η is the limits of error.

For an even sound wave medium, advantage of the present invention is described below.Velocity of longitudinal wave is 1500m/s, and time sampling interval is 1ms, x direction sampling interval be 10m, z direction sampling interval is 5m, operator length M ∈ [2,10].

Regulation b=2.74 and M=8, the inventive method (i.e. non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method) and traditional finite difference method is adopted to obtain the Changing Pattern of graph of errors with the direction of propagation, respectively as shown in Figures 2 and 3.Regulation η=10 ^-2.5≈ 0.003, adopts the inventive method and traditional finite difference method to obtain the Changing Pattern of graph of errors with operator length, respectively as shown in Figures 4 and 5.From Fig. 2 to Fig. 5, compared with the time-space domain finite difference method (traditional finite difference method) based on Taylor, the inventive method has wider effective band and less numerical solidification.Regulation η=10 ^-2.5≈ 0.003, adopts traditional finite difference method and the inventive method to obtain the wave field snapshot comparison diagram of even sound wave medium, as shown in Figure 6.Wherein, the left side figure in Fig. 6 is the wave field snapshot adopting traditional finite difference method to obtain even sound wave medium, and the right figure is the wave field snapshot of the even sound wave medium adopting the inventive method to obtain.Fig. 7 is the comparison of wave shape figure of dotted line position in Fig. 6.From Fig. 6 and Fig. 7, when operator length is identical, the inventive method simulation precision is higher.

In order to further illustrate validity of the present invention, we just drill with Marmousi model.Time sampling interval is 1ms, x direction sampling interval be 12m, z direction sampling interval is 9m.The operator length M=9 of traditional method of finite difference.Time-space domain finite difference adopts the algorithm becoming operator length, selects the operator of different length, M ∈ [2,9] according to different speed.Fig. 8 is the wave field snapshot comparison diagram of the Marmousi model obtained respectively by traditional finite difference method and the inventive method that provides of the embodiment of the present invention in the 1.0s moment, upper graph be the Marmousi model that obtained by traditional finite difference method at the wave field snapshot in 1.0s moment, lower edge graph is the Marmousi model that obtained respectively by the inventive method wave field snapshot in the 1.0s moment; As shown in Figure 8, the inventive method is applicable to the forward simulation of complex dielectrics, and its simulate effect is obviously better than traditional method of finite difference.

In sum, the inventive method has the following advantages: 1. the numerical solidification reducing medium-high frequency section.2. more realistic based on the dispersion relation of time-space domain, simulation precision is higher.3. be applicable to rectangle and rectangular parallelepiped grid, the demand of specific accuracy requirement and saving calculated amount in actual production can be met.

The foregoing is only the preferred embodiments of the present invention, be not limited to the present invention, for a person skilled in the art, the embodiment of the present invention can have various modifications and variations.Within the spirit and principles in the present invention all, any amendment done, equivalent replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (10)

1. a non-equilateral fourdrinier wire lattice wave equation finite difference optimum formwork design method, is characterized in that, comprising:

According to time sampling interval and spatial sampling interval, mesh generation is carried out to the simulated domain of actual geologic model;

According to the given limits of error and wave-number range, corresponding operator length is obtained to different acoustic velocity;

Based on the time-space domain method of finite difference of Least-squares minimization and the operator length of described correspondence, obtain the optimization finite difference coefficient of grid;

The optimization finite difference coefficient obtained is brought into the wave equation of difference scheme, carry out Wave equation forward modeling.

2. the method for claim 1, it is characterized in that, describedly according to time sampling interval and spatial sampling interval, difference gridding subdivision is carried out to the simulated domain of actual geologic model, comprising: the simulated domain of actual geologic model is split into rectangular parallelepiped grid or rectangular node;

Based on the time-space domain method of finite difference of Least-squares minimization and the operator length of described correspondence, obtain the optimization finite difference coefficient of grid, comprising: obtain the optimization finite difference coefficient of rectangular parallelepiped grid or the optimization finite difference coefficient of rectangular node.

3. method as claimed in claim 2, is characterized in that, the optimization finite difference coefficient of described acquisition rectangular parallelepiped grid, according to following formulae discovery:

Wherein,

Wherein, b is wave number, and M is operator length, a _mfor the finite difference coefficient after optimization, θ is the angle of plane wave propagation direction and surface level, θ ∈ [0, π]; φ is the position angle of plane wave propagation, φ ∈ [0,2 π]; v is acoustic velocity, and τ is time sampling interval, and h is x, y direction sampling interval; C=Δ z/h, c are parametric variable, and Δ z is z direction sampling interval; β=kh, k are parametric variable, and β is wave-number range, β ∈ [0, b]; M is parametric variable, and m is integer, m ∈ [1, M]; N is parametric variable, and n is integer, n ∈ [1, M].

4. method as claimed in claim 3, it is characterized in that, the maximum error that the optimization finite difference coefficient of described rectangular parallelepiped grid is corresponding meets following constraint condition:

ξ _1max<η；

Wherein, ξ _1maxfor the maximum error that the optimization finite difference coefficient of rectangular parallelepiped grid is corresponding; η is the limits of error.

5. method as claimed in claim 4, is characterized in that, the maximum error ξ that the optimization finite difference coefficient of described rectangular parallelepiped grid is corresponding _1maxbe calculated as follows:

${\mathrm{\&xi;}}_{1\max }=\underset{\mathrm{\&beta;}\&Element;[0,b],\mathrm{\&theta;}\&Element;\left[\mathrm{0,2}\mathrm{\&pi;}\right]}{\max }\left|{\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})\right|;$

Wherein,

${\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{\mathrm{r\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\left[\begin{array}{c}{\sin }^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\right)-{\sin }^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\right)\\ -\frac{1}{c^2}{\sin }^2\left(\frac{\mathrm{mc\&beta;}}{2}\sin \mathrm{\&theta;}\right)\end{array}\right]}-1.$

6. method as claimed in claim 5, it is characterized in that, the described wave equation optimization finite difference coefficient obtained being brought into difference scheme, carry out Wave equation forward modeling, comprise: the optimization finite difference coefficient of the rectangular parallelepiped grid obtained is substituted in the three-dimensional acoustic wave wave equation of difference scheme, carries out three-dimensional acoustic wave Wave equation forward modeling; The three-dimensional acoustic wave wave equation of described difference scheme is:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,\mathrm{yz}}^{t-1}-{2P}_{x,y,z}^t+P_{x,y,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,y,z-m}^t+P_{x,y,z+m}^t)]\\ +\frac{1}{h^2}[{2a}_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,y,z}^t+P_{x+m,y,z}^t+P_{x,y-m,z}^t+P_{x,y+m,z}^t)]\end{array};$

Wherein, P is acoustic pressure.

7. method as claimed in claim 2, is characterized in that, the optimization finite difference coefficient of described acquisition rectangular node, according to following formulae discovery:

Wherein,

Wherein, b is wave number, and M is operator length, a _mfor the finite difference coefficient after optimization, θ is the angle of plane wave propagation direction and surface level, θ ∈ [0, π]; v is acoustic velocity, and τ is time sampling interval, and h is x direction sampling interval; C=Δ z/h, c are parametric variable, and Δ z is z direction sampling interval; β=kh, k are parametric variable, and β is wave-number range, β ∈ [0, b]; M is parametric variable, and m is integer, m ∈ [1, M]; N is parametric variable, and n is integer, n ∈ [1, M].

8. method as claimed in claim 7, it is characterized in that, the maximum error that the optimization finite difference coefficient of described rectangular node is corresponding meets following constraint condition:

ξ _2max<η；

Wherein, ξ _2maxfor the maximum error that the optimization finite difference coefficient of rectangular node is corresponding; η is the limits of error.

9. method as claimed in claim 8, is characterized in that, the maximum error ξ that the optimization finite difference coefficient of described rectangular node is corresponding _2maxbe calculated as follows:

${\mathrm{\&xi;}}_{2\max }=\underset{\mathrm{\&beta;}\&Element;[0,b],\mathrm{\&theta;}\&Element;\left[\mathrm{0,2}\mathrm{\&pi;}\right]}{\max }\left|{\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})\right|;$

Wherein, ${\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{\mathrm{r\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m[{\sin }^2\left(\frac{\mathrm{m\&beta;}}{2}\cos \mathrm{\&theta;}\right)+\frac{1}{c^2}{\sin }^2\left(\frac{\mathrm{mc\&beta;}}{2}\sin \mathrm{\&theta;}\right)]}-1.$

10. method as claimed in claim 9, it is characterized in that, described the optimization finite difference coefficient obtained is brought in the wave equation of difference scheme, carry out Wave equation forward modeling, also comprise: the optimization finite difference coefficient of the rectangular node obtained is substituted in the two-dimentional Acoustic Wave-equation of difference scheme, carries out two-dimentional Acoustic Wave-equation forward simulation; The two-dimentional Acoustic Wave-equation of described difference scheme is:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,z}^{t-1}-{2P}_{x,z}^t+P_{x,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,z-m}^t+P_{x,z+m}^t)]\\ +\frac{1}{h^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,z}^t+P_{x+m,z}^t)]\end{array};$

Wherein, P is acoustic pressure.