CN108279437A - Variable density ACOUSTIC WAVE EQUATION time higher order accuracy staggering mesh finite-difference method - Google Patents

Abstract

The invention discloses a time high-order precision staggered grid finite difference method for a variable density acoustic wave equation. Firstly, a new difference template is designed; the new difference template is used to carry out staggered grid finite difference numerical discretization of the variable density acoustic wave equation; the dispersion relation of the variable density acoustic wave equation under the new difference template is derived; Taylor series expansion or optimization method is used to obtain Take new finite difference coefficients; enter the new difference coefficients into the corresponding discrete equations to realize the numerical simulation of the variable density acoustic wave equation; analyze the accuracy and stability of the new difference method. The beneficial effect of the present invention is that it has accurate space (2M) order and time 10-order simulation precision. Computational efficiency can be improved by increasing the time step size. In the long-distance simulation, it can better guarantee that the waveform will not be distorted.

Description

Staggered-grid finite difference method for high-order temporal accuracy of variable-density acoustic wave equations

technical field

The invention belongs to the technical field of geophysics, and relates to a new wave equation numerical simulation method, which can be applied to geophysical problems such as seismic wave forward modeling, reverse time migration and full waveform inversion.

Background technique

Numerical simulation of the wave equation is an effective means to understand the propagation of seismic waves in the subsurface and help interpret observation data, and it is also the basic unit of reverse time migration imaging and full waveform inversion. So far, wave equation numerical simulation has been widely used in various links in the field of geophysical exploration. The commonly used numerical solutions of partial differential equations include finite element method, finite difference method, pseudospectral method and so on. The finite difference method is easy to implement and has high computational efficiency, and is the most commonly used wave equation numerical simulation algorithm. Geophysicists have done a lot of research on the finite difference method, and successively proposed explicit finite difference and implicit finite difference; regular grid finite difference, staggered grid finite difference and rotated staggered grid finite difference; based on Taylor series expansion Finite difference and optimization-based finite difference; spatial domain finite difference and time-space domain finite difference, etc. Most of the above-mentioned finite difference methods are designed based on the conventional "cross" type difference template, and the time accuracy is always only second order. When the simulation time is long or the time step is large, these difference methods will produce strong time dispersion, and the simulation accuracy is not enough.

Compared with the conventional method, the time high-order finite difference method can better suppress the time dispersion and has better stability. The Lax-Wendroff method is the earliest time high-order finite difference method, which replaces high-order time derivatives with spatial derivatives to improve time accuracy. However, the Lax-Wendroff method is difficult to extend to any even-order precision, and the calculation amount is large. The study found that the time simulation accuracy of the finite difference method is closely related to the type of difference template used. Subsequently, a finite-difference method based on a "rhombic" template for a regular grid of acoustic wave equations appeared, which can simultaneously obtain any even-order space and time precision. Geophysicists also combine conventional "cross-shaped" templates and "diamond-shaped" templates to obtain space (2M) order and time (2N) order accuracy. In order to further solve the variable-density acoustic wave equation, the researchers proposed a time-higher-order staggered grid finite difference method, and extended it to the case of rectangular or cubic grids. However, most of these methods only have fourth-order or sixth-order time precision.

So far, there are relatively few time high-order staggered grid finite difference methods for variable density acoustic wave equations, and the time accuracy is still not enough. When the time step is large or the simulation distance is long, the numerical dispersion is serious, which has a great impact on the quality of subsequent imaging and inversion.

Contents of the invention

The purpose of the present invention is to provide a time high-order precision staggered grid finite difference method for variable density acoustic wave equations.

The technical scheme adopted in the present invention is to carry out according to the following steps:

(1) Design a new differential template;

(2) Using a new difference template to numerically discretize the variable density acoustic wave equation with staggered grid finite difference;

(3) Deduce the dispersion relation of the variable density acoustic wave equation under the new difference template;

(4) Using Taylor series expansion or optimization method to obtain new finite difference coefficients;

(5) Enter the new differential coefficient into the corresponding discrete equation to realize the numerical simulation of the variable density acoustic wave equation;

(6) Analyze the precision and stability of the new difference method.

Further, the method of designing a new differential template in step (1) is as follows:

Establish the two-dimensional variable density acoustic wave equation:

Where: K=ρv ² is bulk modulus, v is velocity, ρ is density, p is pressure, v=[v _x , v _z ] ^T is polarization velocity vector, ▽ and represent the gradient and divergence operators respectively; the new difference template used in equation (1) calculates the spatial partial derivative by adding grid points on both sides of the coordinate axis, M and N in the new difference template are the operator lengths, d _i,j is the difference coefficient, M controls the number of grid points on the coordinate axis, and N controls the number of grid points outside the coordinate axis.

Further, in step (2), the new difference template is used to carry out the staggered grid finite difference numerical discretization method for the variable density acoustic wave equation as follows:

Among them: Δt is the time sampling interval, h space sampling interval, d _i,j is the differential coefficient, M and N are the operator length, known from the symmetry d _m,n =d _n,m , when N is even and odd , there are only N ² /4 and (N ² -1)/4 independent differential coefficients respectively.

Further, in step (3), the method for deriving the dispersion relation of the variable density acoustic wave equation under the new differential template is as follows:

Under the assumption of plane wave theory, the wave field is expressed as:

in: k is the wave number, θ is the propagation angle, k _x ＝kcos(θ), k _z ＝ksin(θ), ω is the angular frequency;

Substitute equation (3) into (2) and simplify to get:

Where: r=vΔt/h,

Equation (4) is the time-space dispersion relationship of the variable-density acoustic wave equation under the new differential template.

Further, the method of obtaining new finite difference coefficients by using Taylor series expansion or optimization method is as follows:

Expand the sine and cosine function TE in equation (4):

in:

With ωΔt=khr and Simplify equation (6) to get:

E+F=G, (8)

and

Comparing the coefficients of (k _x h) ^2l and (k _z h) ^2l (1≤l≤M) at both ends of equation (8), we get:

Equation (11) is equivalent to:

Comparing both sides of equation (8) The coefficients get:

Equation (12) and equation (13) determine the space and time precision of new differential method respectively, when N=2-5 is provided below, the equivalent equation of equation (13);

(a) N=2, equation (13) becomes

From equation (11), it can be seen that:

Comparing equations (14) and (15) gives:

b _1,1 = a ₂ (16)

(b) N=3, equation (13) becomes

Substitute Equation (16) into Equation (17) to get:

From equation (11) we know

Comparing equations (18) and (19) gives:

b _1,2 +b _2,1 = 3a ₃ (20)

(c) N=4, equation (13) is equivalent to:

Substitute equations (16) and (20) into equations (21) and (22) to get:

From equation (11) we know

Compare (23)–(25) to get:

b _1,3 +b _3,1 = 4a ₄ , b _2,2 = 3a ₄ (26)

(d) N=5, equation (16) is equivalent to:

Equations (16) and (26) are equal to equations (27) and (28):

From equation (11) we know

Comparing equations (29)–(31), we get:

b _1,4 +b _4,1 =5a ₅ , b _2,3 +b _3,2 =10a ₅ , (32) and

In order to satisfy equations (20) and (33) at the same time, the following additional conditions are required:

b _1,2 = a ₃ (34)

When N=5, to make equation (13) valid, it needs to satisfy:

b _1,1 =a ₂ ,b _1,2 +b _2,1 =3a ₃ ,b _1,3 +b _3,1 =4a ₄ ,b _2,2 =3a ₄ ,

b _1,4 +b _4,1 =5a ₅ , b _2,3 +b _3,2 =10a ₅ , b _1,2 =a ₃ . (35)

The beneficial effect of the present invention is that it has accurate space (2M) order and time 10-order simulation precision. Computational efficiency can be improved by increasing the time step size. In the long-distance simulation, it can better guarantee that the waveform will not be distorted.

Description of drawings

Fig. 1 Flowchart of time high-order precision finite difference numerical simulation of variable density acoustic wave equation;

Fig. 2 Temporal high-order precision staggered grid finite difference template;

Figure 3 New staggered grid finite difference template;

Fig. 4v _FD /v variation curve with N;

Fig. 5 The change curve of stability factor s with operator length M and N;

Figure 6. Wave field snapshots of different difference methods for the homogeneous medium model at 1.5s;

Figure 7 Marmousi model;

Fig. 8 Surface seismic records of different difference methods of Marmousi model.

Detailed ways

The present invention will be described in detail below in combination with specific embodiments.

As shown in Figure 1, it is a flow chart of implementing the variable density acoustic wave equation high-precision finite difference numerical simulation for the present invention, specifically including:

(1) Design a new differential template.

Based on the existing regular grid finite difference template, combined with the "cross" and "diamond" templates, a new staggered grid finite difference template is designed.

(2) Using a new difference template to numerically discretize the variable density acoustic wave equation.

The staggered grid finite-difference numerical discretization of the variable-density acoustic wave equation is carried out using a new difference template.

(3) Deduce the dispersion relation of the variable density acoustic wave equation under the new difference template.

Based on the plane wave theory, the time-space dispersion relationship of the variable-density acoustic wave equation under the new differential template is derived.

(4) Calculate the finite difference coefficient.

Use Taylor series expansion or optimization method to find new finite difference coefficients.

(5) Solve the variable density acoustic wave equation numerically.

Put the new difference coefficient into the corresponding discrete equation to realize the numerical simulation of the variable density acoustic wave equation.

(6) Analyze the precision and stability of the new difference method.

Compared with the existing differential method, the simulation accuracy and stability of the new differential method are analyzed.

Specifically, the method of designing a new finite difference template in step 1 is as follows:

(1) Establish a two-dimensional variable density acoustic wave equation:

Where: K=ρv ² is the bulk modulus, v is the velocity, ρ is the density, p is the pressure, v=[v _x ,v _z ] ^T is the polarization velocity vector. ▽ and denote the gradient and divergence operators, respectively.

Equation (1) can be solved by the conventional staggered grid finite difference method, but this method is based on the "cross" type difference template, and the time accuracy is only second order. The present invention adopts a new difference template, as shown in Fig. 2, the temporal high-order precision staggered grid finite difference templates M and N are operator lengths, d _{i, j} are difference coefficients. The new template is an extension of the existing "cross" type difference template, and the spatial partial derivative is calculated by adding grid points on both sides of the coordinate axis. In the new template, M and N are operator lengths, d _{i, j} are differential coefficients, M controls the number of grid points on the coordinate axis, and N controls the number of grid points outside the coordinate axis. When N=1, the grid points on both sides of the coordinate axis disappear, and the new template degenerates into a conventional "cross" type difference template.

In step 2, the method of numerically discretizing the variable-density acoustic wave equation using a new differential template is as follows:

(2) the staggered grid finite difference template shown in Figure 2 is adopted in the present invention to discretize the variable density acoustic wave equation:

Among them: Δt is the time sampling interval, h is the space sampling interval, d _{i, j} are the difference coefficients, M and N are the operator lengths. From the symmetry (d _m,n =d _n,m ), it can be seen that when N is even and odd, there are only N ² /4 and (N ² -1)/4 independent differential coefficients respectively.

In step 3, the method of deriving the dispersion relation of the variable density acoustic wave equation under the new differential template is as follows:

(3) Under the assumption of plane wave theory, the wave field can be expressed as:

in: k is the wave number, θ is the propagation angle (the included angle with the positive semi-axis of the x-axis), k _x =kcos(θ), k _z =ksin(θ), and ω is the angular frequency.

Substitute equation (3) into (2) and simplify to get:

Where: r=vΔt/h,

Equation (4) is the time-space dispersion relationship of the variable-density acoustic wave equation under the new differential template.

The method of obtaining the finite difference coefficient in step 4 is as follows:

(4) Next, use the Taylor series expansion (TE) method to find the difference coefficient. Expand the sine and cosine function TE in equation (4):

in:

With ωΔt=khr and Simplify equation (6) to get:

E+F=G, (8)

and

Comparing the coefficients of (k _x h) ^2l and (k _z h) ^2l (1≤l≤M) at both ends of equation (8), we get:

Equation (11) is equivalent to:

Comparing both sides of equation (8) The coefficients get:

Equation (12) and equation (13) determine the spatial and temporal precision of the new difference method, respectively. As long as the obtained finite difference coefficients can guarantee the establishment of equation (13), the difference method can obtain (2N) order time accuracy. The equivalent equation of equation (13) is given below when N=2-5.

(a) N=2, equation (13) becomes

From equation (11), it can be seen that:

Comparing equations (14) and (15) gives:

b _1,1 = a ₂ . (16)

(b) N=3, equation (13) becomes

Substitute Equation (16) into Equation (17) to get:

From equation (11) we know

Comparing equations (18) and (19) gives:

b _1,2 +b _2,1 ＝3a ₃ . (20)

(c) N=4, equation (13) is equivalent to:

Substitute equations (16) and (20) into equations (21) and (22) to get:

From equation (11) we know

Compare (23)–(25) to get:

b _1,3 +b _3,1 = 4a ₄ , b _2,2 = 3a ₄ . (26)

(d) N=5, equation (16) is equivalent to:

Equations (16) and (26) are equal to equations (27) and (28) to get:

From equation (11) we know

Comparing equations (29)–(31), we get:

b _1,4 +b _4,1 =5a ₅ ,b _2,3 +b _3,2 =10a ₅ , (32)

and

In order to satisfy equations (20) and (33) at the same time, the following additional conditions are required:

b _1,2 = a ₃ . (34)

When N=5, to make equation (13) valid, it needs to satisfy:

b _1,1 =a ₂ ,b _1,2 +b _2,1 =3a ₃ ,b _1,3 +b _3,1 =4a ₄ ,b _2,2 =3a ₄ ,

b _1,4 +b _4,1 =5a ₅ ,b _2,3 +b _3,2 =10a ₅ ,b _1,2 =a ₃ . (35)

In order to obtain space (2M) order and time 10th order accuracy, the differential coefficients need to make M+7 linear equations (Equations 12 and 35) hold simultaneously. However, when N=5, the number of differential coefficients is only M+6. Therefore, the problem is an overdetermined problem, that is, there is no set of solutions (difference coefficients) that satisfy the condition.

A new differential template (as shown in FIG. 3 ) is given below. Comparing Figure 2 and Figure 3, it is found that the new template is very close to that when N=5 in Figure 2, except that there are eight more points (red). The number of differential coefficients in the new template is also changed to M+7 (one more d _5,1 ). The new template is also equivalent to the case of d _4,2 =0, d _3,3 =0 when N=6 in FIG. 2 . In order to avoid re-deriving the dispersion relationship, we set d _4,2 =0, d _3,3 =0 when N=6, and solve Equation 12 and Equation 35 to obtain the difference coefficients. At this time, there are M+7 equations and differential coefficients, which is a well-posed problem. From the above derivation, it can be seen that the newly obtained differential coefficients can ensure that Equation (12) and Equation (13) (N≤5) hold. Therefore, the new difference method has a spatial (2M) order, temporal 10th order simulation accuracy.

In step 5, the numerical solution to the variable density acoustic wave equation is as follows:

(5) Substituting the obtained differential coefficients (d _{i, j} ) into equation (2) for recursion, the time high-order precision staggered grid finite difference numerical simulation of the variable density acoustic wave equation is realized. Based on the simulation accuracy and stability requirements, select appropriate parameters such as spatial sampling interval, time step, and main frequency of the source. In order to avoid the influence of artificially truncated boundaries, the present invention adopts perfectly matched layer (PML) absorption conditions to eliminate boundary reflections.

The method of analyzing the accuracy and stability of the new differential method in step 6 is as follows:

(6) The dispersion formula and the stability formula are respectively:

and

Among them: δ is the relative phase velocity, the closer it is to 1, the smaller the dispersion; s is the stability factor, the larger the value, the better the stability. The accuracy and stability of the new difference method can be analyzed using formulas (36) and (37).

The present invention is a new finite-difference numerical solution method of seismic wave equation, which can simultaneously obtain the simulation accuracy of space (2M) order and time 10th order along any propagation direction; it can effectively suppress the numerical dispersion error, and when the time step is large or the simulation Accurate waveforms can still be obtained at a long distance; it can provide a reliable forward modeling engine for variable density acoustic wave equation reverse time migration and full waveform inversion, thereby greatly improving the accuracy of seismic wave imaging and inversion.

First, a homogeneous medium model is taken as an example to illustrate the advantages of the present invention. Fig. 4 is the dispersion curve when the operator length N=1-5, M=16, r=0.45. When N=1, it is equivalent to the conventional staggered grid finite difference method. The variation curve of v _FD /v with N, M=16, r=0.45, θ=π/4. It can be seen from the figure that the simulation accuracy of the newly proposed finite difference method is higher than that of the conventional staggered grid finite difference method. And with the increase of N, the time dispersion decreases gradually. Figure 5 is the stability curve of N=1-5, r=0.45, the change curve of the stability factor s with the operator length M and N, it can be seen from the figure that the new time high-order finite difference method is more stable than the conventional method it is good. And the smaller N is, the stricter the stability condition is. Figure 6 is a 1.5s wavefield snapshot of different differential methods. Wave field snapshots of different difference methods for homogeneous medium model at 1.5s, M=8, h=10m, Δt=1.5ms. (a) conventional finite difference method (b) new difference method, N=2 (c) new difference method, N=5. The source wavelet is the 20Hz Lake wavelet, which is located in the center of the model. It can be found that the time high-order precision finite difference method can effectively suppress the numerical dispersion. The temporal 10-order precision method (N=5) proposed in the present invention has higher simulation precision.

Next, the complex Marmousi model (as shown in Fig. 7) is used to test the newly proposed differential method. Fig. 8 shows the ground seismic records of different difference methods of Marmousi model, M=8, h=10m, Δt=1.0ms. (a) Conventional finite difference method. (b) New difference method, N=2. (c) New difference method, N=5. The source wavelet is a 15Hz Reckel wavelet, which is located in the middle of the earth's surface. It can be seen from the figure that the time high-order precision staggered grid finite difference method proposed in the present invention can well suppress numerical dispersion, and can be applied to forward modeling of complex media middle.

The above description is only a preferred embodiment of the present invention, and does not limit the present invention in any form. Any simple modifications made to the above embodiments according to the technical essence of the present invention, equivalent changes and modifications, all belong to this invention. within the scope of the technical solution of the invention.

Claims (5)

1. The time high-order precision staggered grid finite difference method of the variable density acoustic wave equation is characterized in that it is carried out according to the following steps: (1) Design a new differential template; (2) Using a new difference template to numerically discretize the variable density acoustic wave equation with staggered grid finite difference; (3) Deduce the dispersion relation of the variable density acoustic wave equation under the new difference template; (4) Using Taylor series expansion or optimization method to obtain new finite difference coefficients; (5) Enter the new differential coefficient into the corresponding discrete equation to realize the numerical simulation of the variable density acoustic wave equation; (6) Analyze the precision and stability of the new difference method. 2. according to the described variable density acoustic wave equation time high-order precision staggered grid finite difference method of claim 1, it is characterized in that: in the described step (1), the new difference template method of designing is as follows: Establish the two-dimensional variable density acoustic wave equation: Where: K=ρv ² is the bulk modulus, v is the velocity, ρ is the density, p is the pressure, v=[v _x , v _z ] ^T is the polarization velocity vector, and represent the gradient and divergence operators respectively; the new difference template used in equation (1) calculates the spatial partial derivative by adding grid points on both sides of the coordinate axis, M and N in the new difference template are the operator lengths, d _i,j is the difference coefficient, M controls the number of grid points on the coordinate axis, and N controls the number of grid points outside the coordinate axis. 3. according to the time high-order precision staggered grid finite difference method of the variable density acoustic wave equation of claim 1, it is characterized in that: adopt new differential template to carry out the staggered grid finite difference to variable density acoustic wave equation in the described step (2) The numerical discretization method is as follows: Among them: Δt is the time sampling interval, h space sampling interval, d _i,j is the differential coefficient, M and N are the operator length, known from the symmetry d _m,n =d _n,m , when N is even and odd , there are only N ² /4 and (N ² -1)/4 independent differential coefficients respectively. 4. according to the described variable density acoustic wave equation time high-order precision staggered grid finite difference method of claim 1, it is characterized in that: in the described step (3), the dispersion relation method of variable density acoustic wave equation under deriving new differential template is as follows: Under the assumption of plane wave theory, the wave field is expressed as: in: k is the wave number, θ is the propagation angle, k _x ＝kcos(θ), k _z ＝ksin(θ), ω is the angular frequency; Substitute equation (3) into (2) and simplify to get: Where: r=vΔt/h, Equation (4) is the time-space dispersion relationship of the variable-density acoustic wave equation under the new differential template. 5. according to the said variable density acoustic wave equation time high-order precision staggered grid finite difference method of claim 4, it is characterized in that: adopt Taylor series expansion or optimization method to seek new finite difference coefficient method as follows: Expand the sine and cosine function TE in equation (4): in: With ωΔt=khr and Simplify equation (6) to get: E+F=G, (8) and Comparing the coefficients of (k _x h) ^2l and (k _z h) ^2l (1≤l≤M) at both ends of equation (8), we get: Equation (11) is equivalent to: Comparing both sides of equation (8) The coefficients get: Equation (12) and equation (13) determine the space and time precision of new differential method respectively, when N=2-5 is provided below, the equivalent equation of equation (13); (a) N=2, equation (13) becomes From equation (11), it can be seen that: Comparing equations (14) and (15) gives: b _1,1 = a ₂ (16) (b) N=3, equation (13) becomes Substitute Equation (16) into Equation (17) to get: From equation (11) we know Comparing equations (18) and (19) gives: b _1,2 +b _2,1 = 3a ₃ (20) (c) N=4, equation (13) is equivalent to: Substitute equations (16) and (20) into equations (21) and (22) to get: From equation (11) we know Compare (23)–(25) to get: b _1,3 +b _3,1 = 4a ₄ , b _2,2 = 3a ₄ (26) (d) N=5, equation (16) is equivalent to: Equations (16) and (26) are equal to equations (27) and (28): From equation (11) we know Comparing equations (29)–(31), we get: b _1,4 +b _4,1 =5a ₅ ,b _2,3 +b _3,2 =10a ₅ , (32) and In order to satisfy equations (20) and (33) at the same time, the following additional conditions are required: b _1,2 = a ₃ (34) When N=5, to make equation (13) valid, it needs to satisfy: b _1,1 =a ₂ ,b _1,2 +b _2,1 =3a ₃ ,b _1,3 +b _3,1 =4a ₄ ,b _2,2 =3a ₄ , b _1,4 +b _4,1 =5a ₅ , b _2,3 +b _3,2 =10a ₅ , b _1,2 =a ₃ . (35)