1 INTRODUCTION

With the increase of exploration depth of oil-producing region, there are increasing needs for the higher precision imaging and inversion, and then as a basis of imaging and inversion methods, numerical modeling technology becomes more and more important. Numerical modeling method based on wave equation contains the kinematics and dynamics characteristics of wavefields, and is widely used(Cheng et al., 2008). The finitedifference method(Alterman and Karal, 1968; Kelly et al., 1976; Dablain, 1986; Igel et al., 1995) and pseudospectral method(Gazdag, 1981; Kosloff and Baysal, 1982; Carcione, 1992)as two numerical modeling method based on wave equation are well known. The finite-difference method was put forward in the late 1960s, Alterman and Karal(1968)first used explicit 2nd order finite-difference operator to simulate the elastic wavefields in layered media. Kelly et al.(1976)proposed the finite-difference method which was suitable for heterogeneous media. Dablain(1986)carried out scalar wave numerical modeling using high-order operators. Liu et al.(1998)derived the finite-difference operators of any even-order accuracy for any even derivative from the Taylor series expansion. Dong et al.(2000)derived a staggered-grid high-order finite-difference method for first-order elastic wave equation. Pei and Mou(2003)derived staggered-grid finite-difference operators of any even-order accuracy. Liu and Sen(2009a, 2010)proposed a new time-space domain high-order finite-difference method for the acoustic wave equation, and then Yan and Liu(2013a, 2013b)extended their methods to visco-acoustic and anisotropic acoustic equation.

Numerical dispersion of the finite-difference method is inevitable because of the approximation of difference operator to differential operator. In order to decrease the numerical dispersion, we have to adopt a finer grid or lower wavelet frequency, however, which will respectively bring huge computational cost and the loss of part of the high frequencies. Some geophysicists have proposed a variety of methods to suppress the numerical dispersion. The first method is FCT technology. Boris and Book(1973)first applied the FCT in reducing the numerical dispersion. Yang and Teng(1997)extended this method to three-component seismic numerical modeling for anisotropic media. The second method is optimization methods. Holberg(1987) and Robertsson(1994)minimized the error using optimization methods to calculate the finite-difference coefficients. Zhang and Yao(2012)used the simulated annealing algorithm to minimize the objective functions constructed by the maximum norm to get the finite-difference coefficients, which has larger wave number coverage and smaller margin of error. Liu(2013)adopted least-squares method to get global optimized finite-difference operators. Yang et al.(2014)extended this method to staggered-grid finite-difference methods. The third method is truncated windows methods. The methods are based on the knowledge that the pseudo-spectral method can be interpreted as the limit of finite-difference schemes of increasing orders(Fornberg, 1987). Zhou and Greenhalgh(1992)chose the generalized powered Hanning window to get optimized finite-difference operators. Igel et al.(1995)used Gaussian window to develop staggered-grid finite-difference operators. Chu and Stoffa(2012)unified the finite-difference method and the pseudo-spectral method using the binomial window, and adopted an improved binomial window to get optimized finite-difference operators, and they thought that this strategy of improved binomial window is the same as that of truncated finite-difference methods proposed by Liu and Sen(2009b)in essence.

In space domain, we can truncate the spatial convolution series of the pseudo-spectral method to get finitedifference method. In terms of digital signal processor, truncating is equal to using a window function which will cause spectrum leakage. In order to suppress the numerical dispersion, choosing some suitable window function is necessary, different window functions have different results. Specifically, the window function which has narrower main lobe and larger attenuation of side lobe in its amplitude response bring higher accuracy of finite-difference approximation, in this paper, we design a window function based on this theory to truncate the spatial convolution series of the pseudo-spectral method to get optimized finite-difference operators which can efficiently suppress the numerical dispersion.

The generalized weighted Hanning window designed by Zhou and Greenhalgh(1992)is one of trigonometric window functions, which has relatively small attenuation of side lobe in amplitude response thus leads to the fluctuation of accuracy error. The fluctuation of accuracy error greatly affects the result of numerical modelling(Zhang and Yao, 2012). The improved binomial window proposed by Chu and Stoffa(2012)has a wide main lobe in its amplitude response, resulting in a narrow wave number coverage range, at the same time, the accuracy error curves of finite-difference operators based on improved binomial window show obvious fluctuation especially for the low-order operators(below twelfth-order). In this paper, we first analyze the influence on the accuracy of finite-difference operators caused by the properties of main lobe and side lobe in the amplitude response of truncated windows, and then introduce the auto-convolution and weighted combination to the design of window functions. Finally, we design Chebyshev auto-convolution combined window which has narrow main lobe and large attenuation of side lobe in its amplitude response. Through adjusting three parameters of auto-convolution combined window, visually and straight forwardly control the properties of main lobe and side lobe to adjust the accuracy of finite-difference approximation.

2 APPROXIMATION ERROR ANALYSIS 2.1 Conventional Finite-Difference Operators

According to the discrete signal sampling theorem, a b and limited continuous signal f(x)can be reconstructed by a uniform sampling signal f_n using the sinc interpolator(Diniz et al., 2012):

where △x is the spatial sampling interval. $\frac{\pi }{\Delta x}$ is Nyquist wave number.

We separately compute the first derivative and the second derivative at x = 0 by differentiating both sides of Eq.(1):

We can truncate the Eqs.(2) and (3)using an(N +1)-point window function to get conventional finite-difference operators, where N is an even number. In Eqs.(4) and (5), w_n^N = $\left( \begin{matrix} N \\ \frac{N}{2}+n \\ \end{matrix} \right)/\left( \begin{matrix} N \\ \frac{N}{2} \\ \end{matrix} \right)$ and $\left( \frac{A}{B} \right)=\frac{A!}{B!\left( A-B \right)!}$. Chu and Stoffa(2012)named this window function as binomial window, and then proposed an improved binomial window w_n^N+M to optimize the finite-difference operators, where M is positive even number. For conventional finite-difference operators, M = 0.

2.2 Influence of the Truncated Windows on the Accuracy

Zhou and Greenhalgh(1992)designed the generalized weighted Hanning window which is one of trigonometric window functions. Igel et al.(1995)used Gaussian window. Chu and Stoffa(2012)proposed the improved binomial window. We compare these four window functions, analyze the influence on the accuracy of finitedifference operators caused by the properties of main lobe and side lobe in the amplitude response of window functions.

Figure 1 and Fig. 2 respectively display these four window functions for N = 24 and their amplitude responses. As shown in Fig. 1 and Fig. 2, the improved binomial window applies more tapering than other windows, accordingly, has wider main lobe in its amplitude response. Diniz et al.(2012)thought that wide main lobe means wide transition b and which reduces the wavenumber coverage and causes numerical dispersion for the part of high wavenumber. Compared with other windows, the improved binomial window has widest main lobe, as a result, the finite-difference operators truncated by the improved binomial window are not accurate for the part of high wavenumber, meanwhile, the improved binomial window has largest attenuation of side lobe, which bring best stability of accuracy error. As for Kaiser and Chebyshev window functions, we can adjust the parameters β and r which are respectively the parameters of these two window functions to get different amplitude responses. Here, we set r = 60 and β = 5.

Figure 3 shows the amplitude responses of Kaiser window function(β = 2, 5, 9) and Chebyshev window function(r = 35, 60, 95)for N = 24.

As demonstrated in Fig. 3, for Kaiser window function, with the increase of β, main lobe becomes wider and attenuation of side lobe becomes larger; for Chebyshev window function, with the increase of r, the variation is the same. Compared with these two window functions, Chebyshev window function has relatively larger attenuation of side lobe, and the finite-difference operators truncated by Chebyshev window function have better stability of accuracy error.

n = 0 is singular point. Lee and Seo(2002)think that Eq.(2) and Eq.(3)can be expressed as follows:

where d_n¹ = $-\frac{1}{n}$ cos(nπ), d_n² = $-\frac{2}{{{n}^{2}}}$ cos(nπ), n = 1, 2, 3, · · ·, ∞, d₀² = $-\sum\limits_{n=1}^{\infty }{\left( d_{n}^{2}+d_{-n}^{2} \right)}=4\times \left[ -\varsigma \left( 2 \right)/2 \right]=-\frac{{{\pi }^{2}}}{3}$, $\varsigma $ is Riemann $\varsigma $ function.

We can get finite-difference operators by truncating Eqs.(6) and (7)using window functions, which can be generalized to the Eqs.(8) and (9):

where b_n = d_n¹ w(n)= - $\frac{1}{n}$ cos(nπ)w(n), c_n = d_n² w(n)= - $\frac{2}{{{n}^{2}}}$ cos(nπ)w(n), n = 1, 2, · · ·, N/2, c₀ = -2 $\sum\limits_{n=1}^{N/2}{\left( {{c}_{n}}+{{c}_{n}} \right)}$. We perform the Fourier transform on Eqs.(8) and (9), which give We use the Eq.(12)to measure the accuracy errors of the finite-difference operators for the second derivative:

Figure 4a shows that, accuracy error curves for second derivative of Hanning and Kaiser window functions have wider wavenumber coverage. The improved binomial window brings the narrowest wavenumber coverage. The finite-difference operators truncated by Chebyshev window function has appropriate wave number coverage. We magnify the accuracy error curves one thous and times; the fluctuation of accuracy error curves of Kaiser and Hanning windows is more obvious than that of the remaining two windows. Considering the wave number coverage and the fluctuation of accuracy error curves, Chebyshev window is a better choice to truncate Eqs.(6) and (7)to derive optimized finite-difference operators for elastic wave modelling.

2.3 Accuracy Error Analysis

We apply one auto-convolution on an(N + 1)-point Chebyshev window function to get an(2N + 1)-point Chebyshev auto-convolution window function, N is an even number. Similarly, we can perform L autoconvolutions on Chebyshev window function to get L-order Chebyshev auto-convolution window function. We first perform the Fourier transform on Chebyshev window function, ω = k_x△x, which gives

The convolution operation in space domain is equal to multiplication operation in frequency domain. As a result, L-order Chebyshev auto-convolution window function can be gotten from multiplication of LW(ω)in frequency domain.

Amplitude responses of Chebyshev auto-convolution window functions after different number of autoconvolutions are shown in Fig. 5. Auto-convolutions obviously increase the attenuation of side lobe and make the main lobe wider at the same time, with the increase of number of auto-convolutions, the impact is greater.

We fix the length of window functions, here N = 24. Two and four auto-convolutions are separately flagged as L = 2 and L = 4. The amplitude responses of Chebyshev auto-convolution window functions after two and four auto-convolutions are shown in Fig. 6, with the increase of number of auto-convolutions, main lobe becomes wider and attenuation of side lobe larger.

We still utilize the Eq.(12)to measure the accuracy errors of the finite-difference operators for the second derivative. Fig. 7 displays the accuracy error curves of the finite-difference operators truncated by three window functions shown in Fig. 6 for N = 24.

As shown in Fig. 7, with the increase of the number of auto-convolutions, the finite-difference operators truncated by the auto-convolution window functions have narrower wavenumber coverage, however, have better stability of accuracy error. Considering the wavenumber coverage and stability of accuracy error, we can get window function which has good properties of truncation through performing two auto-convolutions.

Based on the previous analysis, our purpose is to design a window function which has narrow main lobe and large attenuation of side lobe, and get the finite-difference operators which have better stability of accuracy error on the premise of maintaining appropriate wave number coverage, thus we introduce a linear weighted function:

where, w_c(n)is the Chebyshev window function, w_c^L (n)is the Chebyshev window function after L autoconvolutions, w_a(n)is combined window function. λ is weight coefficient.

We can choose different auto-convolution combined window functions through adjusting three parameters to get the optimized finite-difference operators. The parameters are ripple rate r, the number of autoconvolutions L and weight coefficient λ. In this paper, we set λ = 0.89, r = 60, L = 2. Fig. 8 displays the accuracy error curves of the finite-difference operators(N=8, 12, 16, 24)truncated by the Chebyshev auto-convolution combined window function for the second derivative. The finite-difference operators based on Chebyshev autoconvolution combined window have much higher accuracy than the conventional finite-difference operators. The accuracy of our optimized eighth-order and twelfth-order finite-difference operators can respectively reach, even exceed that of the conventional twelfth-order and twenty-fourth-order finite-difference operators. The maximum deviation of absolute accuracy error curves is within a range of [0, 0.0004].

Figure 9 displays the comparison between the finite-difference operators based on the improved binomial window and the operators based on Chebyshev auto-convolution combined window function. As see from Fig. 9, apparently, for low-order operators(especially eighth-order in Fig. 9), the accuracy error curve of the operator based on the improved binomial window has obvious fluctuation, and for high-order operators, in the condition of same order, the accuracy of operators based on the improved binomial window is lower than that of the operators based on Chebyshev auto-convolution combined window function.

Finite-difference weights truncated by the Chebyshev auto-convolution combined window function for the second derivative with N=4, 8, 12, 16 are given in Table 1.

Similarly, we also use Eq.(12)to measure the accuracy errors of the finite-difference operators for the first derivative:

Figure 10 displays the accuracy error curves of the finite-difference operators(N=8, 12, 16, 24)truncated by the Chebyshev auto-convolution combined window function for the first derivative. As shown in Fig. 10, compared with the conventional finite-difference operators, the finite-difference operators based on Chebyshev auto-convolution combined window function have wider wavenumber coverage and better statability of accuracy error for the same order.

Finite-difference weights truncated by the Chebyshev auto-convolution combined window function for the first derivative with N=4, 8, 12, 16 are given in Table 2.

3 NUMERICAL EXAMPLES

The linear theory of elasticity establishes a relationship between vector wavefield displacement and the dimensionless linear strain tensor:

where ε_kl stands for linear strain tensor, u_l and u_k st and for vector wavefield displacements.

The linear strain tensor ε_kl can be related to the stress tensor σ_ij using a constitutive relationship:

where c_ijkl is stiffness tensor.

Substituting Eqs.(16) and (17)into the equations of motion derived from Newton’s second law, that describe wave propagation for anisotropic elastic media, ignoring the body force per unit volume:

In Eq.(18), ∂_tt² is the second-order temporal derivative.

3.1 Homogeneous Isotropic Media

A 2D homogeneous medium is defined on a grid of 255×300 with grid spacing of 5 m. The P-wave velocity is 2000 m·s^-1, The S-wave velocity is 1500 m·s^-1, a point source is located at the center of the media. Fig. 11 displays the snapshots of impulse responses generated by the finite-difference operators(N=8)truncated by the Chebyshev auto-convolution combined window function and the conventional finite-difference operators(N=8, 12)at nt=500, the dominant frequency of a Ricker wavelet is 30 Hz.

Figure 11a is the snapshots(X and Z components)generated by the eighth-order conventional operator, which have obvious numerical dispersion. The eighth-order operator based on Chebshev auto-convolution combined window function can effectively suppress numerical dispersion as shown in Fig. 11b. Using the twelfthorder conventional operator can also significantly suppresses the dispersions(Fig. 11c)but it doubles the computational cost. In conclusion, Compared with the conventional eighth-order and twelfth-order finite-difference operators, our optimized eighth-order finite-difference operator based on Chebyshev auto-convolution combined window function has a better result.

3.2 Marmousi Model

To further examine the quality of our optimized finite-difference operators based on combined window, we test on a Marmousi model. The grid is 737×750, △x=12.5 m, △z=4 m, the Marmousi P-wave velocity model is shown in Fig. 12, the S-wave velocity is a scaled version of the P-wave velocity with v_S = v_P/$\sqrt{3}$, the density is 1000 kg·m^-3, the dominant frequency of a Ricker wavelet is 30 Hz, positioned at x=231.25 m, z=25 m. The record-length is 5 s with a time step of 0.0005 s.

We ran the code on a PC with a Intel(R)Core i5 2300 @2.8GHz CPU and a 192-core NVIDIA GTX550Ti GPU card. The version of CUDA is 4.2. We position the source at x=0.0 m, 2303.1 m, 4606.3 m, 6909.4 m, 9212.5 m.

Computation cost using the conventional finite-difference operators and the finite-difference operators based on Chebyshev auto-convolution window function for different orders are given in Table 3 and Table 4.

Figure 13 shows the comparison of a shot record for the Marmousi model computed by the conventional finite-difference operators(N=8, 12) and the finite-difference operators based on Chebyshev auto-convolution combined window function(N=8)(Z component). Figs. 13(b, c, d)is separately the zoomed in view of block 1, 2, 3 of Fig. 13a.

The left figures of Figs. 13(b, c, d)generated by the eighth-order conventional operator have obvious numerical dispersion, and the middle figures of Figs. 13(b, c, d)generated by the eighth-order operators based on Chebshev auto-convolution combined window function have a better result, significantly suppress the dispersion. Using the twelfth-order conventional operator can also effectively suppresses the numerical dispersion(the right figures of Figs. 13(b, c, d)), however, compared with the shot record generated by the eighth-order operators based on Chebshev auto-convolution combined window function, it still has relatively obvious dispersion and doubles the computational cost. The numerical modelling results in the complex medium show that, the accuracy of the finite-difference operators based on Chebshev auto-convolution combined window function is better than that of the conventional finite-difference operators.

4 CONCLUSION

We propose optimized finite-difference operators based on Chebyshev auto-convolution combined window function. We first analyze the influence on the accuracy of finite-difference approximation caused by the properties of amplitude responses of window function, narrower main lobe and larger attenuation of side lobe bring higher accuracy of finite-difference approximation, and significantly suppress numerical dispersion. Compared with the generalized powered Hanning window and improved binomial window, the Chebyshev auto-convolution combined window has better attributes of amplitude response to get optimized finite-difference operators, and can visually and more intuitively adjust the accuracy of finite-difference approximation through adjusting three parameters.

The finite-difference operators based on Chebyshev auto-convolution combined window have higher accuracy than that of the conventional finite-difference operators. In this paper, we set λ=0.89, r=60, L=2. The eighth-order and twelfth-order operators based on Chebyshev auto-convolution combined window can respectively reach, even exceed the accuracy of the twelfth-order and twenty-fourth-order conventional finite-difference operators, and the maximum deviation of absolute error is within [0, 0.0004]. For higher-order operators, the accuracy increase becomes more obvious. Numerical dispersion analysis and elastic wave numerical modelling demonstrate that our finite-difference operators have higher accuracy than the conventional finite-difference operators and can more efficiently suppress numerical dispersion under the same condition.

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