The numerical modeling of frequency domain scalar wave equation based on the rhombus stencil

Journal of Applied Geophysics 166 (2019) 42–56

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

The numerical modeling of frequency domain scalar wave equation based on the rhombus stencil Aman Li a,b,c,⁎, Hong Liu a,b,c, Yuxin Yuan a,b,c, Ting Hu a,b,c, Xiaona Cui a,b,c a b c

Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China Institutions of Earth Science, Chinese Academy of Science, Beijing 100029, China University of Chinese Academy of Sciences, Beijing 100049, China

a r t i c l e

i n f o

Article history: Received 7 November 2018 Received in revised form 28 April 2019 Accepted 1 May 2019 Available online 8 May 2019 Keywords: Numerical modeling Frequency domain Scalar wave Rhombus stencil

a b s t r a c t Frequency domain numerical modeling is the basis of frequency domain full waveform inversion. The accuracy of forward modeling seriously influences the result of inversion. In this paper, we develope a new method for frequency domain scalar wave equation to improve numerical simulation accuracy. The new method uses a rhombus stencil to approximate the Laplacian term and uses the weighted average of the difference scheme stencil points to approximate the mass acceleration term. Compared with the conventional stencil and the optimal square stencil, the numerical dispersion analysis shows that the new rhombus stencil can significantly decrease the grid points of per wavelength and reduce the numerical dispersion. Accuracy analysis demonstrates that the results of the new rhombus stencil are in best agreement with the analytical solutions relative to the other stencils in frequency domain and time domain. Two numerical examples are used to demonstrate the effectiveness of the new rhombus stencil for frequency domain scalar wave equation. © 2019 Elsevier B.V. All rights reserved.

1. Introduction In recent years, due to the improvement of computing ability, full waveform inversion (FWI) has developed rapidly and played almost the most important role in the research of exploration geophysics. FWI is used to iteratively improve a starting subsurface model by minimizing the difference between the modeled seismic data and the recorded seismic data (Tarantola, 1984). Because FWI can use all types of waves information from seismograms to rebuild high resolution imaging (Virieux and Operto, 2009). Many efforts have been focused on it. FWI can be basically classified into three categories: time domain (Virieux and Operto, 2009; Tarantola, 1984; Gauthier et al., 1986; Boonyasiriwat et al., 2009), frequency domain (Shin, 2006; Pratt and Worthington, 1990; Pratt, 1990; Plessix, 2009; Brossier et al., 2010) and hybrid domain (Sirgue et al., 2014; Liu et al., 2015; Kim et al., 2013), all of which have made important progress. FWI in frequency domain has some inherent advantages. First, the frequency domain FWI from low frequency to high frequency is a multiscale method itself. Second, only a few frequency components need to be selected for inversion (Pratt and Worthington, 1990), which can greatly improve the inversion efficiency. Moreover, when the impedance matrix is decomposed by the direct solver, it is suitable ⁎ Corresponding author at: Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China. E-mail address: [email protected] (A. Li).

https://doi.org/10.1016/j.jappgeo.2019.05.002 0926-9851/© 2019 Elsevier B.V. All rights reserved.

for multishot parallel computation. We all know that forward modeling is the basis of the inversion and it can largely affect the quality of inversion results. However, frequency domain modeling is problematic. When using the iterative method, we need find a good preconditioner to accelerate its convergence. When the direct method is used, the computational time and the amount of required memory become too large, especially for 3D. Here, we focus on designing more accuracy and efficiency stencil for direct solvers. Frequency domain forward modeling was pioneered by Lysmer and Drake (1972) using finite-element method to calculate the wave propagation in the earth. This method was developed further by Marfurt (1984) and Marfurt and Shin (1989). The finite difference method became more and more popular for solving geophysical problems due to its ease of implementation, parallelism and high efficiency (Virieux, 1986; Virieux, 1984; Alterman and Karal, 1968; Alford et al., 1974). Pratt and Worthington (1990) and Pratt (1990) applied finite difference forward modeling in frequency domain to the inversion and seismic imaging. They used the second order central finite difference operator to approximate the spatial derivative, so it required more grid points per wavelength compared with simulating in time domain for the same accuracy. By combining a rotated coordinate and the technique of mass acceleration term, Jo et al. (1996) developed an optimal 9-point scheme for the scalar wave equation. They used the steepest descent method to optimize the normalized phase curves to get the optimal coefficients. This scheme successfully reduces the grid points per wavelength to less than 4. Using the same method, Shin and Sohn (1998) combined

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

43

Fig. 1. The conventional cross stencil. (a) The conventional 5-point stencil (M = 1), (b) the conventional 9-point stencil (M = 2), (c) the conventional 13-point stencil (M = 3).

four rotated finite difference operators and developed an optimal 25-point scheme for the scalar wave equation. This new scheme reduces the grid points per wavelength to 2.5. Hustedt et al. (2004) and Operto et al. (2007) generalized the rotated coordinate method to variable density case and 3D case, respectively. It can be argued that those optimal schemes based on the rotated coordinate method have achieved considerable success. However, they require the equal spatial sampling intervals, which greatly affects their applications in practice. To overcome this limitation, Chen (2012) presented an optimal 9-point scheme based on average-derivate method (ADM) for acoustic wave equation. In this method, the finite difference of spatial derivative is approximated by the weighted average of three grid points in the orthogonal direction, which can be applied to different spatial sampling intervals, increasing the flexibility and widening its application range. Compared with the rotated coordinate method, average-derivate method has higher accuracy, so it has been widely used. Considering the calculation and precision, Tang et al. (2015) proposed an optimal ADM 17-point scheme and Zhang et al. (2014) developed an optimal ADM 25-point scheme. Chen (2014) directly extended average-derivate method to 3D case and proposed an optimal ADM 27-point scheme for 3D acoustic wave equation. Chen and Cao (2016) and Chen and Cao (2018) introduced

average-derivate method to 2D and 3D frequency domain elastic wave modeling, respectively. Chu and Stoffa (2012) used implicit finite difference operator for frequency domain numerical simulation to improve the accuracy of the simulation results. Liu et al. (2018) applied an optimization method to the implicit finite difference operator and developed a generic expression of it. In this paper, we present a new method for frequency domain scalar wave equation based on the rhombus stencil. The rhombus stencil was firstly introduced in Finkelstein and Kastner (2009), as a part case of their methodology which could achieve the spectral order of accuracy. Then Liu and Sen (2013) used the rhombus stencil for time space domain dispersion-relation-based finite difference. They demonstrated that new stencils could reach the same arbitrary even-order accuracy along all directions and were more accurate and more stable than conventional stencils. I introduce this stencil into the frequency domain numerical simulation to form a high precision method. The outline of the paper is as follows. We firstly present a new method based on rhombus stencil for frequency domain scalar wave equation. Then, we perform optimization of the coefficients for the rhombus stencil scheme and make dispersion analysis for some cases. This is followed by accuracy analysis and the tests of computational

44

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

Fig. 2. The optimal square stencil. (a) The optimal 9-point square stencil, (b) the optimal 25-point square stencil.

Fig. 3. The new rhombus stencil. (a) The 5-point rhombus stencil (M = 1), (b) the 13-point rhombus stencil (M = 2), (c) the 25-point rhombus stencil (M = 3).

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56 Table 1 The optimization coefficients for some value of M. M

1

2

3

a0, 0 a1, 0 a2, 0 a3, 0 a1, 1 a1, 2 a2, 1 b0, 0 b1, 0 b2, 0 b3, 0 b1, 1 b1, 2 b2, 1

−3.9992 0.9998 0 0 0 0 0 0.7356 0.0661 0 0 0 0 0

−1.9952 0.1375 0.0702 0 0.2911 0 0 0.4576 0.1195 −0.0010 0 0.0171 0 0

−0.7176 −0.1011 0.0741 0.0002 0.0766 0.0648 0.0648 0.2504 0.1176 0.0114 0.0004 0.0538 0.0021 0.0021

In order to improve accuracy and reduce frequency dispersion, researchers have studied the classical method in depth and proposed many optimization methods. To our knowledge, these optimization methods are basically based on the square scheme, as shown in Fig. 2. For the details about those methods, you can refer to the literatures (Shin and Sohn, 1998; Jo et al., 1996; Chen, 2012). Liu and Sen (2013) used a rhombus stencil for time space domain dispersion-relationbased finite difference and demonstrated that the rhombus stencil was more accurate and more stable than conventional stencil. Now, we introduce the rhombus stencil for frequency domain scalar wave equation shown in Fig. 3. The finite difference approximation formula for the Laplacian term can be expressed as ∇ P¼

2. Theory

In a 2D Cartesian coordinate system, the frequency domain acoustic wave equation in a homogeneous medium is ∇2 P þ

2

ω ∂ P ∂ P ω P ¼ 2 þ 2 þ 2 P ¼ 0; v2 v ∂x ∂z 2

ð1Þ

where Pis the pressure wavefield, v is the velocity of the propagation, ω is the circular frequency. For the conventional finite difference discretization shown in Fig. 1, the finite difference approximation formula to Eq. (1) is 1 2

h

" a0;0 P m;n þ ¼ 0;

M X i¼1

#   ω2 ai;0 P m−i;n þ P mþi;n þ P m;n−i þ P m;nþi þ 2 P m;n v

−i XM X 1 M−1 2

h

M X

#   ai;0 P m−i;n þ P mþi;n þ P m;n−i þ P m;nþi

i¼1

  ai; j P m−i;n− j þ P m−i;nþ j þ P mþi;n− j þ P mþi;nþ j ;

ð3Þ

i¼1 j¼1

where ai, j = aj, i. Converting Eq. (3) into the wavenumber domain and letting wavenumber be zero, we obtain the following relationship a0;0 þ 4

2

a0;0 P m;n þ

2

h

2.1. A new scheme based on rhombus stencil

2

"

1

2

þ efficiency. Finally, two numerical modeling examples are used to demonstrate the theoretical analysis.

45

M X

ai;0 þ 4

i¼1

M−1 −i X XM

ai; j ¼ 0:

ð4Þ

i¼1 j¼1

The mass acceleration term (ω2/v2)Pm, n is distributed at the collocation grid and its neighboring grids " M X   ω2 ω2 P ¼ b0;0 P m;n þ bi;0 P m−i;n þ P mþi;n þ P m;n−i þ P m;nþi m;n 2 2 v v i¼1 3 ð5Þ M−1 −i XM X   þ bi; j P m−i;n− j þ P m−i;nþ j þ P mþi;n− j þ P mþi;nþ j 5; i¼1 j¼1

where bi, j = bi, j, bi, jis the weighting coefficients of the mass acceleration term and satisfies the following relationship ð2Þ

where Pm, n = P(mh, nh), ai, 0is the coefficients of central finite difference, h is the grid size, square grid spacing is used in the discretization.

b0;0 þ 4

M X i¼1

bi;0 þ 4

M−1 −i X XM

bi; j ¼ 1:

i¼1 j¼1

Fig. 4. Normalized phase velocity curves of (a) the conventional 5-point stencil and (b) the 5-point rhombus stencil for different propagation anglesθ.

ð6Þ

46

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

Fig. 5. Normalized phase velocity curves of (a) the conventional 9-point stencil, (b) the Jo's optimal 9-point stencil, (c) the Chen's average 9-point stencil and (d) the 13-point rhombus stencil for different propagation anglesθ.

Substituting Eq. (3) and Eq. (5) into Eq. (1), we obtain the rhombus stencil scheme: " M X   ω2 b0;0 P m;n þ bi;0 P m−i;n þ P mþi;n þ P m;n−i þ P m;nþi v2 i¼1 þ 1 2

h

"

M−1 −i X XM

a0;0 P m;n þ

  bi; j P m−i;n− j þ P m−i;nþ j þ P mþi;n− j þ P mþi;nþ j 5þ

M X i¼1

þ

M−1 −i X XM

Substituting M = 1 into Eq. (7), we obtain the following scheme

3

i¼1 j¼1

ai; j



  ai;0 P m−i;n þ P mþi;n þ P m;n−i þ P m;nþi

Next, we examine some special cases. (1) M = 1

ð7Þ

3  P m−i;n− j þ P m−i;nþ j þ P mþi;n− j þ P mþi;nþ j 5 ¼ 0:

i¼1 j¼1

WhenMtakes different values in Eq. (7), it represents different schemes.

  ω2  b0;0 P m;n þ b1;0 P m−1;n þ P mþ1;n þ P m;n−1 þ P m;nþ1 þ v2   1 a0;0 P m;n þ a1;0 P m−1;n þ P mþ1;n þ P m;n−1 þ P m;nþ1 ¼ 0: 2 h

ð8Þ

Eq. (8) represents the 5-point rhombus stencil scheme, as shown in Fig. 3(a). It looks the same as the conventional 5-point scheme, but the mass acceleration term is different.

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

47

Fig. 6. Normalized phase velocity curves of (a) the conventional 13-point stencil, (b) the Shin's optimal 25-point stencil, (c) the Zhang's average 25-point stencil and (d) the 25-point rhombus stencil for different propagation anglesθ.

Substituting M = 3 into Eq. (7), we obtain the following scheme

(2) M ¼ 2 Substituting M = 2 into Eq. (7), we obtain the following scheme   ω2  b0;0 P m;n þ b1;0 P m−1;n þ P mþ1;n þ P m;n−1 þ P m;nþ1 v2   þ b2;0 P m−2;n þ P mþ2;n þ P m;n−2 þ P m;nþ2  þ b1;1 P m−1;n−1 þ P m−1;nþ1 þ P mþ1;n−1 þ P mþ1;nþ1 þ   1 a0;0 P m;n þ a1;0 P m−1;n þ P mþ1;n þ P m;n−1 þ P m;nþ1 2 h   þ a2;0 P m−2;n þ P mþ2;n þ P m;n−2 þ P m;nþ2  þ a1;1 P m−1;n−1 þ P m−1;nþ1 þ P mþ1;n−1 þ P mþ1;nþ1 ¼ 0:

ð9Þ

Eq. (9) represents the 13-point rhombus stencil scheme, as shown in Fig. 3(b). (3) M ¼ 3

  ω2  b0;0 P m;n þ b1;0 P m−1;n þ P mþ1;n þ P m;n−1 þ P m;nþ1 v2   þ b2;0 P m−2;n þ P mþ2;n þ P m;n−2 þ P m;nþ2  þ b3;0 P m−3;n þ P mþ3;n þ P m;n−3 þ P m;nþ3  þ b1;1 P m−1;n−1 þ P m−1;nþ1 þ P mþ1;n−1 þ P mþ1;nþ1  þ b1;2 P m−1;n‐2 þ P mþ1;n‐2 þ P mþ1;n−2 þ P mþ1;nþ2 þ b2;1 P m−2;n‐1 þ P mþ2;n‐1 þ P m‐2;nþ1 þ P mþ2;nþ1 þ   1 a0;0 P m;n þ a1;0 P m−1;n þ P mþ1;n þ P m;n−1 þ P m;nþ1 2 h   þ a2;0 P m−2;n þ P mþ2;n þ P m;n−2 þ P m;nþ2  þ a3;0 P m−3;n þ P mþ3;n þ P m;n−3 þ P m;nþ3  þ a1;1 P m−1;n−1 þ P m−1;nþ1 þ P mþ1;n−1 þ P mþ1;nþ1  þ a1;2 P m−1;n‐2 þ P mþ1;n‐2 þ P mþ1;n−2 þ P mþ1;nþ2 þ a2;1 P m−2;n‐1 þ P mþ2;n‐1 þ P m‐2;nþ1 þ P mþ2;nþ1 ¼ 0:

ð10Þ

48

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

Fig. 7. The geometry of the homogeneous model for (a) frequency-domain seismograms and (b) time-domain seismograms.

Fig. 8. The frequency-domain normalized seismograms of the conventional 9-point stencil, the Jo's optimal 9-point stencil, the Chen's average 9-point stencil the 13-point rhombus stencil and the analytical solutions in the homogeneous model for different frequencies. (a) f = 5HZ, (b) f = 10HZ, (c) f = 15HZ, (d) f = 20HZ.

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

49

Fig. 9. The time-domain normalized seismograms of the conventional 9-point stencil, the Jo's optimal 9-point stencil, the Chen's average 9-point stencil, the 13-point rhombus stencil and the analytical solutions in the homogeneous model for different geophone locations. (a) geophone 1, (b) geophone 2 and (c) geophone 3 of Fig. 7(b).

Eq. (10) represents the 25-point rhombus stencil scheme, as shown in Fig. 3(c). 2.2. Optimization and dispersion analysis In this section, we perform optimization of the coefficients for the rhombus stencil optimal scheme (7) and then make dispersion analysis for some special cases. Using the plane wave theory, we now perform the classical dispersion analysis by assuming a plane harmonic wave solution expressed as P ¼ P 0 eiðkx xþkz zÞ ;

ð11Þ

kx ¼ k sinθ; kz ¼ k cosθ;

ð12Þ

Where k is wavenumber, k x andk z are the horizontal and vertical wavenumber, and θ is the propagation direction angle of the plane wave. Substituting Eq. (11) into Eq. (7) and simplifying it,

we2have

3 M M−1 −i X XM X ω2 4 b0;0 þ 2 bi;0 ½ cosðikx hÞ þ cosðikz hÞ þ 4 bi; j ½ cosðikx hÞ cosðikz hÞ5þ v2 i¼1 i¼1 j¼1 2 3 M M−1 X X X M−i 14 a0;0 þ 2 ai;0 ½ cosðikx hÞ þ cosðikz hÞ þ 4 ai; j ½ cosðikx hÞ cosðikz hÞ5 ¼ 0: 2 h i¼1 i¼1 j¼1

ð13Þ The phase velocity vph is defined as ω : k

vph ¼

ð14Þ

Then, substituting Eq. (14) into Eq. (13), we can obtain the normalized phase velocities vph ¼ v vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XM XM−1 XM−i u u a0;0 þ 2 a ½ cosðikx hÞ þ cosðikz hÞ þ 4 a ½ cosðikx hÞ cosðikz hÞ G u i¼1 i;0 i¼1 j¼1 i; j t− ; XM XM−1 XM−i 2π b0;0 þ 2 b ½ cos ð ik h Þ þ cos ð ik h Þ  þ 4 b ½ cosðikx hÞ cosðikz hÞ x z i¼1 i;0 i¼1 j¼1 i; j

ð15Þ

50

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

Fig. 10. The normalized computation time of (A) the conventional 9-point stencil, (B) the Jo's optimal 9-point stencil, (C) the Chen's average 9-point stencil and (D) the 13-point rhombus stencil for the same accuracy.

where G = 2π/khis the number of grid points per wavelength. In order to get the optimization coefficients, we replace the expressions of kxh and kzh in Eq. (15) with kx h ¼

2π sinθ 2π cosθ ; kz h ¼ ; G G

ð16Þ

Then we can get ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XM XM−1 XM−i u u a þ 2 a A þ 4 a B 0;0 i;0 i; j vph G u i¼1 i¼1 j¼1 t− ¼ XM XM−1 XM−i 2π v b A þ 4 b B b0;0 þ 2 i;0 i¼1 j¼1 i; j     i¼1 2π sinθ 2π cosθ þ cos i A ¼ cos i G  G   2π sinθ 2π cosθ cos i B ¼ cos i G G

ð17Þ

The coefficients in Eq. (17) are determined by minimizing the normalized phase velocity error " E¼∬

1−

vph v

2 # d~ kdθ;

ð18Þ

where ~k ¼ 1=G. The range of the propagation direction angle θ is ~ is determined by the value of Mfor taken as[0, π/2]. The range of k different rhombus stencil scheme. There are many methods to solve the problem of (18). Here, we adopt a constrained nonlinear optimization program fmincon in MATLAB to determine the optimization coefficients. The optimization coefficients for some different value of Mare listed in Table 1. After getting the optimization coefficients, we make numerical dispersion analysis according to some most commonly used values ofM. Fig. 4 shows normalized phase velocity curves of the conventional 5-point stencil (M = 1) and the 5-point rhombus stencil (M = 1) for different propagation anglesθ. Within the phase velocity errors of 1%, the conventional 5-point stencil requires G ≈ 12.8, while the 5-point rhombus stencil only requires G ≈ 6.9. Fig. 5 shows normalized phase velocity curves of the conventional 9-point stencil (M = 2) and the 13-point rhombus stencil (M = 2) for different propagation anglesθ. For comparison, we also list the normalized phase velocity curves of the Jo's optimal 9-point stencil and the Chen's average 9-point stencil in Fig. 5. To keep the phase velocity errors in 1%, the conventional 9-point stencil requiresG ≈ 5.2, the Jo's optimal 9-point stencil requires G ≈ 3.2, the Chen's average 9-point stencil requires G ≈ 3.5, while the 13-point rhombus stencil only requires G ≈ 2.6. Similarly, Fig. 6 shows normalized phase velocity curves of the conventional 13-point stencil (M = 3), the Shin's optimal 25-point stencil, the Zhang's average 25-point stencil and the 25-point rhombus stencil (M = 3) for different propagation anglesθ. Maintaining the same standard that the phase velocity errors in 1%, the conventional 13-point stencil requiresG ≈ 3.9, the Shin's optimal 25-point stencil requires G ≈ 2.5, the Zhang's average 25-point stencil requiresG ≈ 2.7, while the 25-point rhombus stencil only requires G ≈ 2.1. From the above discussion, we find that the rhombus stencil method can reduce G by almost half compared with the conventional stencil method. We can also find that the rhombus stencil method has a better effect on reducing the phase velocity errors compared with the optimal square scheme. This means that our rhombus stencil method can reduce the computation cost with a coarse grid for the same accuracy.

Fig. 11. (a) The two-layer model. (b) 20HZ monochromatic wavefield computed by the 13-point rhombus stencil in frequency domain.

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

51

Fig. 12. Time domain seismograms of the two-layer model. Time domain seismograms computed by (a) the conventional 9-point stencil frequency domain modeling, (b) the Jo's optimal 9-point stencil frequency domain modeling, (c) the Chen's optimal 9-point stencil frequency domain modeling, (d) the 13-point rhombus stencil frequency domain modeling, and (f) the 2th-order in time and the 12th-order in space time domain modeling.

52

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

Fig. 13. single trace computed by different stencils at 300 m of the two-layer model.

2.3. Accuracy analysis The analytical solution of the frequency domain scalar wave equation (Morse et al., 1954) in a homogeneous medium is ð2Þ

P ðx; z; ωÞ ¼ −iπH 0

ω

r F ðωÞ; v

ð19Þ

where H(2) 0 is the second kind of Hankel function of zero order, F(ω) is the frequency domain source, i is the imaginary unit and r is the distance from the source to the receiver. Taking the inverse Fourier transform of Eq. (19), Alford et al. (1974) obtained the time domain analytical solution P ðx; z; t Þ ¼

1 2π

Z

þ∞

−∞

ð2Þ

−iπH 0

ω

r F ðωÞdω: v

ð20Þ

Then, we use the two formulas to make accuracy analysis. Here, the conventional 9-point stencil, the Jo's optimal 9-point stencil, the Chen's average 9-point stencil and the 13-point rhombus stencil are used as the comparison examples.

We perform 2D acoustic wave modeling in a homogeneous model with a velocity of 2000 m/s as shown in Fig. 7. The number of horizontal and vertical grid points are both 100. A Ricker wavelet with a peak frequency of 20HZ is placed at the center of the model as a source. Satisfied sampling theorem, the sampling interval is taken as h = 15m. When we calculate the frequency domain numerical solutions, we use the geometry of the homogeneous model as shown in Fig. 7(a). The geophone array is located at the surface. When we calculate the time domain numerical solutions, we change the geometry as shown in Fig. 7(b). The three geophones are separately placed at the left of the source, above the source, and along the diagonal direction. Firstly, let's look at the effect of single frequency. Fig. 8 shows the frequency domain normalized seismograms of the conventional 9-point stencil, the Jo's optimal 9-point stencil, the Chen's average 9-point stencil, the 13-point rhombus stencil and the analytical solutions in the homogeneous model for a series of different frequencies (5HZ, 10HZ, 15HZ, 20HZ). As shown in Fig. 8(a), when the frequency is 5HZ, all the above method can correctly do numerical simulations. When the frequency becomes 10HZ or 15HZ, all the above methods can almost correctly do numerical simulations, but the 13-point rhombus stencil has the highest precision, which can be seen from Fig. 8(b) or Fig. 8(c). When the frequency becomes 20HZ, the conventional 9-point stencil can't get the right numerical solution and the numerical solutions of the Jo's optimal 9-point stencil and the Chen's average 9-point stencil have big errors, however the numerical solution of the 13-point rhombus stencil is still in good agreement with the analytical solution. These phenomena can be explained by the following formula v ¼ λf ;

ð21Þ

where λ is the wavelength and f is the frequency. When the velocity is constant, the wavelength decreases as the frequency increases. Because the sampling interval is fixed, the decreasing wavelength is equivalent to the decreasing sampling point. This means that the sampling point decreases as the frequency increases. When the frequency is 5HZ, the sampling point can meet the requirement of all the above methods, so all methods have a good agreement with the analytical solution. As the frequency increases, the sampling point decreases, and only the method requiring less sampling points can have a better agreement with

Fig. 14. (a) The section of the SEG/EAGE salt model. (b) 30HZ monochromatic wavefield computed by the 13-point rhombus stencil in frequency domain.

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

53

Fig. 15. Time domain seismograms of the salt model. Time domain seismograms computed by (a) the conventional 9-point stencil frequency domain modeling, (b) the Jo's optimal 9-point stencil frequency domain modeling, (c) the Chen's optimal 9-point stencil frequency domain modeling, (d) the 13-point rhombus stencil frequency domain modeling, and (f) the 2thorder in time and the 12th-order in space time domain modeling.

54

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

3. Numerical examples We examine the validity of the rhombus stencil method on two models through numerical simulations. For consistency, we still use the 13-point rhombus stencil for modeling. For comparison, we do numerical simulation of the conventional 9-point stencil, the Jo's optimal 9-point stencil and the Chen's average 9-point stencil in frequency domain on two models. As reference, we also do time domain forward modeling with the 2th-order in time and the 12th-order in space on two models. 3.1. Two-layer model

Fig. 16. Single trace computed by different stencils at 1000 m of the salt model.

the analytical solution. For the methods that can't meet the sampling requirement, the discrepancy is related to the propagation angle. As shown in Fig. 8, the middle of the curves have a big discrepancy, which can be explained by the dispersion curves in Fig. 5. Then we have a look at the performance of the time domain. Fig. 9 shows the time domain normalized seismograms of the conventional 9-point stencil, the Jo's optimal 9-point stencil, the Chen's average 9-point stencil, the 13-point rhombus stencil and the analytical solutions in the homogeneous model for three geophones. It's obvious that the results of the 13-point rhombus stencil have the best agreement with the analytical solution compared with other methods. For the same method, Fig. 9(c) has a better agreement with the analytical solution compared with Fig. 9 (a) and Fig. 9(b). It's the fact that numerical dispersion is dependent on the propagation angles. As the Fig. 5 shows, the smaller the propagation angle is, the larger the discrepancy is. These phenomena in frequency domain and in time domain are all consistent with our above dispersion analysis, which proves the effectiveness of our method.

2.4. Computational efficiency Because the 13-point rhombus stencil needs more grid points around the discrete center point than other optimal methods, we need to test the computational efficiency. We use a homogeneous model with a velocity of 2000 m/s. The size of the model is 3000 m long and 3000 m wide. A Ricker wavelet with a peak frequency of 20HZ is used as a source. We test the computational time of different stencils while letting different stencils maintaining the same accuracy. Fig. 10 shows the normalized computation time for different stencils. Although the 13-point rhombus stencil needs more grid points around the discrete center point than other optimal methods, the computational cost is relatively smaller because it can reduce the grid points of per wavelength. To keep the phase velocity errors in 1%, the 13-point rhombus stencil reduces the total number of grid points to %48.8 of the conventional 9-point stencil, to %16.8 of the Jo's optimal 9-point stencil, and to %24.0 of the Chen's average 9-point stencil. So as shown in Fig. 10, the 13-point rhombus stencil has the smallest total computational time.

Firstly, we consider a two-layer model. The velocity structure of the two-layer model is shown in Fig. 11(a). The number of spatial grid points for horizontal and vertical direction are nx = 101 and nz = 101, respectively. We take the sampling interval to be 15 m. The velocity of the first layer in the two-layer model is 2000 m/s and the velocity of the second layer is 3500 m/s. The time sampling interval is 0.002 s and the recorded seismogram time is 1.0 s. Numerical simulation is often seriously influenced by artificial boundary reflection. We exploit the PML boundary condition (Berenger, 1994) to absorb the unwanted reflections. The thickness of the PML is 30. We place a Ricker wavelet with a dominant frequency of 20HZ at the middle of the two-layer model surface as a source. The geophone array is located at the third level of the model . Fig. 11(b) shows the 20HZ monochromatic wavefield computed by the 13-point rhombus stencil in frequency domain. It can be clearly seen the position of the reflection interface from the Fig. 11(b). The incident wave, reflected wave and transmitted wave are clearly visible near the interface of medium. The perturbation occurring above the interface is the result of mutual interference of the reflected wave and incident wave. Fig. 12 shows the time domain seismograms computed by the conventional 9-point stencil frequency domain modeling, the Jo's optimal 9-point stencil frequency domain modeling, the Chen's optimal 9-point stencil frequency domain modeling, the13-point rhombus stencil frequency domain modeling, and the 2th-order in time and the 12thorder in space time domain modeling for the two-layer model. From Fig. 12(a), we can find that the conventional 9-point stencil has a serious numerical dispersion. As shown in Fig. 12(b) and Fig. 12(c), compared with the conventional 9-point stencil, the Jo's optimal 9-point stencil and the Chen's optimal 9-point stencil frequency domain modeling can reduce the numerical dispersion. Fig. 12(d) shows that compared with the two above optimal stencil, the13-point rhombus stencil can more significantly reduce the numerical dispersion. Comparing Fig. 12 (d) with Fig. 12(f), you can see that the 13-point rhombus stencil frequency domain modeling and the 2th-order in time and the 12thorder in space time domain modeling have almost the same precision. Fig. 13 shows some single trace computed by different stencils at 300 m of the two-layer model, which can visibly illustrate the effectiveness of our method. 3.2. Salt model To further verify the validity of our method, we consider a two dimensional section of the SEG/EAGE salt model, shown in Fig. 14 (a). The velocity range of the salt model is from 1800 m/s to 4482 m/s. The number of spatial grid points for horizontal and vertical direction are nx = 400 and nz = 160, respectively. The spatial sampling interval is 13 m. The time sampling interval is 0.002 s and the recorded seismogram time is 3.2 s. We still use the PML boundary condition and the thickness of the PML is set to be 20. When

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

55

Fig. 17. Snapshots at 0.7 s. Snapshot computed by (a) the conventional 9-point stencil frequency domain modeling, (b) the Jo's optimal 9-point stencil frequency domain modeling, (c) the Chen's optimal 9-point stencil frequency domain modeling, (d) the 13-point rhombus stencil frequency domain modeling, and (f) the 2th-order in time and the 12th-order in space time domain modeling at 0.7 s.

we do numerical simulation, a Ricker wavelet with a dominant frequency of 20HZ is placed at nx = 200andnz = 2 as a source. The geophone array is located at nz = 1. Fig. 14(b) shows the 30HZ monochromatic wavefield computed by the 13-point rhombus stencil in frequency domain. Due to the irregular variable velocity interface, the mutual interference of the incident wave and reflected wave lead to a complicated frequency domain wavefields.

Fig. 15 shows the time domain seismograms computed by the conventional 9-point stencil frequency domain modeling, the Jo's optimal 9-point stencil frequency domain modeling, the Chen's optimal 9-point stencil frequency domain modeling, the 13-point rhombus stencil frequency domain modeling, and the 2th-order in time and the 12th-order in space time domain modeling for the salt model. Fig. 16 shows some single trace computed by different stencils at 1000 m of

56

A. Li et al. / Journal of Applied Geophysics 166 (2019) 42–56

the salt model. In order to see the effect of the numerical simulation more clearly, Fig. 17 shows the time domain snapshots corresponding to the above methods at 0.7 s, respectively. The irregular variable velocity interface also make the time domain seismograms and snapshots become very complex. From the results of the salt model, using the same above analysis method, we obtain the conclusion that our method can still effectively improve the accuracy and reduce the numerical dispersion for the complex medium. 4. Conclusions We have developed a new method based on the rhombus stencil for frequency domain scalar wave equation. We perform optimization to get the coefficients of the rhombus stencil optimal scheme. The dispersion analysis shows that compared with the conventional schemes, the rhombus stencil method can reduce the number grid points of per wavelength by almost half. Compared with the optimal square scheme, the rhombus stencil method has a better effect on reducing the phase velocity errors. The accuracy analysis shows that the results of the rhombus stencil scheme in frequency domain and time domain are both in best agreement with the analytical solutions. Because our method can significantly decrease the grid points of per wavelength, it can reduce the total computational time. Finally, two numerical modeling examples demonstrate that our method can effectively improve the accuracy and reduce the numerical dispersion. We mainly focus on the method study on the 2D frequency domain scalar wave equation in this paper. Due to the effectiveness of this method, we will extend it to the 3D frequency domain scalar wave equation in our future work. Acknowledgments This research is supported by Major State Research Development Program of China (Grant No. 2016YFC0601101), National Natural Science Foundation of China (Grant No. 41630319) and National Key Scientific Instrument and Equipment Development Project (Grant No. 2018YFF01013503). References Alford, R.M., Kelly, K.R., Boore, D.M., 1974. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics 39, 834–842. Alterman, Z., Karal, F.C., 1968. Propagation of elastic waves in layered media by finite difference methods. Bull. Seismol. Soc. Am. 58, 367–398. Berenger, J.P., 1994. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 114, 185–200. Boonyasiriwat, C., Valasek, P., Routh, P., Cao, W., Schuster, G.T., Macy, B., 2009. An efficient multiscale method for time-domain waveform tomography. Geophysics 74, WCC59. Brossier, R., Operto, S., Virieux, J., 2010. Which data residual norm for robust elastic frequency-domain full waveform inversion? Geophysics 75, R37–R46. Chen, J.B., 2012. An average-derivative optimal scheme for frequency-domain scalar wave equation. Geophysics 77, T201–T210. Chen, J.B., 2014. A 27-point scheme for a 3D frequency-domain scalar wave equation based on an average-derivative method. Geophys. Prospect. 62, 258–277.

Chen, J.B., Cao, J., 2016. Modeling of frequency-domain elastic-wave equation with an average-derivative optimal method. Geophysics 81, T339–T356. Chen, J.B., Cao, J., 2018. An average-derivative optimal scheme for modeling of the frequency-domain 3D elastic wave equation. Geophysics 83, T209–T234. Chu, C., Stoffa, P.L., 2012. An implicit finite-difference operator for the helmholtz equation. Geophysics 77, T97–T107. Finkelstein, B., Kastner, R., 2009. The spectral order of accuracy: a new unified tool in the design methodology of excitation-adaptive wave equation FDTD schemes. J. Comput. Phys. 228, 8958–8984. Gauthier, O., Virieux, J., Tarantola, A., 1986. Two-dimensional nonlinear inversion of seismic waveforms: numerical results. Geophysics 51, 1387. Hustedt, B., Operto, S., Virieux, J., 2004. Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling. Geophys. J. R. Astron. Soc. 157, 1269–1296. Jo, C.H., Shin, C., Suh, J.H., 1996. An optimal 9-point, finite-difference, frequency-space, 2D scalar wave extrapolator. Geophysics 61, 529–537. Kim, Y., Shin, C., Calandra, H., Min, D.J., 2013. An algorithm for 3D acoustic time-LaplaceFourier-domain hybrid full waveform inversion. Geophysics 78, R151–R166. Liu, Y., Sen, M.K., 2013. Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation. J. Comput. Phys. 232, 327–345. Liu, L., Ding, R., Liu, H., Liu, H., 2015. 3D hybrid-domain full waveform inversion on GPU. Comput. Geosci. 83, 27–36. Liu, Z., Song, P., Li, J., Li, J., Zhang, X., 2018. An optimized implicit finite-difference scheme for the two-dimensional Helmholtz equation. Geophys. J. Int. 202, 1805–1826. Lysmer, J., Drake, L.A., 1972. A finite element method for seismology. Methods Comput. Phys. Adv. Res. Appl. 11, 181–216. Marfurt, K.J., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics 49, 533–549. Marfurt, K.J., Shin, C.S., 1989. The future of iterative modeling in geophysical exploration. Handbook of Geophysical Exploration Seismic Exploration. vol. 21, pp. 203–228. Morse, P.M., Feshbach, H., Condon, E.U., 1954. Methods of theoretical physics parts I & II. Phys. Today 7, 15–16. Operto, S., Virieux, J., Amestoy, P., Giraud, L., L'Excellent, J.Y., 2007. 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: a feasibility study. Geophysics 72, 195. Plessix, R.E., 2009. Three-dimensional frequency-domain full-waveform inversion with an iterative solver. Geophysics 74, WCC149–WCC157. Pratt, R.G., 1990. Frequency-domain elastic wave modeling by finite differences: a tool for crosshole seismic imaging. Geophysics 55, 626–632. Pratt, R.G., Worthington, M.H., 1990. Inverse theory applied to multi-source cross-hole tomography, part I : acoustic wave-equation method. Geophys. Prospect. 38, 287–310. Shin, C., 2006. Waveform inversion using a logarithmic wavefield. Geophysics 71, 13–22. Shin, C., Sohn, H., 1998. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator. Geophysics 63, 289–296. Sirgue, L., Etgen, J.T., Albertin, U., 2014. 3D Frequency domain waveform inversion using time domain finite difference methods. Eage Conference and Exhibition Incorporating Spe Europec. Tang, X., Liu, H., Zhang, H., Liu, L., Wang, Z., 2015. An adaptable 17-point scheme for highaccuracy frequency-domain acoustic wave modeling in 2D constant density media. Geophysics 80, T211–T221. Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics 49, 1259–1266. Virieux, J., 1984. SH-wave propagation in heterogeneous media: velocity-stress finitedifference method. Geophysics 49, 1933–1942. Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method. Geophysics 51, 889–901. Virieux, J., Operto, S., 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics 74, WCC1–WCC26. Zhang, H., Liu, H., Liu, L., Jin, W., Shi, X., 2014. Frequency domain acoustic equation highorder modeling based on an average-derivative method (in Chinese). Chin. J. Geophys. 57, 1599–1611.