1 INTRODUCTION

Due to its simplicity, flexibility, and cost-efficiency, the finite-difference (FD) algorithm, especially the staggered-grid finite-difference algorithm, has been widely used in the area of seismic modeling (Virieux, 1984, 1986; Pei, 2005; Yan and Liu, 2012), reverse time migration (RTM) (Baysal et al., 1983; McMechan, 1983; Yan and Liu, 2013) and full waveform inversion (FWI) (Warner et al., 2013; Virieux and Operto, 2009). At present, the finite-difference, especially the staggered-grid finite-difference scheme, has gained general acceptance because of its high accuracy and stability. Normally, the finite-difference scheme can be divided into two parts, the explicit one and the implicit one.

Traditionally, the explicit finite-difference schemes are derived from the Taylor’s theorem or the optimized algorithms (Liu et al., 1998; Dong et al., 2000; Liu, 2014). In order to achieve highly accurate numerical solution, they usually utilize longer stencils to achieve high stability and numerical accuracy. What’s more, according to Kosloff et al. (2010), there exists a saturation effect whereby the rate of improvement per order of the difference operator decreases as the increasing of orders. And experience shows that the most economical order lies somewhere between 4 and 8. Therefore, in recent years, much attention has been paid to the implicit finite-difference operators, because the implicit operators could utilize smaller stencils to achieve high stability and numerical accuracy (Kosloff et al., 2010).

Unlike the explicit finite-difference schemes, the implicit ones are based on the rational expansions, which converge faster than power series expansions (Fornberg, 1998). Although implicit difference operators had been widely used in many fields (Lele, 1992; Ekaterinaris, 1999; Shang, 1999), they were just introduced to solve the problem of seismic wave propagation in recent years (Liu and Sen, 2009a, 2009b; Kosloff et al., 2010; Chu and Stoffa, 2010, 2012). Kosloff et al. (2008) introduced the implicit finite-difference algorithm into the calculation of the second-order spatial derivatives of time-domain acoustic wave equation for constant density media. Afterwards, Kosloff et al. (2010) derived the first-order implicit staggered-grid finite-difference operators and applied the modeling results to RTM. The modeling and RTM performances prove that this operator can decrease numerical dispersion effectively. Chu and Stoffa (2010) proposed an implicit algorithm to calculate the second-order spatial derivatives of frequency-space domain acoustic wave equation. Afterwards, Chu and Stoffa (2012) derived a new modeling method which combines the second-order implicit split-time integration and the second-order spatial derivatives implicit finite-difference method.

However, the conventional implicit finite-difference operators mentioned previously generally have to obtain the high accuracy solutions from the inversion of either pentadiagonal or more diagonals matrices using the LU decomposition method, which greatly reduces the computational efficiency. Zhang et al. (2007) divided the spatial derivative operator into three directions of x, y, z and derived split implicit finite-difference scheme for second-order spatial derivatives. This method can not only give accurate point-source response, but also can be applied to the condition of velocity lateral variation. This algorithm is unconditionally stable if given the appropriate parameters. Zhou and Zhang (2011) proposed a prefactored strategy which divides the original diagonal matrix into two independent low triangular matrix and upper triangular matrix. This prefactored optimized implicit finite-difference scheme replaces the traditional LU decomposition method by using L+U prefactored method. Zhou and Zhang (2011) utilized Taylor series expansion and least-squares method to obtain optimized coefficients of implicit algorithm and improved the modeling accuracy significantly. Whereas these two methods are only suitable for second-order acoustic wave equation and cannot be extended to the condition of first-order elastic wave equation. Liu and Sen(2009a, b) derived new regular-grid and staggeredgrid implicit finite-difference operators for spatial derivatives of first-order and second-order based on the work of Claerbout (1985). Afterwards, Liu (2014) improved the computational accuracy and efficiency remarkably by using least-squares method to acquire the optimize coefficients of this finite-difference scheme.

In this paper, inspired by the idea of Claerbout (1985) and Liu and Sen(2009a, b), we propose a new kind of globally optimized implicit staggered-grid finite-difference (OISFD) scheme for first-order derivatives. Then we calculate the optimized coefficients by non-linear optimal method. Specifically, we firstly construct the objective function based on the two norm theory and then use simulated annealing method to solve optimized coefficients. Afterwards, we introduce this operator into the first-order velocity-stress elastic wave equation and test the numerical simulation results with the theoretical model. Meanwhile, we analyse the modeling accuracy and efficiency, and compare the optimized operator with the conventional explicit staggered-grid finite-difference method (ESFDM) and implicit staggered-grid finite-difference method (ISFDM).

2 FIRST-ORDER ELASTIC VELOCITY-STRESS EQUATION AND HIGH-ORDER STAGGERED-GRID FINITE-DIFFERENCE 2.1 First-Order Spatial Derivative Implicit Staggered-Grid Finite-Difference Scheme

Considering the first-order spatial derivatives of function u(x), its (2M + 2N)th-order ISFDM (function provides 2M order accuracy and the derivative of function provides 2N order accuracy) can be represented as follows (Chu, 2009; Kosloff et al., 2010):

(1)

where , M and N are positive integers, a_m and b_n are corresponding finite-difference coefficients, h is horizontal grid interval.

Using the plane wave theory, we let

(2)

where u₀ is a constant, , k is the wavenumber. Substituting formula (2) into formula (1), and setting β=kh/2, we have

(3)

Expanding the cosine and sine term of formula (3) by using Taylor series expansion, we can directly calculate the coefficients by matrix inversion. The calculation coefficients can be seen in Table 1.

2.2 First-Order Spatial Derivative Optimized Implicit Staggered-Grid Finite-Difference Scheme

From the derivation of ISFDM, we can see that 1D linear equation matrix inverse algorithm should be used in solving spatial derivative :

(4)

Thus

(5)

where matrix B_x, A_x is determined by a_m, b_n(m=1, …, M, n=1, 2, …, N) of formula (1) respectively, matrix size is (NX − 1) × (NX − 1) or (NX − 2) × (NX − 2), NX is the number of horizontal samples.

Formula (1) usually needs to meet the conditions of N > 1 in order to achieve the needed high precision. When N=1, we have

(6)

where B_x is a triangular matrix. We can solve it by chasing algorithm or LU decomposition algorithm. When N=2, we have

(7)

where B_x is a pentadiagonal matrix, thus chasing algorithm is not applicable under this situation. We must solve this linear system through LU decomposition algorithm.

We can notice that the chasing algorithm is not applicable under this situation, when N≥2. Therefore, LU decomposition algorithm becomes an optimal algorithm.

The computational efficiency of matrix inverse algorithm is a crucial factor for the implicit finite-difference modeling. Independent tests have shown that for same size matrix inversion the average CPU time of the chasing method is much less than that of LU decomposition method (Table 2). Therefore, constructing a new kind of implicit finite-difference algorithm which can be solved by the chasing method is an efficient way of improving the computational efficiency of the implicit finite-difference modeling.

According to Claerbout (1985) and Liu and Sen (2009b)’s idea, the second-order FD operator for function u(x) can improve its accuracy by adding one higherorder item. To be programed easily, the expression is suggested to be changed as follows (Claerbout 1985), where b is an adjustable constant.

(8)

Liu and Sen (2009b) proposed implicit staggered-grid finite-difference scheme which improves the computational accuracy by adding a third-order derivative item:

(9)

Inspired by this idea, we construct an OISFD operator expressed as

(10)

that is

(11)

where a, c are adjustable constants, h is the grid intervals. For simplicity, we define as follows:

(12)

(13)

We write the formula (12) as the form of 2M order when we need high-order accuracy:

(14)

we can rewrite formula (10) into the following matrix form:

(15)

that is

(16)

where

(17a)

(17b)

The size of matrix C_x and A_x is (NX−1)×(NX−1) or (NX−2)×(NX−2). Similarly, The size of matrix C_x and A_x is (NZ −1)×(NZ −1) or (NZ −2)×(NZ −2) in the z direction. NX and NZ is the horizontal and vertical grid numbers respectively. From formula (16) we can see that during the computation of the OISFD operator, we use the chasing method to obtain twice matrix inversions, thus avoid the time-consuming LU decomposition procedure and improve the computational efficiency.

2.3 Computing Strategy of Tridiagonal Matrix in Implicit Algorithm

The formula (16) requires twice tridiagonal matrix inversions, when solving one spatial derivative. Therefore we can calculate some of the repeated inverse process in advance, thus improving the OISFDM’s computing efficiency (Liu and Sen, 2009a, b). We can observe that both the tridiagonal matrix in formula (17a) and (17b) contain only one parameter, i.e. a or c. The value of each element on the diagonal is equal, and the same kind of situation applies equally to the two relative up and down diagonals. Therefore, we don’t need to establish the real matrix C_x and A_x, thus directly assigning the corresponding parameters during the calculation of chasing algorithm. What’s more, we can calculate the process that decomposes the matrix into the lower and the upper matrices and complete the matrix inversion using the chasing method before time iterative calculation. This means we can greatly save the memory and computational cost when using our optimized OISFD methods.

3 CALCULATE THE COEFFICIENTS OF THE OPTIMIZED IMPLICIT STAGGERED-GRID FINITE-DIFFERENCE FORMULAS BY SIMULATED ANNEALING ALGORITHM

Using the plane wave theory, we let u(x)=u₀e^ikx, β=kh/2 and substitute the formula (13) and (14) into formula (10), then we have

(18)

According to formula (18), we construct the object function of the 2M-th-order by means of the two norm theory, and get

(19)

where

(20)

And in this paper we choose β_l=0, β_r=0.58π.

We should note that the optimized operator defined in formula (11) has a multiplication of FD coefficients a and c. Therefore it is a nonlinear FD operator, and cannot be obtained by conventional method in solving linear matrix inversion. Commonly used non-linear global optimization methods are conjugate gradient method, gauss-newton method, simulated annealing method, etc. So far, the simulated annealing method has been widely recognized, because of its ability to overcome local minimum, the merit to overcome the dependence of initial value, the advantage of constructing flexible objective function. Thus in this paper we choose the simulated annealing algorithm and acquire the optimized coefficients in formula (19).

The idea of simulated annealing algorithm was put forward by Metropolis etc. (1953). Fig. 1 is the flow chart of the ISFD scheme using the simulated annealing algorithm. The algorithm is trying to find a set of optimal coefficients under a given Tolerance limit. Function obj is defined by the formula (19) and (20). The initial-values of the parameters in the algorithm are as follows: the length of MARKOV chain is N=1000, the initial temperature is T=1000, the final temperature is T₀=100, the step factor is Stepfactor=0.02, the tolerance is Tolerance=5×10⁻⁸, the maximum search value is XMAX=[1; 1], PreX=XMAX×rand, PreBest X=PreX, Best X=–XMAX×rand, select the next point randomly NextX=PreX+Stepfactor×XMAX×(rand–0.5).

According to the principle of thermodynamics, at the temperature of T, the cooling probability of the energy difference dE is P (dE), and it could be expressed as

(21)

where

(22)

K is the Boltzmann constant in physics; exp is the natural index, and dE > 0. The higher the temperature is, the greater the probability of reaching an energy difference of dE, and vice versa. Essentially, temperature T controls the perturbation range of the optimal solution. High temperature is more likely to make it achieve a wider range (Zhang and Yao, 2013). Because −dE/KT < 0, the range of function P (dE) is (0, 1). In fact, under the given tolerance limit Tolerrance, there will be many different groups of solutions whereas other optimization methods have only one set of solution. Therefore, we can take a further step to choose a set of optimal solutions by balancing between the error and the accurate wave number coverage (Zhang and Yao, 2013).

Based on the assumption of discretization, we use the optimized first-order operator to approach the real wave number, and get the coefficients array X=[a; c]. Through the calculation process above, we get the array X=[0.09737101; −0.05866972]. In other words, a=0.09737101, c=−0.05866972, while M=1. The higher order difference coefficients are shown in Table 3. In fact, due to the flexibility of simulated annealing method, we can also use the maximum norm or the relative error method for setting the objective function (Zhang and Yao, 2013).

4 ACCURACY COMPARISONS BETWEEN THE IMPLICIT SCHEME AND THE EXPLICIT SCHEME

Conventional ESFDM accuracy analysis can be expressed as (Liu and Sen, 2009b)

(23)

where N is a positive integer, cn is the corresponding differential coefficient, the values are shown in Table 1 from Liu and Sen (2009b).

According to formula (1), we can acquire the conventional ESFDM accuracy analysis expression:

(24)

where M and N are positive integers, a_m and b_n are the corresponding differential coefficients shown in Table 1. The ISFDM proposed by Liu and Sen (2009b) could be described as

(25)

where N is a positive integer, c_n is the corresponding differential coefficient, b is a constant, and their values can be found in Table 2 from Liu and Sen (2009b). The ISFDM proposed by Liu and Sen (2014) uses the same expression as formula (25), and the coefficients are calculated using the least squares method, the values are listed in Table 9 from Liu and Sen (2014).

Here we define fOISFDM(β) as follows to describe the numerical dispersion of OISFD operator:

(26)

where a and c are constants, c_m is the corresponding differential coefficient shown in Table 3.

According to the accuracy curves of approximation from Fig. 2 to Fig. 5, we can indicate that the OISFD operators have wider frequency coverage and smaller accuracy error fluctuation than those of the conventional ESFD operators and ISFD operators when using the same order spatial difference operator. What’s more, the curves demonstrate that the accuracy of the (2+4)th-order OISFD operator is much higher than those of the conventional 6th-order and the 8th-order ESFD operators, and it can even approach the accuracy of conventional 10th-order ESFD operator. In contrast with the conventional ISFD operator, the accuracy of the (2+4)th-order OISFD operator is between the conventional (4+2)th-order and the (6+2)th-order ISFD operators. In addition, the (2+4)th-order OISFD operator’s accuracy is higher than Liu and Sen (2009 b)’s (4+2)th-order ISFD operator, lower than Liu and Sen (2009 b)’s (6+2)th-order ISFD operator, and close to the accuracy of Liu (2014)’s (4+2)th-order ISFD operator. With the increase of difference order, the accuracy of OISFD operator got remarkable improvement.

5 NUMERICAL SIMULATION AND ANALYSIS

First-order velocity-stress wave equation for 2D elastic media can be represented as follows (Virieux et al., 1986):

(27)

where t is time, x, z is coordinate direction, (v_x, v_z) is velocity vector, (τ_xx, τ_zz, τ_xz) is a stress vector containing three components, λ, µ is Lame constant, ρ is the density.

6 NUMERICAL MODELING EXAMPLES 6.1 Homogeneous Model

To demonstrate the ability of OISFDM in elastic wave modelling, we first design a homogeneous model with the parameters listed in Table 4, and test the results of the conventional ESFDM, the conventional ISFDM, Liu and Sen (2009b)’s ISFDM, and OISFDM of numerical simulation.

Figure 6 is the snapshots computed by different methods respectively, at the moment of 0.55 s. It demonstrates that the dispersion of FD modelling decreases with the increase of order. The comparison between snapshots shows that the conventional 6th-order ESFDM leads to a strongly distorted waveform, while the dispersion caused by the (2+4)th-order OISFDM is quite weak and could achieve nearly identical and sometimes even higher accuracy than that of the conventional 10th-order ESFDM under the same modeling parameters. Fig. 6 also shows that (2+4)th-order OISFDM has better modelling waveform than that of the conventional (4+2)th-order ISFDM, and its simulation accuracy is similar with the conventional (4+4)th-order ISFDM, between those of the Liu and Sen (2009b)’s (4+2)th-order ISFDM and (6+2)th-order ISFDM.

6.2 Complex Sigsbee2A Model

To further examine the advantage of the OISFDM, we also perform numerical modelling on the modified 2D Sigsbee2A model in Fig. 7, and Table 5 gives modeling parameters used in this simulation. Fig. 8 and Fig. 10 show the seismograms of horizontal and vertical component computed by the OISFDM, the conventional ESFDM and the conventional ISFDM, respectively. In order to further compare the numerical modeling results clearly, we make local magnification of Fig. 8 and Fig. 10, and get the corresponding Fig. 9 and Fig. 11. Note that the result obtained by the (2+4)th-order OISFDM can be more effective in suppressing the numerical dispersion compared with the results of conventional 6th-order ESFDM and (4+2)th-order ISFDM, and can lead to almost the same or better modelling accuracy than that of the conventional (4+4)th-order ISFDM, thus between the accuracy of Liu and Sen (2009b)’s (4+2)thorder and (6+2)th-order ISFDM.

Figure 12 gives the corresponding snapshots at t=0.92 s. We observe that the (2+4)th-order OISFDM can obviously depress the waveform distortion, and it could achieve almost the same simulation accuracy as the 10th-order ESFDM. The conventional (4+2)th-order and even the (4+4)th-order ISFDM can cause greater numerical dispersion than that of the (2+4)th-order OISFDM. When compared with Liu and Sen (2009b)’s ISFDM, the accuracy of the (2+4)th-order OISFDM is between those of Liu and Sen (2009b)’s (4+2)th-order and (6+2)th-order ISFDM. Therefore, the results are fully consistent with the theoretical accuracy comparison in the previous section.

Additionally, according to the computation costs of 4000 time iterations in Fig. 13, the (2+4)th-order OISFDM requires the least CPU time, and thus achieves the highest computational efficiency. What’s more, according to the formula (1), as the value N increases, the simulation accuracy improvement of the conventional ISFDM is limited, under the condition of the same value M. It couldn’t significantly improve the calculation accuracy, because of the increase in the inverse matrix bandwidth, which will significantly add the calculation cost and reduce the computational efficiency.

7 CONCLUSIONS

We have presented an OISFD operator with any order of accuracy for calculating first spatial derivatives on a staggered-grid, and developed the objective function by means of the two norm theory. Different from traditional global optimized methods, we select the simulated annealing algorithm to get the optimized coefficients. Finally, after introducing the optimized scheme to the homogeneous and complex medium and conducting elastic wave-field simulation, we can indicate by the seismograms, snapshots and computational efficiency that the OISFD operators can achieve higher accuracy comparing with those using the conventional ESFD operators, ISFD operators and Liu and Sen (2009b)’s ISFD operators. On the premise of the same simulation parameters, the (2+4)th-order OISFDM can be more effective than conventional 6th-order ESFDM in suppressing numerical dispersion and can achieve the same precision as the conventional 10th-order ESFDM; in contrast with the conventional ISFDM, we can see that the (2+4)th-order OISFDM can provide better numerical results than conventional (4+2)th-order ISFDM, and can obtain slightly higher precision than conventional (4+4)th-order ISFDM; furthermore, the (2+4)th-order OISFDM can reach the simulation accuracy between that of Liu and Sen (2009b)’s (4+2)th-order ISFDM and (6+2)th-order ISFDM. Therefore, the OISFDM has higher simulation accuracy and computational efficiency, thus providing a new numerical algorithm for the elastic wave forward modeling.

In addition, the computation amount and memory requirement will have a sharp increase, when using the OISFDM to the 3D complex medium, which greatly affects the efficiency of the numerical simulation. So we can consider introducing the GPU/CPU collaborative parallel architecture, divide the model in the slow block in parallel direction, by the GPU in each wave field numerical calculation of each time step. This means we can greatly save the memory demand and computational cost when using our optimized OISFD methods. But this idea still needs to be discussed in detail in further step.

ACKNOWLEDGMENTS

We acknowledge Doctor Hongyong Yan for his guidance and help for this paper. This work was supported by the National 863 Program of China (2012AA061202), the National Science and Technology of Major Projects (2011ZX05008-006-50, 2011ZX05003-003, 2011ZX05023-005-002).