Open access peer-reviewed chapter

# Optimized Finite Difference Methods for Seismic Acoustic Wave Modeling

Written By

Yanfei Wang and Wenquan Liang

Reviewed: October 12th, 2017 Published: December 20th, 2017

DOI: 10.5772/intechopen.71647

From the Edited Volume

## Computational and Experimental Studies of Acoustic Waves

Edited by Mahmut Reyhanoglu

Chapter metrics overview

View Full Metrics

## 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.

### Keywords

• finite difference scheme
• optimized finite difference coefficient
• staggered grid
• regularization
• wave equation

## 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 [17].

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

P t = v 2 v z z + v x x , E1
v x t = P x , E2
v z t = P z . E3

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

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

2 p x 2 + 2 p z 2 = 1 v 2 2 p t 2 E4

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

v x t = τ xx x + τ xz z , E5
v z t = τ xz x + τ zz z , E6
τ xx t = α 2 v x x + α 2 2 β 2 v z z , E7
τ zz t = α 2 v z z + α 2 2 β 2 v x x , E8
τ xz t = β 2 v x z + v z x . E9

where (vx , vz ) 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

2 v x t 2 = α 2 2 v x x 2 + α 2 β 2 2 v z x z + β 2 2 v x z 2 , E10
2 v z t 2 = α 2 2 v z z 2 + α 2 β 2 2 v x x z + β 2 2 v z x 2 . E11

## 3. Finite difference operators

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

P t = v 2 h m = 1 M 1 c m v z 0 , m 1 / 2 v z 0 , m + 1 / 2 + m = 1 M 1 c m v x m 1 / 2 , 0 v x m + 1 / 2 , 0 , E12
v x t = 1 h m = 1 M 2 c m P m 1 / 2 , 0 0 P m + 1 / 2 , 0 0 , E13
v z t = 1 h m = 1 M 2 c m P 0 , m 1 / 2 0 P 0 , m + 1 / 2 0 , E14
Q m , j n = Q x + mh z + jh t + ; Q = v x , v z , P E15

where M 1 and M 2 are the length of the FD operators, cm 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:

P t = 1 Δt P 0 , 0 1 P 0 , 0 0 . E16

where Δt is the time step.

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

v x t = 1 h m = 1 M 1 c m [ τ x x m 1 / 2 , 0 0 τ x x m + 1 / 2 , 0 0 ] + 1 h m = 1 M 1 c m [ τ x z 0 , m 1 / 2 0 τ x z 0 , m + 1 / 2 0 ] E17
v z t = 1 h m = 1 M 1 c m [ τ x z m 1 / 2 , 0 0 τ x z m + 1 / 2 , 0 0 ] + 1 h m = 1 M 1 c m [ τ z z 0 , m 1 / 2 0 τ z z 0 , m + 1 / 2 0 ] E18
τ x x t = α 2 h m = 1 M 2 c m [ v x m 1 / 2 , 0 0 v x m + 1 / 2 , 0 0 ] + α 2 2 β 2 h m = 1 M 2 c m [ v z 0 , m 1 / 2 0 v z 0 , m + 1 / 2 0 ] E19
τ z z t = α 2 2 β 2 h m = 1 M 2 c m [ v x m 1 / 2 , 0 0 v x m + 1 / 2 , 0 0 ] + α 2 h m = 1 M 2 c m [ v z 0 , m 1 / 2 0 v z 0 , m + 1 / 2 0 ] E20
τ x z t = β 2 h m = 1 M 2 c m [ v z m 1 / 2 , 0 0 v z m + 1 / 2 , 0 0 ] + β 2 h m = 1 M 2 c m [ v x 0 , m 1 / 2 0 v x 0 , m + 1 / 2 0 ] E21
Q m , j n = Q x + mh z + jh t + ; Q = v x , v z , τ xx , τ zz , τ xz E22

where M 1 and M 2 are the length of the FD operators, cm 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 k x x + mh + k z z + jh ω t + . E23

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

m = 1 M c m sin m 0.5 k x h 2 + m = 1 M c m sin m 0.5 k z h 2 1 2 r 2 1 cos kvτ . E24

where r = vΔt/h, M1 = M2 = M, and(kx , kz ) = 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 c = θ = 0 2 π m = 1 M c m sin m 0.5 k x h 2 + m = 1 M c m sin m 0.5 k z h 2 E25

and the right side of Eq. 24 by

d = θ = 0 2 π 1 2 r 2 1 cos kvτ . E26

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

Φ c = k = 0 K F c d 2 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]:

Ratio = K K total = 2 π f max / v π / h = f max v / 2 h . 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 α c = Φ c + 1 2 α 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,

P t = v 2 h m = 1 M c m v z 0 , m 1 / 2 v z 0 , m + 1 / 2 + m = 1 M c m v x m 1 / 2 , 0 v x m + 1 / 2 , 0 , E30
v x t = P 1 / 2 , 0 0 P 1 / 2 , 0 0 h , E31
v z t = P 0 , 1 / 2 0 P 0 , 1 / 2 0 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 vx and vz 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

( e i k z h 2 e i k z h 2 ) m=1 M c m ( e i k z ( m 1 2 )h e i k z ( m 1 2 )h ) +( e i k x h 2 e i k x h 2 ) m=1 M c m ( e i k x ( m 1 2 )h e i k x ( m 1 2 )h )= h 2 v 2 e iωτ + e iωτ 2 Δ t 2 E33

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

2 sin k z h 2 m = 1 M c m sin m 0.5 k z h 2 sin k x h 2 m = 1 M c m sin m 0.5 k x h = r 2 cos ωτ 1 . E34

Using the basic trigonometric function

sin α sin β = cos α + β cos α β 2 , E35

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

m = 1 M c m cos mk x h cos m 1 k x h + cos mk z h cos m 1 k z h = r 2 cos ωτ 1 . E36

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

m = 1 M c m cos mk x h cos m 1 k x h + cos mk y h cos m 1 k y h + cos mk z h cos m 1 k z h = r 2 cos ωτ 1 E37

where k x k y k z = k sin θ cos φ sin θ sin φ cos θ .

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]:

φ = 0 2 π θ = 0 2 π a k 1 , x , 1 h + a k 1 , y , 1 h + a k 1 , z , 1 h 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 a k M , x , M h + a k M , y , M h + a k M , x , M h c 1 c M = 1 r 2 φ = 0 2 π θ = 0 2 π cos k 1 1 cos k M 1 E38

where a k l , m h = cos mk l h cos m 1 k l h , the ith component of kl (l = x, y, z) is represented as k i, l , kx  = k cos θ cos φ, ky  = k cos θ sin φ, kz  = 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

δ = v FD v = 1 rkh arccos 1 + r 2 q 2 . E39

where

q 2 = m = 1 M c m cos mk x h cos m 1 k x h + cos mk z h cos m 1 k z h E40

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

ε = h v FD h v = h v v v FD 1 = h v 1 δ 1 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−5 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).

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

g = m = 1 M c m cos mk x h cos m 1 k x h + cos mk z h cos m 1 k z h . E42

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

g < 0 ; r 2 g 2 . E43

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

r 2 g 2 m = 1 M 4 c m 1 m = 1 2 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.

### 4.2. Optimizing finite difference operators for the elastic wave equation

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

α 2 D xx + β 2 D zz D tt α 2 β 2 D xz α 2 β 2 D xz α 2 D zz + β 2 D xx D tt v x v z = 0 . E45

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

D tt = 1 2 α 2 + β 2 D xx + D zz ± 1 2 α 2 β 2 D xx + D zz 2 4 D xx D zz D xz D xz . E46

Suppose D xx Dzz = D xz D xz , then two equations can be obtained from Eq. (46):

α 2 D xx + α 2 D zz D tt = 0 . E47
β 2 D xx + β 2 D zz D 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):

m = 1 M c m sin m 0.5 k x h 2 + m = 1 M c m sin m 0.5 k z h 2 1 2 r 2 1 cos kvτ . E49

where r = βΔt/h, M1 = M2 = M, and (kx , kz ) = 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):

v x t = 1 h m = 1 M 1 c m [ τ x x m 1 / 2 , 0 0 τ x x m + 1 / 2 , 0 0 ] + 1 h m = 1 M 1 c m [ τ x z 0 , m 1 / 2 0 τ x z 0 , m + 1 / 2 0 ] E50
v z t = 1 h m = 1 M 1 c m [ τ x z m 1 / 2 , 0 0 τ x z m + 1 / 2 , 0 0 ] + 1 h m = 1 M 1 c m [ τ z z 0 , m 1 / 2 0 τ z z 0 , m + 1 / 2 0 ] E51
τ x x t = α 2 h [ v x 1 / 2 , 0 0 v x 1 / 2 , 0 0 ] + α 2 2 β 2 h [ v z 0 , 1 / 2 0 v z 0 , 1 / 2 0 ] E52
τ x x t = α 2 2 β 2 h [ v x 1 / 2 , 0 0 v x 1 / 2 , 0 0 ] + α 2 h [ v z 0 , 1 / 2 0 v z 0 , 1 / 2 0 ] E53
τ x z t = β 2 h [ v z 1 / 2 , 0 0 v z 1 / 2 , 0 0 ] + β 2 h [ v x 0 , 1 / 2 0 v x 0 , 1 / 2 0 ] 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:

m = 1 M 1 c m cos mk x h cos m 1 k x h + cos mk z h cos m 1 k z h = r 2 cos ωτ 1 . 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 .

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).

#### 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 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 9 compares snapshots of particle velocity vx 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.

### 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.

#### 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 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) .

## 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.

## References

1. 1. Alterman Z, Karal FC. Propagation of elastic waves in layered media by finite difference methods. Bulletin of the Seismological Society of America. 1968;58(1):367-398
2. 2. Virieux J. SH-wave propagation in heterogeneous media: Velocity-stress finite- difference method. Geophysics. 1984;49(11):1933-1942
3. 3. Virieux J. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics. 1986;51(4):889-901
4. 4. Dablain MA. The application of high-order differencing to the scalar wave equation. Geophysics. 1986;51(1):54-66
5. 5. Robertsson JO, Blanch JO, Symes WW. Viscoelastic finite-difference modeling. Geophysics. 1994;59(9):1444-1456
6. 6. Bohlen T, Wittkamp F. Three-dimensional viscoelastic time-domain finite-difference seismic modeling using the staggered Adams-Bashforth time integrator. Geophysical Journal International. 2016;204(3):1781-1788
7. 7. Etemadsaeed L, Moczo P, Kristek J, Ansari A, Kristekova M. A no-cost improved velocity-stress staggered-grid finite-difference scheme for modelling seismic wave propagation. Geophysical Journal International. 2016;207(1):481-511
8. 8. Alford RM, Kelly KR, Boore DM. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics. 1974;39(6):834-842
9. 9. Kelly KR, Ward RW, Treitel S, Alford RM. Synthetic seismograms: A finite-difference approach. Geophysics. 2012;41:2-27
10. 10. Madariaga R. Dynamics of an expanding circular fault. Bulletin of the Seismological Society of America. 1976;66(3):639-666
11. 11. Levander AR. Fourth-order finite-difference p-sv seismograms. Geophysics. 1988;53(11):1425-1436
12. 12. Yang D, Teng J, Zhang Z, Liu E. A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media. Bulletin of the Seismological Society of America. 2003;93(2):882-890
13. 13. Yang D, Lu M, Wu R, Peng J. An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations. Bulletin of the Seismological Society of America. 2004;94(5):1982-1992
14. 14. Chen JB. High-order time discretizations in seismic modeling. Geophysics. 2007;72(5):SM115-SM122
15. 15. Chen JB. A stability formula for lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation. Geophysics. 2011;76(2):T37-T42
16. 16. Song X, Fomel S. Fourier finite-difference wave propagation. Geophysics. 2011;76(5):T123-T129
17. 17. Chu C, Stoffa PL. Determination of finite-difference weights using scaled binomial windows. Geophysics. 2012;77(3):W17-W26
18. 18. Fomel S, Ying L, Song X. Seismic wave extrapolation using low-rank symbol approximation. Geophysical Prospecting. 2013;61(3):526-536
19. 19. Finkelstein B, Kastner R. Finite difference time domain dispersion reduction schemes. Journal of Computational Physics. 2007;221:422-438
20. 20. Finkelstein B, Kastner R. A comprehensive new methodology for formulating FDTD schemes with controlled order of accuracy and dispersion. IEEE Transactions on Antennas and Propagation. 2008;56:3516-3525
21. 21. Etgen J. T. A tutorial on optimizing time domain finite-difference scheme: Beyond Holberg: Stanford Exploration Project Report. 2007;129:33-43
22. 22. Liu Y, Sen MK. A new time-space domain high-order finite difference method for the acoustic wave equation. Journal of Computational Physics. 2009;228:8779-8806
23. 23. Liu Y, Sen MK. Scalar wave equation modeling with time-space domain dispersion-relation-based staggered-grid finite-difference schemes. Bulletin of the Seismological Society of America. 2011;101(1):141-159
24. 24. Zhang JH, Yao ZX. Optimized finite-difference operator for broadband seismic wave modeling. Geophysics. 2013;78:A13-A18
25. 25. Liang WQ, Yang CC, Wang YF, Liu HW. Acoustic wave equation modeling with new time-space domain finite difference operators. Chinese Journal of Geophysics. 2013;56(6):840-850
26. 26. Ren Z, Liu Y. Acoustic and elastic modeling by optimal time-space-domain staggered-grid finite-difference schemes. Geophysics. 2014;80(1):T17-T40
27. 27. Wang Y, Liang W, Nashed Z, Li X, Liang G, Yang C. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion- relationship-preserving method. Geophysics. 2014;79(5):T277-T285
28. 28. Chen H, Zhou H, Zhang Q, Chen Y. Modeling elastic wave propagation using K space operator-based temporal high-order staggered-grid finite-difference method. IEEE Transactions on Geoscience and Remote Sensing. 2017;55(2):801-815
29. 29. Yong P, Huang J, Li Z, Liao W, Qu L, Li Q, Liu P. Optimized equivalent staggered-grid FD method for elastic wave modelling based on plane wave solutions. Geophysical Journal International. 2017;208(2):1157-1172
30. 30. Liu Y, Sen MK. Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation. Journal of Computational Physics. 2013;232(1):327-345
31. 31. Liu H, Dai N, Niu F, Wu W. An explicit time evolution method for acoustic wave propagation. Geophysics. 2014;79(3):T117-T124
32. 32. Tan S, Huang L. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation. Geophysical Journal International. 2014a;197(2):1250-1267
33. 33. Tan S, Huang L. A staggered-grid finite-difference scheme optimized in the time-space domain for modeling scalar-wave propagation in geophysical problems. Journal of Computational Physics. 2014b;276:613-634
34. 34. Margrave GF. Numerical Methods of Exploration Seismology with Algorithms in Matlab. Calgary: Department of Geology and Geophysics, the University of Calgary; 2001
35. 35. Liu Y. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling. Geophysical Journal International. 2014;197(2):1033-1047

Written By

Yanfei Wang and Wenquan Liang

Reviewed: October 12th, 2017 Published: December 20th, 2017