1 INTRODUCTION

Seismic wave simulation plays an extremely important role in investigating seismic wave propagation in complex media. Various forward modeling methods have been developed in the literature,and generally they can be classified into three categories: geometry ray tracing methods^{[1, 2]},integral equation methods^{[3, 4]},and numerical methods^{[5, 6, 7, 8, 9]}. Geometry ray tracing methods are also called ray tracing as well,which are based on asymptotic solution of wave equation and mainly consider the high-frequent soltion. Integral equation methods treat the wave field of the scatters as the first or second kind unknown functions,and numerical methods are implemented to solve these function. Direct methods can obtain entire information of wave field and are still valid in various complex media. The most frequently used numerical methods include the finite difference method,Fourier pseudo-spectral method,finite element method,and spectral element method. All of these methods have their own advantages and disadvantages. Carcione et al.^[10] and Yang et al.^[11] discussed these methods in detail. Here we focus on the numerical methods.

The traditional finite difference method (FDM) is an effective tool for seismic wave simulations,but finite difference operators can not cover the high wave number domain,which is not ideal in simulation of high frequency seismic waves. The pseudo spectral method (PSM),which is based on fast Fourier transform,can cover the wave numbers including the Nyquist wave number,but lacks causality and does not agree with the characteristic of differentiation,such as the locality^[15],and its efficiency is relatively low as well. To circumvent the drawbacks of PSM,Mora^[16] and Holberg^[17] designed implicit a convolutional differentiator operator which preserves features of both PSM and FDM. Soon Zhou et al.^[18] developed an explicit convolutional differentiator operator and imposed window function on this operator to truncate the long operator to the short operator,and found that the finite operator length can achieve the same accuracy as PSM. Zhang et al.^[19] applied the truncated convolutional differentiator operator to simulation of seismic and acoustic wave fields and got reliable results. Dai et al.^[20] presented operator coefficients in power and Hamming windows. Long et al.^[21] proposed a method to choose optimal coefficients for window function. They first transformed derivative operator into wave number domain,and then compared analytical and numerical wave numbers,and finally obtained various optimal coefficients. Due to pursuing approximation to high accuracy in the high wave number domain,the error of the optimized convolutional differentiator operator is obvious in the low wave number domain.

Compared with the Shannon convolutional differentiator operator (SCD),the generalized Shannon convolutional differentiator (GCD) selects Gauss envelope to damp faster as the length of GCD increases,which means that GCD is extremely beneficial for embodying the locality of differentiation. We can use a shorter operator and a larger spatial increment to calculate spatial derivative to decrease computation costs. Wei^[22] first discussed the relationship between GCD and SCD,and carried out numerous numerical experiments and verified the accuracy of GCD. Qian^[23] deduced the truncation error of GCD and demonstrated the high accuracy and stability of this operator.

Here we use GCD to compute spatial derivatives of seismic wave equations,and exploit an effective method to optimize GCD. We transform the seismic wave equations into the Hamilton system,and introduce Lie operators^{[24, 25]} to construct exponential operators,and then use BCH formula^{[24, 25, 26]} for exponential operators expansion. Finally we obtain an optimal two-stage and second-order symplectic scheme for temporal evolution. For each time step,the generalized momentum and coordinate need to be updated only twice,which means that it is more efficient than the three-degree and third-order symplectic scheme. For symplectic schemes with one degree of freedom,Ruth^[27],Iwatsu^[28],and Li et al.^{[29, 30]} assigned specific value to one symplectic coefficient or the sum of two symplectic coefficients,and they did not adopt the optimized method to choose symplectic coefficients. Wang et al.^[31] proposed a method to obtain optimal coefficients based on one-dimensional oscillation equation,but the method may not be valid for wave equation. Here the optimized method based on the minimum error truncation has general significance. The optimal symplectic scheme is selected for spatial discretization and optiaml symplectic scheme is used for temporal discretization,and then the operator for seismic simulation is called as the optimal symplectic algorithm and generalized discrete convolutional differentiator (OSGCD). Compared with other schemes,such as PRK Lobatto III A-III B scheme and Iwatsu three-stage and third-order scheme,our scheme is more accurate than common second-order schemes,and is more efficient than third-order schemes. Various numerical results show that OSGCD is a good candidate for high precision seismic simulations.

2 CONSTRUCTION OF THE OSGCD METHOD 2.1 Structure-Preserving Temporal Discretization Scheme

Ignoring external force,seismic wave propagation in isotropic and elastic media can be described by the following equations in terms of displacements:

where (u,w,v) is three components of displacement; u,w are horizontal components; v is vertical component; ρ is mass density; and λ,μ are Lam′e constant.

Only considering elastic waves in two dimensions,Eq.(1) can be rewritten as

Eq.(2) can be written in the following simple form:

where is the second order derivative of displacement vector with respect to time.

Introducing generalized momentum p = $\dot U$ and generalized coordinate q = U,the 2n-dimensional phase space Z = (p,q) in terms of generalized momentum and coordinate advances in time domain can be expressed in the following Hamilton canonical equation:

where J is 2n-dimensional anti-symmetry matrix; J^-1 is the inversion of J; I_n is n-order identity matrix; H is Hamilton function which can be expressed as where T(p),V (q) are generalized kinetic and potential energy in the Hamilton system,respectively.

H represents the total energy of the dynamic system. If H is not related to time,system energy is conservative. Under the symplectic transform,the canonical form of Eq.(4) preserves. The basic theorem of Hamilton dynamics is that the phase at time t can be obtained by canonical transform in terms of H and phase at original time,which can be written as

where Z(0) is the initial phase; and Z(t) is the phase at time t.

Canonical transform Eq.(6) can preserves symplectic structure of the Hamilton system,that is to say,the 2-differential form does not change at each time step. From Eq.(5),we find that seismic wave equation is a separable Hamilton system. For a separable Hamilton system,it is easy to construct symplectic schemes with accuracy of any order^[38]. To construct symplectic schemes simply,we define Lie operators^{[24, 25]} A and B in terms of T(p) and V (q),respectively,which can be written as

where {} represents Poisson bracket.

Third or higher order schemes can improve accuracy,but huge computational cost hiders the implementation of high-order schemes in seismic simulations. In order to meet the need of practical applications,we design a second-order symplectic scheme as follows:

where a₁,a₂,b₁,b₂ are symplectic coefficients needed to be determined.

In Eq.(8) four coefficients need to satisfy three equations,thus Eq.(8) has one degree of freedom. Using BCH formula^{[24, 25, 26]} for exponential operator expansion,we obtain optimal coefficients based on minimum error truncation; the process is discussed in the appendix A. The optimal coefficients are

2.2 GCD Operator for Spatial Difference

The traditional FDM has a high efficiency in seismic wave simulations,but suffers from numerical dispersion when spatial increment is large^{[39, 40]}. Preserving the advantages of the FDM and suppressing its drawbacks is a focus of this paper. We introduce the generalized Shannon convolutional differential operator (GCD) to calculate spatial derivatives.

The generalized Shannon singular function can be defined as^[22]

where △ is spatial increment; and r is an arbitrary constant.

When △ → 0,Eq.(10) is a good candidate of Delta function. The generalized Shannon interpolation function can be written as

For an arbitrary discrete sequence {u(xi)},the spatial derivate of u(x) can be written as

Solving Eq.(11) numerically involves calculation of the first-order,second-order and mixed spatial derivatives. The first-order operator is

If the spatial increment is x = m△,Eq.(13) can be expressed as

We can obtain the second-order operator in a similar way,which is

Eqs.(14) and (15) are long operators,which involve the values at all spatial nodes. In practical calculation, we need to truncate long operators into short operators. To avoid the Gibbs phenomenon,we choose the Hanning window to truncate the operators,and the window function can be written as

where W is the length of an operator.

The truncated operator is

From Eq.(17),we can find that the first-order operator is symmetric,while the second-order operator is of anti-symmetry. To ensure the minimum error of group velocity and energy conservation,the first-order operator needs no change,while we must revise the value of the middle node of the second-order operator to be equal to the opposite value of the sum of the value of the other nodes^[20]. Three coefficients α,β,r need to be determined in Eq.(17). Upon the derivative approximating principle,we obtain optimal coefficients for Eq.(17),which are presented in Table1 and 2. The process of derivation can be found in Appendix B. Table1 and 2 show that the optimal coefficients are determined by the minimal sample rate and the operator length. In practical seismic simulations,we can choose different coefficients under specific situations to obtain high accuracy. In practical application,the coefficients of N_min = 5 can satisfy most cases. Pursuing the high accuracy in the high wave number domain and choosing smaller N_min will lead to larger errors in the low wave number domain. To balance the efficiency and accuracy,we choose the operator with W = 4,N_min = 5 in seismic wave simulations.

We use the first-order derivative of f(x) to analyze the accuracy of GCD in calculating spatial derivatives. f'(x) can be transformed into ik_xF(k_x) by Fourier transform. F(k_x) is phase of f(x) under Fourier transform,and the filter response of the spatial derivative under Fourier transform is ik_x. The filter response of GCD is

where C_n^(W) is the coefficients of GCD operator; and G = △/λ is the inverse of sampling rate in each wavelength.

The filter response rate of GCD is R_k =$\sum\limits_{n = 1}^W {\frac{{C_n^{(W)}}}{{\pi G}}} \sin (2\pi nG)$.Fig.1 shows the filter response rates of FDM and GCD,which indicates that GCD with radius equal to 4 (GCD4) is more accurate than FDM with radius equal to 5 (FDM5).

2.3 Stability Analysis of OSGCD

Using OSGCD for seismic simulations,the time interval and spatial increment should satisfy some proportional relation to keep the algorithm stable. Following literature ^[30],we define the stability condition as

where △t is time step; and a and C are constants need to be determined.

In Appendix C,derivation process is present in a 2D case,and the result is a = 0.6986,C = 0.8138.

3 NUMERICAL EXPERIMENTS 3.1 2D Infinite Homogeneous Model

We design a 2000 m×2000 m homogeneous model (Fig.2) to test the efficiency and accuracy of OSGCD in practical seismic wave simulations. Four methods are employed to compare with our method,which are (1) PRK Lobatto III A-III B^{[32, 33, 34, 35]} in combination with GCD (PRK-GCD),(2) Iwatsu three-stage and third-order symplectic scheme^{[15, 28, 29, 30, 36]} in combination with GCD (TO-GCD),(3) optimal two stage and second order symplectic scheme in combination with SCD^{[15, 28, 29, 41]} (OS-SCD) and (4) optimal second-degree and secondorder symplectic scheme in combination with FDM (OS-FEM). In this numerical experiment,we choose V_P=2000 m/s,V_S=1250 m/s and ρ=2000 kg/m3. A vertical point-focus source is used,which is Ricker wavelet with a peak amplitude spectrum at 30 Hz,and located at the center of the model. One receiver R₁ is located at (1200 m,1200 m). The spatial increment is 5 m,and the time increment is 0.5 ms. The numerical tests are implemented on the same computer with Intel (R) Core (TM) i5-2320 CPU and 16 GB RAM. The results are shown in Fig.3,Fig.4,and Table3, respectively.

From Fig.3a,one can find that the waveforms generated by OSGCD,PRK-GCD,and TO-GCD are very close to the analytical waveform,which means that numerical solutions modeled by these three methods can obtain credible results,while from Fig.3b,it is clear that the absolute errors of both OSGCD and TO-GCD are extremely close to each other,and very smaller than PRK-GCD; Table3 shows that TO-GCD consumes largest memory and CPU time to achieve high accuracy,which means that TO-GCD is less efficient than OSGCD in seismic simulations. We also compare the waveforms calculated by OSGCD,OS-SCD,and OS-FDM as shown in Fig.4. Fig.4a shows that three derivative methods are valid,but from absolute errors,as shown in Fig.4b, OSGCD is more accurate than the other two methods,while the computation costs are almost the same among these three methods. The comparisons demonstrate that OSGCD has great advantages over the other four methods in computational costs and numerical accuracy.

3.2 Long-Time Simulations

To investigate the ability of OSGCD in long-time simulations,a 2D 6000 m×6000 m homogeneous model with rigid boundaries is constructed. We use OSGCD and the popular Runge-Kutta FDM (RK-FDM) to simulate SH wave propagation. A shear source,which is a Ricker wavelet with a peak spectrum at 20 Hz,is located at the center of the model. The velocity of P-wave is 3000 m/s; the velocity of S-wave is 2000 m/s; the mass density is 2000 kg/m3. The spatial increment is 10 m and the time increment is 1 ms.

Figure5 displays the snapshots of the wavefield simulated by these two methods at different time instants. From Fig.5(a,b),it can be found that in short-time simulations these two methods result are in the largely same snapshots. But as the time increases,the wavefront in Fig.5c is blurred severely,while the wavefronts in Fig.5d are still quite clear,which means that OSGCD is quite suitable for long-time seismic simulations.

3.3 Wave Propagating in the Marmousi Model

The media of the real earth are strongly heterogeneous. It is necessary to investigate the ability of our methods for simulating seismic wave in strongly heterogeneous media. We choose the 3840 m×1210 m Marmousi model,which contains continuous and discontinuous elastic parameter changes,and the velocity profile of Pwave is shown in Fig.6. The velocity ranges from 1500 m/s to 5500 m/s; the S-wave velocity is V_S = V_P/$\sqrt 3 $ the mass density is obtained by Gardner formula ρ = 310 × V_P^{0.25 [43]}. A point-explosive source is located at centre of the model,which is a Ricker wavelet with an amplitude spectrum peak at 30 Hz. A receiver is located at (192,91) node point. The spatial interval is 10 m and the temporal interval is 1.0 ms. The PML boundary condition^[44] is adopted to avoid boundary reflection. We compare OSGCD with PSM in this experiment.

The snapshots of wave field at t=0.25 s in terms of the vertical displacement component are shown in Fig.7,and seismic records are presented in Fig.8. It can be found that there is no obvious difference between OSGCD and PSM in seismic wave simulation in the heterogeneous media. Because of the dramatic disturbance of elastic parameters of the model in Fig.6,we can see clearly the diffraction waves and scattered waves in Fig.7. This experiment indicates that OSGCD can generate creditable results when simulating seismic waves in strongly heterogeneous media.

4 CONCLUSIONS

Transforming the seismic wave equation in terms of displacement into the Hamiltonian system,we design a new symplectic scheme by means of Lie operators. The optimal coefficients are obtained based on minimum error truncation. We develop an optimal convolutional differential operator to calculate spatial derivatives. The stability conditions of OSGCD are presented. In numerical experiments,we compare our methods with traditional methods in various aspects. The results demonstrate that the performance of OSGCD is more efficient and can reduce computing costs dramatically. Our method is very suitable for long-time computation for its structure preserving properties. In simulating seismic waves in heterogeneous media,the results obtained by the new method are nearly the same as the PSM; the seismograms contain lots of information which is especially useful for determining geological structure.

5 Appendix A

Defining reciprocal relations^[45] [A,B] = AB - BA and [A,B,A] = [A,[B,A]],we have

By Baker-Campell-Hausdroff (BCH) formula, we obtain the expansion of W in the following form:

where a₁,a₂,b₁,b₂ must satisfy the following equivalent equations: where Minimum represents the sum of squares of η₁ and η₂,Minimum is minimum by choosing certain coefficients.

From (A3),we obtain ${a_1} = 1 - \frac{1}{{\sqrt 2 }},{a_2} = \frac{1}{{\sqrt 2 }},{b_1} = 1 - \frac{1}{{\sqrt 2 }},{b_2} = \frac{1}{{\sqrt 2 }}$.

6 Appendix B

Considering the following 2D harmonic plane wave:

where k_x = k sin θ,k_z = k cos θ; k =$\frac{w}{v}$ is wave number. θ denotes the plane wave incident angle to z axis.

The numerical derivative of (B1) calculated by GCD is ${\hat d_1}$ * P(x,z,w),and the analytic derivative of (B1) is P'(x,z,w). The numerical derivative should approach as close as possible to the analytic derivative,and we have the following object function:

where G = k△/2π = △/λ is the inverse of sample rate per wave length.

Eq.(B2) is nolinear optimization problem,and we construct the following iteration based on the steepest descent method:

where

7 Appendix C

We set the stability conditions of M1,M2,and M3 to be

where the spatial interval h = △x = △z; a and C are constants to be determined. The error of generalized coordinate and momentum at the nth time step is where E = exp[i(k_xx + k_zz)]; k_x,k_z are the wave number in x- and z- direction respectively. From Eq.(6),we can write the transfer error as where ; I₂ is two level identity matrix,Λ₂ is two level matrix and can be expressed as Λ₂ = (a_ij)_2×2,a₁₁ = V_P² k_x² + V_S ² k_z² ,a₁₂ = (V_P² - V_S ² )k_xk_z,a₂₁ = a₁₂,a₂₂ = V_S ² k_x² + V_P²k_z² ,k_x = k sin θ,k_z = k cos θ,k = $\frac{\pi }{\Delta }$ and 0 ≤ θ ≤ 2π. We define r =$\frac{{{V_P}}}{{{V_S}}}$,α = $\frac{{{V_P}\Delta t}}{h}$and G is a function of the variables r,α. If and only if the spectral radius of matrix G is less than or equal to 1,(C3) converges; so G must satisfy the following equation: We can solve (C4) numerically. When r increases from 1 to 4 with an increment of 0.01,the corresponding maximum α can be calculated. Finally,the constants a and C in (C1) are determined by the nonlinear fitting method. The result is

8 ACKNOWLEDGMENTS

We would like to thank the two anonymous reviewers for their careful reading and valuable comments on the manuscript. This work was supported by the National Natural Science Foundation of China (41174047, 40874024 and 41204041).