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.
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 [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
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
where
Substituting Eqs. (2) and (3) into Eq. (1), the second-order acoustic wave equation can be written as
The first-order elastic wave equations in 2D heterogeneous media are [3]
where (
Substituting Eqs. (7)–(9) into Eqs. (5)–(6), the second-order elastic wave equation can be written as
3. Finite difference operators
The commonly used staggered grid FD scheme for the first-order acoustic wave equation is as follows:
where
The second-order FD operator is usually used for the first-order time derivative:
where Δ
The commonly used staggered grid FD scheme for the first-order elastic wave equation is as follows:
where
4. Optimizing finite difference operators
4.1. Optimizing finite difference operators for the acoustic wave equation
Using the plane wave theory, let
The following dispersion relation can be obtained by substituting Eqs. (13)–(14) into Eq. (12) [23, 26, 27]:
where
Let
and the right side of Eq. 24 by
The aim is to minimize the dispersion error for a fixed range of wavenumbers:
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]:
The direct minimization of the objective function
where
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,
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
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
From Eq. (33), the following dispersion relation can be obtained in the frequency-wavenumber domain (it is a special case of Eq. (24)):
Using the basic trigonometric function
we obtain Eq. (36) from Eq. (34):
Similarly, the new dispersion relation for the 3D first-order acoustic wave equation in the frequency-wavenumber domain is
where
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
where
The 2D dispersion error
where
The difference between the FD propagation time and the exact propagation time through one grid is defined as [23]
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
Let the left part of Eq. (36) as.
From dispersion relation Eq. (36), it is obvious that.
Then, the stability condition of the new staggered grid FD scheme is (from Eq. (36) with
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
4.2. Optimizing finite difference operators for the elastic wave equation
Eqs. (10) and (11) can be written as [11].
The two roots give the following dispersion relation [11]:
Suppose
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):
where
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):
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 (
Then, the new dispersion relation can be obtained from Eq. (49) in the frequency-wavenumber domain:
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
c_{1} | c_{2} | c_{3} | c_{4} | c_{5} | c_{6} | c_{7} | |
---|---|---|---|---|---|---|---|
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 |
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
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
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
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
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.
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.
Virieux J. SH-wave propagation in heterogeneous media: Velocity-stress finite- difference method. Geophysics. 1984; 49 (11):1933-1942 - 3.
Virieux J. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics. 1986; 51 (4):889-901 - 4.
Dablain MA. The application of high-order differencing to the scalar wave equation. Geophysics. 1986; 51 (1):54-66 - 5.
Robertsson JO, Blanch JO, Symes WW. Viscoelastic finite-difference modeling. Geophysics. 1994; 59 (9):1444-1456 - 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.
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.
Alford RM, Kelly KR, Boore DM. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics. 1974; 39 (6):834-842 - 9.
Kelly KR, Ward RW, Treitel S, Alford RM. Synthetic seismograms: A finite-difference approach. Geophysics. 2012; 41 :2-27 - 10.
Madariaga R. Dynamics of an expanding circular fault. Bulletin of the Seismological Society of America. 1976; 66 (3):639-666 - 11.
Levander AR. Fourth-order finite-difference p-sv seismograms. Geophysics. 1988; 53 (11):1425-1436 - 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.
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.
Chen JB. High-order time discretizations in seismic modeling. Geophysics. 2007; 72 (5):SM115-SM122 - 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.
Song X, Fomel S. Fourier finite-difference wave propagation. Geophysics. 2011; 76 (5):T123-T129 - 17.
Chu C, Stoffa PL. Determination of finite-difference weights using scaled binomial windows. Geophysics. 2012; 77 (3):W17-W26 - 18.
Fomel S, Ying L, Song X. Seismic wave extrapolation using low-rank symbol approximation. Geophysical Prospecting. 2013; 61 (3):526-536 - 19.
Finkelstein B, Kastner R. Finite difference time domain dispersion reduction schemes. Journal of Computational Physics. 2007; 221 :422-438 - 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.
Etgen J. T. A tutorial on optimizing time domain finite-difference scheme: Beyond Holberg: Stanford Exploration Project Report. 2007; 129 :33-43 - 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.
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.
Zhang JH, Yao ZX. Optimized finite-difference operator for broadband seismic wave modeling. Geophysics. 2013; 78 :A13-A18 - 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.
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.
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.
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.
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.
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.
Liu H, Dai N, Niu F, Wu W. An explicit time evolution method for acoustic wave propagation. Geophysics. 2014; 79 (3):T117-T124 - 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.
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.
Margrave GF. Numerical Methods of Exploration Seismology with Algorithms in Matlab. Calgary: Department of Geology and Geophysics, the University of Calgary; 2001 - 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