Optimized Finite Difference Methods for Seismic Acoustic Wave Modeling

Abstract

The finite difference (FD) methods are widely used for approximating the partial derivatives in the acoustic/elastic wave equation. Grid dispersion is one of the key numerical problems and will directly influence the accuracy of the result because of the discretization of the partial derivatives in the wave equation. Therefore, it is of great importance to suppress the grid dispersion by optimizing the FD coefficient. Various optimized methods are introduced in this chapter to determine the FD coefficient. Usually, the identical staggered grid finite difference operator is used for all of the first-order spatial derivatives in the first-order wave equation. In this chapter, we introduce a new staggered grid FD scheme which can improve the efficiency while still preserving high accuracy for the first-order acoustic/elastic wave equation modeling. It uses different staggered grid FD operators for different spatial derivatives in the first-order wave equation. The staggered grid FD coefficients of the new FD scheme can be obtained with a linear method. At last, numerical experiments were done to demonstrate the effectiveness of the introduced method.

1. Introduction

The propagation of seismic waves through the Earth’s subsurface is described by the wave equation, one of the partial differential equations (PDEs), which describe many of the fundamental natural laws. When the subsurface earth structure is complex, it is difficult to obtain the analytic results. The finite difference (FD) method is one of most widely used numerical methods for wave equation modeling because of its high efficiency, smaller memory requirement, and easy implementation [1–7].

The first application of the FD method to wave equation modeling can be possibly traced back to Alterman and Karal [1]. Alford et al. took the grid dispersion analysis for the second-order and fourth-order FD operators and stated that it is necessary to use high-order FD operators [8]. Kelly et al. further illustrated the grid dispersion, absorbing boundary condition, and other implementation aspects of the FD method [9].

Madariaga developed a staggered grid FD scheme to solve a rupture propagation problem [10]. Virieux adapted this scheme to elastic SH waves and P-SV waves in a 2D Cartesian system [2, 3]. Levander introduced a fourth-order staggered grid FD operator in the space domain to improve accuracy [11].

Grid dispersion is one of the key numerical problems affecting practical usage when utilizing the FD method. Since the traditional FD coefficient obtained in the space domain with the Taylor expansion method is only accurate for a very limited wavenumber range [4], many efforts are paid to reducing the grid dispersion with optimized FD coefficient. Yang et al. proposed the nearly analytic discrete method for wave equation and later improved this method [12, 13]. Chen proposed high-order time discretization method to reduce the dispersion caused by the temporal discretization [14, 15]. The Fourier FD was introduced by Song and Fomel with the combination of fast Fourier transform and finite difference operators [16]. Chu and Stoffa improved the FD methods with a scaled binomial windowed FD scheme that leads to more precise discrete operators [17]. Fomel et al. introduced low-rank approximation of the wave propagator matrix to reduce the cost of wave extrapolation [18].

Generally, the FD coefficients of the spatial derivative are determined only in the spatial domain. However, wave equations are solved in the temporal and spatial domains simultaneously. Finkelstein and Kastner propose a systematic design methodology for obtaining FD coefficients to reduce dispersion, which allows the exact phase velocity or (and) group velocity dispersion relationship to be satisfied at some designated frequencies in the temporal-spatial domain [19, 20]. Etgen proposed minimizing the phase velocity error using the least squares (LS) method [21]. Liu and Sen propose a new time-space domain method to determine the higher order FD coefficients for 1D, 2D, and 3D wave equations [22], and then they use this method to get the staggered grid FD coefficients [23]. Zhang and Yao proposed the use of the simulated annealing algorithm and gave an error limitation for determining the FD coefficients in the space or the time-space domain [24]. Liang et al. proposed utilizing the linear method to determine the FD coefficient in the time-space domain [25]. Ren and Liu developed a novel optimal time-space domain staggered grid FD scheme and used least squares method to get the FD coefficients [26]. Wang et al. proposed the regularized optimization method to get the staggered grid FD coefficient in the time-space domain [27]. Chen et al. used K space operator-based high-order staggered grid FD method to improve accuracy [28]. Yong et al. proposed using the optimized equivalent staggered grid FD method with three sets of FD coefficients to improve the simulation accuracy [29]. Compared with the traditional high-order staggered grid FD coefficient obtained by the Taylor expansion method, these methods greatly improved the accuracy with the optimized FD coefficient.

Another way to improve the accuracy and efficiency of the FD methods is using new FD stencil. Liu and Sen studied the rhombus stencil and found that it can reach high-order accuracy along all directions [30]. Liu et al. formulated an explicit time evolution scheme with high temporal accuracy by using a new FD stencil for the second-order wave equation [31]. Tan and Huang propose a staggered grid FD stencil with added points in the diagonal direction for the first-order wave equation [32, 33]. Compared with the traditional staggered grid FD stencil, these methods improved the efficiency by using a larger time step while still preserving high accuracy.

2. Acoustic/elastic wave equations

The first-order velocity-stress acoustic wave equation can be described as

$$\frac{\partial P}{\partial t}=-v^2\left(\frac{\partial v_z}{\partial z}+\frac{\partial v_x}{\partial x}\right),$$ E1

$$\frac{\partial v_x}{\partial t}=-\frac{\partial P}{\partial x},$$ E2

$$\frac{\partial v_z}{\partial t}=-\frac{\partial P}{\partial z}.$$ E3

where P is the acoustic pressure fluctuation, v is the wave propagation speed, and v_x and v_z are the particle velocities.

Substituting Eqs. (2) and (3) into Eq. (1), the second-order acoustic wave equation can be written as

$$\frac{{\partial }^{2}p}{\partial x^2}+\frac{{\partial }^{2}p}{\partial z^2}=\frac{1}{v^2}\frac{{\partial }^{2}p}{\partial t^2}$$ E4

The first-order elastic wave equations in 2D heterogeneous media are [3]

$$\frac{\partial v_x}{\partial t}=\frac{\partial {\tau }_{\mathit{xx}}}{\partial x}+\frac{\partial {\tau }_{\mathit{xz}}}{\partial z},$$ E5

$$\frac{\partial v_z}{\partial t}=\frac{\partial {\tau }_{\mathit{xz}}}{\partial x}+\frac{\partial {\tau }_{\mathit{zz}}}{\partial z},$$ E6

$$\frac{\partial {\tau }_{\mathit{xx}}}{\partial t}={\alpha }^2\frac{\partial v_x}{\partial x}+\left({\alpha }^2-2{\beta }^2\right)\frac{\partial v_z}{\partial z},$$ E7

$$\frac{\partial {\tau }_{\mathit{zz}}}{\partial t}={\alpha }^2\frac{\partial v_z}{\partial z}+\left({\alpha }^2-2{\beta }^2\right)\frac{\partial v_x}{\partial x},$$ E8

$$\frac{\partial {\tau }_{\mathit{xz}}}{\partial t}={\beta }^2\left(\frac{\partial v_x}{\partial z}+\frac{\partial v_z}{\partial x}\right).$$ E9

where (v_x , v_z ) is the velocity vector,(τ_xx , τ_zz , τ_xz ) is the stress vector, and α and β are the P- and S wave propagation speeds, respectively.

Substituting Eqs. (7)–(9) into Eqs. (5)–(6), the second-order elastic wave equation can be written as

$$\frac{{\partial }^2{v}_x}{\partial t^2}={\alpha }^2\frac{{\partial }^2{v}_x}{\partial x^2}+\left({\alpha }^2-{\beta }^2\right)\frac{{\partial }^2{v}_z}{\partial x\partial z}+{\beta }^2\frac{{\partial }^2{v}_x}{\partial z^2},$$ E10

$$\frac{{\partial }^2{v}_z}{\partial t^2}={\alpha }^2\frac{{\partial }^2{v}_z}{\partial z^2}+\left({\alpha }^2-{\beta }^2\right)\frac{{\partial }^2{v}_x}{\partial x\partial z}+{\beta }^2\frac{{\partial }^2{v}_z}{\partial x^2}.$$ E11

3. Finite difference operators

The commonly used staggered grid FD scheme for the first-order acoustic wave equation is as follows:

$$\frac{\partial P}{\partial t}=\frac{-v^2}{h}\left(\sum _{m=1}^{M_1}{c}_m\left(v_{z_{0,m-1/2}}-v_{z_{0,-m+1/2}}\right)+\sum _{m=1}^{M_1}{c}_m\left(v_{x_{m-1/2,0}}-v_{x_{-m+1/2,0}}\right)\right),$$ E12

$$\frac{\partial v_x}{\partial t}=\frac{-1}{h}\sum _{m=1}^{M_2}{c}_m\left(P_{m-1/2,0}^0-P_{-m+1/2,0}^0\right),$$ E13

$$\frac{\partial v_z}{\partial t}=\frac{-1}{h}\sum _{m=1}^{M_2}{c}_m\left(P_{0,m-1/2}^0-P_{0,-m+1/2}^0\right),$$ E14

$$Q_{m,j}^n=Q\left(x+\mathit{mh},z+\mathit{jh},t+\mathit{n\tau }\right);Q=v_x,v_z,P$$ E15

where M ₁ and M ₂ are the length of the FD operators, c_m is the staggered grid FD coefficients, and h is the spatial grid interval.

The second-order FD operator is usually used for the first-order time derivative:

$$\frac{\partial P}{\partial t}=\frac{1}{\Delta t}\left(P_{0,0}^1-P_{0,0}^0\right).$$ E16

where Δt is the time step.

The commonly used staggered grid FD scheme for the first-order elastic wave equation is as follows:

$$\frac{\partial v_x}{\partial t}=\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{xx}{}_{{}_{m-1/2,0}}^0-{\tau }_{xx}{}_{{}_{-m+1/2,0}}^0\right]}+\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{xz}{}_{{}_{0,m-1/2}}^0-{\tau }_{xz}{}_{{}_{0,-m+1/2}}^0\right]}$$ E17

$$\frac{\partial v_z}{\partial t}=\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{xz}{}_{{}_{m-1/2,0}}^0-{\tau }_{xz}{}_{{}_{-m+1/2,0}}^0\right]}+\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{zz}{}_{{}_{0,m-1/2}}^0-{\tau }_{zz}{}_{{}_{0,-m+1/2}}^0\right]}$$ E18

$$\frac{\partial {\tau }_{xx}}{\partial t}=\frac{{\alpha }^2}{h}{\displaystyle \sum _{m=1}^{M_2}{c}_m\left[v_x{}_{{}_{m-1/2,0}}^0-v_x{}_{{}_{-m+1/2,0}}^0\right]}+\frac{{\alpha }^2-2{\beta }^2}{h}{\displaystyle \sum _{m=1}^{M_2}{c}_m\left[v_z{}_{{}_{0,m-1/2}}^0-v_z{}_{{}_{0,-m+1/2}}^0\right]}$$ E19

$$\frac{\partial {\tau }_{zz}}{\partial t}=\frac{{\alpha }^2-2{\beta }^2}{h}{\displaystyle \sum _{m=1}^{M_2}{c}_m\left[v_x{}_{{}_{m-1/2,0}}^0-v_x{}_{{}_{-m+1/2,0}}^0\right]}+\frac{{\alpha }^2}{h}{\displaystyle \sum _{m=1}^{M_2}{c}_m\left[v_z{}_{{}_{0,m-1/2}}^0-v_z{}_{{}_{0,-m+1/2}}^0\right]}$$ E20

$$\frac{\partial {\tau }_{xz}}{\partial t}=\frac{{\beta }^2}{h}{\displaystyle \sum _{m=1}^{M_2}{c}_m\left[v_z{}_{{}_{m-1/2,0}}^0-v_z{}_{{}_{-m+1/2,0}}^0\right]}+\frac{{\beta }^2}{h}{\displaystyle \sum _{m=1}^{M_2}{c}_m\left[v_x{}_{{}_{0,m-1/2}}^0-v_x{}_{{}_{0,-m+1/2}}^0\right]}$$ E21

$$Q_{m,j}^n=Q\left(x+\mathit{mh},z+\mathit{jh},t+\mathit{n\tau }\right);Q=v_x,v_z,{\tau }_{\mathit{xx}},{\tau }_{\mathit{zz}},{\tau }_{\mathit{xz}}$$ E22

where M ₁ and M ₂ are the length of the FD operators, c_m is the staggered grid FD coefficient to be determined, and h is the spatial grid interval.

4. Optimizing finite difference operators

4.1. Optimizing finite difference operators for the acoustic wave equation

Using the plane wave theory, let

$$P_{m,n}^j=e^{i\left[k_x\left(x+\mathit{mh}\right)+k_z\left(z+\mathit{jh}\right)-\omega \left(t+\mathit{n\tau }\right)\right]}.$$ E23

The following dispersion relation can be obtained by substituting Eqs. (13)–(14) into Eq. (12) [23, 26, 27]:

$${\left[\sum _{m=1}^M{c}_{m}sin\left(\left(m-0.5\right)k_{x}h\right)\right]}^2+{\left[\sum _{m=1}^M{c}_{m}sin\left(\left(m-0.5\right)k_{z}h\right)\right]}^2\approx \frac{1}{2r^2}\left[1-cos\left(\mathit{kv\tau }\right)\right].$$ E24

where r = vΔt/h, M₁ = M₂ = M, and(k_x , k_z ) = k(cosθ, sinθ). It can be observed from Eq. (24) that the dispersion relation is complex and optimized methods are needed to address this problem.

Let c be the vector form of the FD coefficients, and denote the left side of Eq. 24 by [27]:

$$F\left(c\right)=\sum _{\theta =0}^{2\pi }\left\{{\left[\sum _{m=1}^M{c}_{m}sin\left(\left(m-0.5\right)k_{x}h\right)\right]}^2+{\left[\sum _{m=1}^M{c}_{m}sin\left(\left(m-0.5\right)k_{z}h\right)\right]}^2\right\}$$ E25

and the right side of Eq. 24 by

$$d=\sum _{\theta =0}^{2\pi }\frac{1}{2r^2}\left[1-cos\left(\mathit{kv\tau }\right)\right].$$ E26

The aim is to minimize the dispersion error for a fixed range of wavenumbers:

$$\mathrm{\Phi }\left(c\right)=\sum _{k=0}^{{\displaystyle \stackrel{\sim }{K}}}{\left[F\left(c\right)-d\right]}^2\to min$$ E27

The upper limit of the wavenumber range used for calculating the FD coefficients is based on the source frequency, the space grid interval, and the wave propagation speed [25]:

$$\mathit{\text{Ratio}}=\frac{{\displaystyle \stackrel{\sim }{K}}}{K_{\mathit{\text{total}}}}=\frac{2\pi f_{\max }/v}{\pi /h}=\frac{f_{\max }}{v/2h}.$$ E28

The direct minimization of the objective function Φ for the FD coefficient may lead to unstable results. Therefore, regularizing technique was applied to restore stability. The regularization model is established as

$$J^{\alpha }\left(c\right)=\Phi \left(c\right)+\frac{1}{2}\alpha {‖\mathit{Dc}‖}^2,$$ E29

where α > 0 is a user-defined regularization parameter and D is a scale operator. The new task is the minimization of J(c), and then the regularized optimization staggered grid FD coefficient can be obtained.

Another way to improve the efficiency and accuracy of the staggered grid FD methods is the utilization of the new staggered grid FD scheme. Different with the previous staggered grid FD scheme, the simplest centered second-order staggered grid FD operator can be used for the spatial derivatives in Eqs. (2) and (3), for example,

$$\frac{\partial P}{\partial t}=\frac{-v^2}{h}\left(\sum _{m=1}^M{c}_m\left(v_{z_{0,m-1/2}}-v_{z_{0,-m+1/2}}\right)+\sum _{m=1}^M{c}_m\left(v_{x_{m-1/2,0}}-v_{x_{-m+1/2,0}}\right)\right),$$ E30

$$\frac{\partial v_x}{\partial t}=\frac{-\left(P_{1/2,0}^0-P_{-1/2,0}^0\right)}{h},$$ E31

$$\frac{\partial v_z}{\partial t}=\frac{-\left(P_{0,1/2}^0-P_{0,-1/2}^0\right)}{h}.$$ E32

The staggered grid FD scheme in Eqs. (30)–(32) can be seen as a new staggered grid FD scheme for the first-order acoustic wave equation. The new staggered grid FD scheme is exactly the same as the traditional staggered grid FD scheme except if the staggered grid FD operator length is shorter for Eqs. (31) and (32). By carefully comparing Eqs. (12) and (14) with Eqs. (30) and (32), we find that the new staggered grid FD scheme is more efficient and can save about 45% of simulation time when M equals 7. It looks like the particle velocities v_x and v_z in Eqs. (31) and (32) are inaccurate since only the second-order staggered grid FD operators are used. However, this is not true since the staggered grid FD coefficient in Eq. (30) is optimized with Eqs. (31) and (32) in consideration. In the following, the huge advantage of the new staggered grid FD scheme will be demonstrated because it can reduce the simulation time while still preserving high accuracy compared with the traditional staggered grid FD scheme.

To get the staggered grid FD coefficient in Eq. (30), we substitute Eqs. (31) and (32) into Eq. (30), using the plane wave theory. Then we get

$$\begin{array}{l}\left(e^{i{k}_z\frac{h}{2}}-e^{-i{k}_z\frac{h}{2}}\right){\displaystyle \sum }_{m=1}^M{c}_m\left(e^{i{k}_z\left(m-\frac{1}{2}\right)h}-e^{-i{k}_z\left(m-\frac{1}{2}\right)h}\right)\hfill \\ +\left(e^{i{k}_x\frac{h}{2}}-e^{-i{k}_x\frac{h}{2}}\right){\displaystyle \sum }_{m=1}^M{c}_m\left(e^{i{k}_x\left(m-\frac{1}{2}\right)h}-e^{-i{k}_x\left(m-\frac{1}{2}\right)h}\right)=\frac{h^2}{v^2}\frac{e^{i\omega \tau }+e^{-i\omega \tau }-2}{\Delta t^2}\hfill \end{array}$$ E33

From Eq. (33), the following dispersion relation can be obtained in the frequency-wavenumber domain (it is a special case of Eq. (24)):

$$-2sin\frac{k_{z}h}{2}\sum _{m=1}^M{c}_m\left(sin\left(m-0.5\right)k_{z}h\right)-2sin\frac{k_{x}h}{2}\sum _{m=1}^M{c}_m\left(sin\left(m-0.5\right)k_{x}h\right)=r^{-2}\left[cos\left(\mathit{\omega \tau }\right)-1\right].$$ E34

Using the basic trigonometric function

$$sin\alpha sin\beta =\frac{-\left[cos\left(\alpha +\beta \right)-cos\left(\alpha -\beta \right)\right]}{2},$$ E35

we obtain Eq. (36) from Eq. (34):

$$\sum _{m=1}^M{c}_m\left[cos\left({\mathit{mk}}_{x}h\right)-cos\left(\left(m-1\right)k_{x}h\right)+cos\left({\mathit{mk}}_{z}h\right)-cos\left(\left(m-1\right)k_{z}h\right)\right]=r^{-2}\left[cos\left(\mathit{\omega \tau }\right)-1\right].$$ E36

Similarly, the new dispersion relation for the 3D first-order acoustic wave equation in the frequency-wavenumber domain is

$$\begin{array}{c}\sum _{m=1}^M{c}_m\left[cos\left({\mathit{mk}}_{x}h\right)-cos\left(\left(m-1\right)k_{x}h\right)+cos\left({\mathit{mk}}_{y}h\right)-cos\left(\left(m-1\right)k_{y}h\right)+cos\left({\mathit{mk}}_{z}h\right)-cos\left(\left(m-1\right)k_{z}h\right)\right]\\ =r^{-2}\left[cos\left(\mathit{\omega \tau }\right)-1\right]\end{array}$$ E37

where $\left(k_x,k_y,k_z\right)=k\left(sin\theta cos\phi ,sin\theta sin\phi ,cos\theta \right).$

Compared with the traditional dispersion relation in Eq. (24), the new dispersion relation in Eqs. (36) and (37) is linear and much simpler.

We assume that there are M equally distributed wavenumber points satisfying the dispersion relation within the wavenumber range specified by Eq. (28). Then, we establish the linear equation from Eq. (37) for the 3D case [25]:

$$\begin{array}{c}\sum _{\phi =0}^{2\pi }\sum _{\theta =0}^{2\pi }\left(\begin{array}{l}{a}_{k_{1,x},1}^h+a_{k_{1,y},1}^h+a_{k_{1,z},1}^h\dots a_{k_{1,x},M}^h+a_{k_{1,y},M}^h+a_{k_{1,x},M}^h\\ ⋮⋮\\ a_{k_{M,x},1}^h+a_{k_{M,y},1}^h+a_{k_{M,z},1}^h\dots a_{k_{M,x},M}^h+a_{k_{M,y},M}^h+a_{k_{M,x},M}^h\end{array}\right)\left(\begin{array}{l}{c}_1\\ ⋮\\ c_M\end{array}\right)\hfill \\ =\frac{1}{r^2}\sum _{\phi =0}^{2\pi }\sum _{\theta =0}^{2\pi }\left(\begin{array}{l}cos\left(k_1\mathit{v\tau }\right)-1\\ ⋮\\ cos\left(k_M\mathit{v\tau }\right)-1\end{array}\right)\hfill \end{array}$$ E38

where $a_{k_l,m}^h=cos\left({\mathit{mk}}_{l}h\right)-cos\left(\left(m-1\right)k_{l}h\right)$ , the ith component of k_l (l = x, y, z) is represented as k _i, l , k_x  = k cos θ cos φ, k_y  = k cos θ sin φ, k_z  = k sin θ, and k(i) for each i = 1,2,…,M + 1 is equally distributed between 0 and Ratio × π / h, where Ratio is determined by Eq. (28). In the following, we will demonstrate that the new staggered grid FD scheme in Eqs. (30)–(32) has similar accuracy compared with the computational intensive traditional staggered grid FD scheme in Eqs. (12)–(14).

The 2D dispersion error δ of the new staggered grid FD scheme is defined as

$$\delta =\frac{v_{\mathit{FD}}}{v}=\frac{1}{\mathit{rkh}}arccos\left(1+r^2{q}_2\right).$$ E39

where

$$q_2=\sum _{m=1}^M{c}_m\left[cos\left({\mathit{mk}}_{x}h\right)-cos\left(\left(m-1\right)k_{x}h\right)+cos\left({\mathit{mk}}_{z}h\right)-cos\left(\left(m-1\right)k_{z}h\right)\right]$$ E40

The difference between the FD propagation time and the exact propagation time through one grid is defined as [23]

$$\epsilon =\frac{h}{v_{\mathit{FD}}}-\frac{h}{v}=\frac{h}{v}\left(\frac{v}{v_{\mathit{FD}}}-1\right)=\frac{h}{v}\left(\frac{1}{\delta }-1\right)$$ E41

Figures 1 and 2 show the dispersion error curves of the traditional and the new staggered grid FD schemes for the homogeneous acoustic model in 2D. All the FD coefficients are determined in the time-space domain with M = 7. From Figures 1 and 2 , we get the conclusion that the new staggered grid FD scheme can also preserve the dispersion relation in a pretty wider range compared with the traditional staggered grid FD methods. For example, with r = 0.0075 in the 2D case, both of them can preserve the dispersion error under 10⁻⁵ within 80% of kh range. However, the new staggered grid FD scheme saves wave equation simulation time because Eqs. (31) and (32) are much simpler than Eqs. (13) and (14).

Figure 1.

Dispersion error curves of the traditional staggered grid FD schemes. (a) r = 0.075 and (b) r = 0.225.

Figure 2.

Dispersion error curves of the new staggered grid FD schemes. (a) r = 0.075 and (b) r = 0.225.

Let the left part of Eq. (36) as.

$$g=\sum _{m=1}^M{c}_m\left[cos\left({\mathit{mk}}_{x}h\right)-cos\left(\left(m-1\right)k_{x}h\right)+cos\left({\mathit{mk}}_{z}h\right)-cos\left(\left(m-1\right)k_{z}h\right)\right].$$ E42

From dispersion relation Eq. (36), it is obvious that.

$$g<0;r^{2}g\ge -2.$$ E43

Then, the stability condition of the new staggered grid FD scheme is (from Eq. (36) with kh = π).

$$r\le \sqrt{\raisebox{1ex}{$-2$}\!\left/ \!\raisebox{-1ex}{$g$}\right.}\le \sqrt{\raisebox{1ex}{$-2$}\!\left/ \!\raisebox{-1ex}{$\sum _{m=1}^{M}4c_m{\left(-1\right)}^m$}\right.}=\frac{1}{\sqrt{2}}{\sqrt{\sum _{m=1}^M{c}_m}}^{-1}.$$ E44

Figure 3 shows the stability condition of the traditional and the new staggered grid FD scheme in 2D. We can see that the stability condition becomes stricter with the increase of the FD operator length. It also shows that the new staggered grid FD scheme’s stability condition is a little bit better than the previous staggered grid FD scheme. For example, the stability conditions are r < 0.54 and r < 0.57, respectively, for the traditional and the new staggered grid FD scheme with M = 7.

Figure 3.

Stability condition in 2D. (a) The traditional staggered grid FD scheme and (b) the new staggered grid FD scheme.

4.2. Optimizing finite difference operators for the elastic wave equation

Eqs. (10) and (11) can be written as [11].

$$\left(\begin{array}{l}{\alpha }^2{D}_{\mathit{xx}}+{\beta }^2{D}_{\mathit{zz}}-D_{\mathit{tt}}\left({\alpha }^2-{\beta }^2\right)D_{\mathit{xz}}\\ \left({\alpha }^2-{\beta }^2\right)D_{\mathit{xz}}{\alpha }^2{D}_{\mathit{zz}}+{\beta }^2{D}_{\mathit{xx}}-D_{\mathit{tt}}\end{array}\right)\left(\begin{array}{l}{v}_x\\ v_z\end{array}\right)=0.$$ E45

The two roots give the following dispersion relation [11]:

$$D_{\mathit{tt}}=\frac{1}{2}\left({\alpha }^2+{\beta }^2\right)\left(D_{\mathit{xx}}+D_{\mathit{zz}}\right)\pm \frac{1}{2}\left({\alpha }^2-{\beta }^2\right)\sqrt{{\left(D_{\mathit{xx}}+D_{\mathit{zz}}\right)}^2-4\left(D_{\mathit{xx}}{D}_{\mathit{zz}}-D_{\mathit{xz}}{D}_{\mathit{xz}}\right)}.$$ E46

Suppose D _xx D_zz = D _xz D _xz , then two equations can be obtained from Eq. (46):

$${\alpha }^2{D}_{\mathit{xx}}+{\alpha }^2{D}_{\mathit{zz}}-D_{\mathit{tt}}=0.$$ E47

$${\beta }^2{D}_{\mathit{xx}}+{\beta }^2{D}_{\mathit{zz}}-D_{\mathit{tt}}=0.$$ E48

Usually, Eq. (48) is used to determine the FD coefficient. For the first-order staggered grid FD scheme, the following dispersion relation can be obtained from Eq. (48):

$${\left[\sum _{m=1}^M{c}_{m}sin\left(\left(m-0.5\right)k_{x}h\right)\right]}^2+{\left[\sum _{m=1}^M{c}_{m}sin\left(\left(m-0.5\right)k_{z}h\right)\right]}^2\approx \frac{1}{2r^2}\left[1-cos\left(\mathit{kv\tau }\right)\right].$$ E49

where r = βΔt/h, M₁ = M₂ = M, and (k_x , k_z ) = k(cosθ, sinθ). It can be observed from Eq. (49) that the dispersion relation is nonlinear and regularized optimized methods can address this problem similarly.

Different with previous staggered grid FD scheme for the first-order elastic wave equation, the simplest centered second-order staggered grid FD operator can be used for the spatial derivatives in Eqs. (7)–(9):

$$\frac{\partial v_x}{\partial t}=\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{xx}{}_{{}_{m-1/2,0}}^0-{\tau }_{xx}{}_{{}_{-m+1/2,0}}^0\right]}+\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{xz}{}_{{}_{0,m-1/2}}^0-{\tau }_{xz}{}_{{}_{0,-m+1/2}}^0\right]}$$ E50

$$\frac{\partial v_z}{\partial t}=\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{xz}{}_{{}_{m-1/2,0}}^0-{\tau }_{xz}{}_{{}_{-m+1/2,0}}^0\right]}+\frac{1}{h}{\displaystyle \sum _{m=1}^{M_1}{c}_m\left[{\tau }_{zz}{}_{{}_{0,m-1/2}}^0-{\tau }_{zz}{}_{{}_{0,-m+1/2}}^0\right]}$$ E51

$$\frac{\partial {\tau }_{xx}}{\partial t}=\frac{{\alpha }^2}{h}\left[v_x{}_{{}_{1/2,0}}^0-v_x{}_{{}_{-1/2,0}}^0\right]+\frac{{\alpha }^2-2{\beta }^2}{h}\left[v_z{}_{{}_{0,1/2}}^0-v_z{}_{{}_{0,-1/2}}^0\right]$$ E52

$$\frac{\partial {\tau }_{xx}}{\partial t}=\frac{{\alpha }^2-2{\beta }^2}{h}\left[v_x{}_{{}_{1/2,0}}^0-v_x{}_{{}_{-1/2,0}}^0\right]+\frac{{\alpha }^2}{h}\left[v_z{}_{{}_{0,1/2}}^0-v_z{}_{{}_{0,-1/2}}^0\right]$$ E53

$$\frac{\partial {\tau }_{xz}}{\partial t}=\frac{{\beta }^2}{h}\left[v_z{}_{{}_{1/2,0}}^0-v_z{}_{{}_{-1/2,0}}^0\right]+\frac{{\beta }^2}{h}\left[v_x{}_{{}_{0,1/2}}^0-v_x{}_{{}_{0,-1/2}}^0\right]$$ E54

The staggered grid FD scheme in Eqs. (52)–(54) is more efficient than the staggered grid FD scheme in Eqs. (19)–(21). It will be demonstrated later that the staggered grid FD scheme in Eqs. (52)–(54) is accurate for the stress vector (τ_xx , τ_zz , τ_xz ) even when only second-order staggered grid FD operator is used.

Then, the new dispersion relation can be obtained from Eq. (49) in the frequency-wavenumber domain:

$$\sum _{m=1}^{M_1}{c}_m\left[cos\left({\mathit{mk}}_{x}h\right)-cos\left(\left(m-1\right)k_{x}h\right)+cos\left({\mathit{mk}}_{z}h\right)-cos\left(\left(m-1\right)k_{z}h\right)\right]=r^{-2}\left[cos\left(\mathit{\omega \tau }\right)-1\right].$$ E55

The staggered grid FD coefficient can be obtained similarly using the linear method.

5. Experiments

5.1. Acoustic wave equation

5.1.1. Numerical modeling in the layered velocity model

We first consider a layered velocity model. The velocity is 1500 m/s for the first layer and 2500 m/s for the second layer as shown in Figure 4 . The sponge boundary code in CREWES Toolbox is used to reduce artificial reflection waves [34]. A Ricker wavelet with the main frequency as 14.3 Hz was used as the seismic source. The seismic source position is denoted as a asterisk, and the receivers A and B are denoted as a circle and a diamond from top to bottom, respectively. The space grid interval is 20m, the FD operator length M is 7 and the time step is 1.5 ms. The staggered grid FD coefficients used in Figure 4 are shown in Table 1 .

Figure 4.

Velocity model.

	c1	c2	c3	c4	c5	c6	c7
	1.22861	−0.102384	0.0204768	−0.00417893	0.000689454	−0.0000769225	0.00000423651
	1.25438	−0.1235307	0.03467231	−0.01192915	0.00405709	−0.001191005	0.0002263204
v = 1500 m/s	1.57866	−0.296598	0.0949307	−0.0344762	0.0120067	−0.00344529	0.000605554

Table 1.

Staggered grid FD coefficient used to obtain the seismograms in Figure 4 with the space grid interval equals 20 m, and the time step equals 1.5 ms. In the first row is the traditional staggered grid FD coefficient obtained from Table 3 of Chu and Stoffa [17]; in the second row is the least squares staggered grid FD coefficient obtained from Table 3 of Liu [35]; and in the last rows are the staggered grid FD coefficients used for Eq. (30). Eqs. (31) and (32) use the simplest second-order staggered grid FD operator.

The seismograms recorded at positions A and B by different methods are presented in Figure 5 . Figure 5(a) is obtained with the traditional staggered grid FD scheme with the FD coefficient obtained in the space domain by Taylor expansion method [17]. The grid dispersion is obvious. Figure 5(b) is obtained with the traditional staggered grid FD scheme with the coefficient obtained in the space domain by the least squares method [35]. The staggered grid FD coefficient provided by Liu is one of the best staggered grid FD coefficient provided in recent years [35]. Figure 5(c) is obtained with the new staggered grid FD scheme with the coefficient determined in the time-space domain by the linear method. Figure 5(d) is obtained with the pseudo-spectrum method with the second-order acoustic wave equation. We observe that the grid dispersion in Figure 5(b) and (c) is similar to each other and is close to the nearly analytic results obtained with the pseudo-spectrum method in Figure 5(d) . However, the required simulation time is reduced by using the new staggered grid FD scheme because Eqs. (31) and (32) are much simpler than Eqs. (13) and (14).

Figure 5.

Seismograms recorded by different simulation methods. (a) The traditional staggered grid FD scheme with FD coefficients determined in the space domain by Taylor expansion method, (b) the traditional staggered grid FD scheme with FD coefficients determined in the space domain by least squares method, (c) the new staggered grid FD scheme with FD coefficients determined in the time-space domain by the linear method, and (d) the pseudo-spectrum method.

5.1.2. Numerical modeling in the salt model

Figure 6 shows the salt model from Society of Exploration of geophysicists with variations of velocities from 1486 to 4790 m/s. The seismic source function is the same as the previous example. The spatial sampling interval is 20 m, temporal step is 1 ms, and M = 7 for the staggered grid FD operators in Figure 7(a) and (b) . In Figure 7(c) , the parameters are M = 7 for the spatial derivatives in Eq. (1), and M = 1 for the spatial derivatives in Eqs. (2) and (3). The pseudo-spectrum method is used for the second-order acoustic wave equation as shown Figure 7(d) .

Figure 6.

SEG BP salt model.

Figure 7.

Seismic records obtained with different methods. (a) The traditional staggered grid FD scheme with FD coefficients determined in the space domain by Taylor expansion method, (b) the traditional staggered grid FD scheme with FD coefficients determined in the time-space domain by least squares method, (c) the new staggered grid FD scheme with FD coefficients determined in the time-space domain by the linear method, and (d) the pseudo-spectrum method.

Figure 7(a) is obtained the with the traditional staggered grid FD scheme with the coefficient obtained in the space domain by Taylor expansion method. The grid dispersion is obvious. Figure 7(b) is obtained with the traditional staggered grid FD scheme with the coefficient obtained in the time-space domain by the least squares method [27]. Most of the grid dispersion is suppressed. Figure 7(c) is obtained the with the new staggered grid FD scheme with the coefficient obtained in the time-space domain by the linear method. The grid dispersion in Figure 7(c) is very similar to the grid dispersion in Figure 7(b) . However, the simulation time to get Figure 7(c) is reduced compared with the simulation time to get Figure 7(b) . Both seismic records in Figure 7(b) and (c) are close to seismic record in Figure 7(d) . We want to mention that the linear method is faster than the LS method to determine the FD coefficients.

Figure 8 Further compares the seismograms in Figure 7 at position x/dx = 400. It is also observed that with the coefficient obtained in the space domain by Taylor expansion method, the grid dispersion is serious in the simulation result. The simulation results are almost overlapped for the traditional staggered grid FD scheme and new staggered grid FD scheme with optimized FD coefficient. However, the required simulation time is reduced by using the new staggered grid FD scheme because Eqs. (31) and (32) are much simpler than Eqs. (13) and (14).

Figure 8.

Seismograms at x/dx = 400 from Figure 4(a)–(d) .

Figure 9 compares snapshots of particle velocity v_x with the different staggered grid FD schemes at 2500 ms. it is also observed that with the coefficient obtained in the space domain by Taylor expansion method, the grid dispersion is most serious. The grid dispersion in Figure 9(c) is very similar to the grid dispersion in Figure 9(b) . It demonstrated that the new staggered grid FD scheme is accurate for the particle velocities in Eqs. (32) and (33) even when only second-order staggered grid FD operator is used.

Figure 9.

Particle velocity snapshots v_x obtained with different methods. (a) The traditional staggered grid FD scheme with the traditional FD coefficients, (b) the traditional staggered grid FD scheme with FD coefficients determined in the time-space domain by the least squares method, and (c) the new staggered grid FD scheme with FD coefficients determined in the time-space domain by the linear method.

5.2. Elastic wave equation

5.2.1. Numerical modeling in the homogeneous media

We first consider a homogeneous model. The P wave propagation speed is 2598 m/s, and the S wave velocity is 1500 m/s. The seismic source position is at the center of the model. The grid space interval is 20 m, the time step is 1 ms, and the operator length M is 7. A Ricker wavelet with the main frequency as 14.3 Hz was used as the seismic source.

The snapshots of the horizontal component obtained by different staggered grid FD methods are presented in Figure 10(a)–(c) . Figure 10(a) is obtained with the traditional staggered grid FD scheme with the traditional FD coefficient. The grid dispersion is obvious. Figure 10(b) is obtained with the traditional staggered grid FD scheme with the new FD coefficient. Compared with Figure 10(a) , the grid dispersion is suppressed. Figure 10(c) is obtained with the new staggered grid FD scheme. The grid dispersion curves in Figure 10(b) and (c) are very similar, which is further demonstrated in Figure 10(d) . However, with the new staggered grid FD scheme, we can save about 45% of the modeling time.

Figure 10.

Snapshots and slices of snapshots of the horizontal component at 698 ms obtained by different simulation methods. (a) The traditional staggered gird FD scheme with traditional FD coefficient, (b) the traditional staggered gird FD scheme with new FD coefficient, (c) the new staggered grid FD scheme with new FD coefficient, and (d) slices of snapshots at x/dx = 125.

5.2.2. Numerical modeling in the homogeneous media

Figure 11 shows the salt model from Society of Exploration of geophysicists. The S wave velocity is obtained from the P wave velocity. The seismic source function is plotted as a red asterisk. The spatial sampling interval is 12.5 m, the temporal step is 1 ms, and M = 7 for staggered grid FD operators.

Figure 11.

SEG BP salt model. (a) P wave velocity and (b) S wave velocity.

Figure 12 displays the seismic records of the horizontal component obtained by different staggered grid FD methods. Figure 12(a) is obtained with the traditional FD scheme with the traditional staggered grid FD coefficient. The grid dispersion is severe. Figure 12(b) is obtained with the traditional FD scheme with the staggered grid FD coefficient obtained by the least squares method. Figure 12(c) is obtained the with the new FD scheme with the staggered grid FD coefficient obtained by the linear method. It is observed that the grid dispersion in Figure 12(b) and 12(c) is smaller than the grid dispersion in Figure 12(a) . Figure 12(d) is seismograms obtained from Figure 12(a)–(c) . It further demonstrated that the grid dispersion in Figure 12(b) and (c) is similar to each other and smaller than the grid dispersion in Figure 12(a) . However, with the new FD scheme, the simulation time is reduced about 45%. In our simulation, there are 525 grids in the z direction and 850 grids in the x direction. With the traditional FD scheme, the simulation time is 920 seconds. With the new staggered FD grid scheme, the simulation time is 530 seconds. The huge reduction in simulation time is due to using the shorter staggered FD operator for the spatial derivatives in Eqs. (7)–(9). Figure 13 is the seismic records of the vertical component obtained by different FD methods. The same pattern can be observed from Figure 13(a)–(d) .

Figure 12.

Seismic records of the horizontal component obtained with different staggered grid FD methods. (a) Seismic records obtained by traditional FD scheme with traditional staggered grid FD coefficient, (b) seismic records obtained by the traditional FD scheme with new staggered grid FD coefficient, (c) seismic records obtained by the new FD scheme with new staggered grid FD coefficient, and (d) seismograms obtained from (a) to (c) at position x/dx = 355.

Figure 13.

Seismic records of the vertical component obtained with different staggered grid FD methods. (a) The traditional FD scheme with traditional staggered grid FD coefficient, (b) the traditional FD scheme with new staggered grid FD coefficient, (c) the new FD scheme with new staggered grid FD coefficient, and (d) seismograms obtained from (a) to (c) at position x/dx = 200.

6. Discussion and conclusion

The FD method is the most commonly used numerical method for wave equation modeling. Suppressing the grid dispersion is an important research area. Optimization method is usually used to determine the FD coefficients which could preserve the dispersion relation in a wider range of wavenumber (Zhang and Yao [24]; Ren and Liu [26]; Tan and Huang [32, 33]). We introduced the regularized optimization method to determine the FD coefficient which would be more robust for extreme conditions. The other way to suppress the grid dispersion is the utilization of the new FD scheme for the spatial derivatives. We introduce to use different FD operators for different spatial derivatives in the first-order wave equation. With the new staggered grid FD scheme, the wave equation modeling speed was accelerated while still preserving high accuracy. Through numerical modeling, we conclude that the introduced methods are more efficient while still preserving high accuracy for the first-order acoustic/elastic wave equation modeling. As a result, the introduced methods can be a substitute for the traditional FD methods used in acoustic/elastic wave equation modeling, which are essential in forward seismic wave modeling and reverse-time migration.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under grant numbers 41325016, 41704120, and 91630202.