## ABSTRACT

The staggered-grid finite-difference (FD) method is widely used in numerical simulation of the wave equation. With stability conditions, grid dispersion often exists because of the discretization of the time and the spatial derivatives in the wave equation. Therefore, suppressing grid dispersion is a key problem for the staggered-grid FD schemes. To reduce the grid dispersion, the traditional method uses high-order staggered-grid schemes in the space domain. However, the wave is propagated in the time and space domain simultaneously. Therefore, some researchers proposed to derive staggered-grid FD schemes based on the time-space domain dispersion relationship. However, such methods were restricted to low frequencies and special angles of propagation. We have developed a regularizing technique to tackle the ill-conditioned property of the symmetric linear system and to stably provide approximate solutions of the FD coefficients for acoustic-wave equations. Dispersion analysis and seismic numerical simulations determined that the proposed method satisfies the dispersion relationship over a much wider range of frequencies and angles of propagation and can ensure FD coefficients being solved via a well-posed linear system and hence improve the forward modeling precision.

## INTRODUCTION

Finite-difference (FD) methods for acoustic-wave equations are widely applied in seismic wave simulations because of the efficiency of computation, over memory requirements, and simplicity of realization. These methods also constitute the basis for reverse-time migration and full-waveform inversion (Alford et al., 1974; Kelley et al., 1976; Alkhalifah, 2000; Basabe and Sen, 2007; Yang et al., 2012). Meanwhile, grid dispersion is one of the key numerical problems when using FD methods. An FD scheme dominated by spatial dispersion delays higher frequencies, whereas an FD scheme dominated by temporal dispersion advances higher frequencies (Dablain, 1986).

To reduce the grid dispersion, two methods could be applied: the traditional method using high-order staggered-grid schemes in the spatial domain and the method of using shorter time and spacing grid intervals. Higher order approximations of spatial derivatives and temporal derivatives are usually applied to increase accuracy and reduce grid dispersion (Chen, 2007, 2011; Chu and Stoffa, 2012; Liang et al., 2013b). Generally, the spatial FD coefficients are determined only in the spatial domain. However, wave equations are solved in the temporal and spatial domains simultaneously (Finkelstein and Kastner, 2007, 2008; Liu and Sen, 2009, 2011; Liang et al., 2013a). Finkelstein and Kastner (2007, 2008) propose a systematic design methodology for obtaining FD coefficients to reduce dispersion, which allows the exact phase velocity or (and) group velocity dispersion relationship to be satisfied at some designated frequencies in the temporal-spatial domain. Liu and Sen (2009) propose a new time-space-domain method to determine the higher order FD coefficients for 1D, 2D, and 3D wave equations, and then they use this method to get variable length FD coefficients (Liu and Sen, 2011). Etgen (2007) proposes minimization of the phase-velocity error in the range of frequencies and prorogation-angles of interest, based on a weighted least-squares method. Gauss-Newton solution method are also considered in Etgen (2007). Zhang and Yao (2013) propose using the simulated annealing algorithm and give an error limitation for determining the FD coefficients in space or time-space domain. Liu et al. (2014) propose an explicit time-evolution method to simulate wave propagation in acoustic media with high temporal accuracy. They find that in constant density acoustic media, including slightly more stencil points, could significantly reduce the inaccuracy at low wavenumbers.

Advantages of the staggered-grid FD methods are that they possess greater accuracy and better stability than traditional FD methods. Many results have been achieved using the staggered-grid FD scheme of spatial partial derivatives thus far. Liu and Sen (2011) introduce the time-space domain method into the staggered-grid FD simulations of the acoustic-wave equations. They use Taylor expansion in the wavenumber direction to establish a system of linear equations and to determine the FD coefficients in a fixed angle of wave propagation. Their method could indeed improve the forward modeling precision.

Numerical experiments of the time-space staggered-grid FD method reveal that (1) Taylor expansion in the wavenumber direction restricts the dispersion relationship to be preserved only in limited wavenumbers, (2) determining the FD coefficients using the dispersion relation in a particular wave-propagation direction may lead to loss of precision in other directions and may also induce numerical anisotropy, and (3) the linear system is symmetric and ill-conditioned. In this paper, we propose a new method to satisfy the dispersion relationship in a much wider range of frequencies and angles of propagation. We also propose a regularizing technique to tackle the ill-conditioned property of the symmetric linear system and to stably provide approximate solutions of the FD coefficients.

## FINITE-DIFFERENCE METHODS USING STAGGERED-GRID SCHEMES IN THE TIME-SPACE DOMAIN

### Staggered-grid finite-difference operators

The time-space domain for the staggered-grid FD scheme can be obtained by substituting $\omega =kv$ into equation 6 and taking the Taylor expansion in $kh$ and assuming a constant angle of wave propagation $\theta =0$ and $\varphi =\pi /8$ (Liu and Sen, 2011).

### Optimizing regularized staggered-grid finite-difference operators

### Regularizing computation of the coefficients

## NUMERICAL DISPERSION-ERROR ANALYSIS

### 1D dispersion analysis

### 2D dispersion analysis

A comparison of the dispersion errors of the FD operator using traditional FD scheme, Liu and Sen’s FD scheme, and our new FD scheme for the 2D FD operator is shown in Figures 4 and 5, for $v=1500\u2009\u2009m/s$ and $v=4500\u2009\u2009m/s$, and the same $M$, $h$, and $\tau $, respectively. It indicates from Figures 4 and 5 that Liu and Sen’s FD method and our new FD method perform better than the traditional FD method. In addition, with increasing of the length of the FD operator, the dispersion relationship preserving frequency range of Liu and Sen’s FD method increases slowly, whereas the new method can preserve the dispersion relation in a larger frequency range.

We remark that when $\alpha =0$, our method reduces to the well-known least-squares method. For short length of the FD operator, our method performs the same as the least-squares method. However, when the length of the FD operator is large, our method will be much more stable than the least-squares method to define the stencil coefficients. The same reason is that the new method possesses the regularizing property when dealing with ill-conditioned linear system and hence provides a stable solution.

### 3D dispersion analysis

Comparison of the dispersion errors of the FD operator using traditional FD scheme and our new FD scheme for the 3D FD operator is shown in Figure 6a and 6b. It indicates from Figure 6 that our new FD method performs better than the traditional FD method and can preserve the dispersion relation in a larger frequency range.

We remark that during our simulations, when the length of the FD operator is long, the least-squares method cannot yield a reasonable solution; instead, using our regularizing least-squares method, the solution is always satisfactory. The reason lies in that the matrix $C$ is extremely ill-conditioned; for example, when $M=20$, the conditioning number can reach $O(1017)$, which is a huge number and may lead to unstable solutions for least-squares method.

## NUMERICAL SIMULATIONS

Because the new method (i.e., the regularizing least-squares method) preserves the numerical dispersion relationship better than the least-squares method, in the following simulations, we only list the comparison results using the traditional method, Liu and Sen’s method, and our new method.

### A homogeneous model

Figure 8a and 8b plots slices of the snapshots of the wavefield. It shows details of differences of the traditional method, the pseudospectrum method, and the new FD method at 1530 ms. Figure 8a is a plot of slices of snapshots when the time step equals 1.5 ms. Comparison of the results from Figure 8a reveals that the pseudospectrum and the new dispersion-relationship-preserving methods can reduce the dispersion and hence provides more accurate wavefield simulation. Figure 8b illustrates slices of snapshots when the time step equals 4.5 ms. Meanwhile, we regard the analytic solution is from the pseudospectrum method with small time step (1.5 ms). It is observed that our new time-space-domain method is effective to reduce temporal dispersion when the time step becomes large.

### Salt model

In the following, we present simulation results using a widely referred salt velocity model from SEG. The velocity model is shown in Figure 9 with variations of velocities from $1486$ to $4790\u2009\u2009m/s$. The spatial sampling interval equals 20 m, the temporal step is 1 ms, and $M=7$ for all FD operators. With the same source function as in former examples and $f0=45\u2009\u2009Hz$, the simulation results using the traditional FD method, Liu and Sen’s method, and the new method can be obtained. To show details of their difference, we also plot slices of the snapshots of the wavefield at 2500 ms ($z=620\u2009\u2009m$) in Figure 10a–10c. It is evident that the new method provides the best simulation results with least dispersion.

## CONCLUSION

Based on analysis of the determination of the FD coefficients in temporal-spatial domain using staggered-grid schemes, we proposed a new method for determining the FD coefficients that has the following features: (1) the method specifies the wavenumber (frequency) upper limit according to the source frequency, the spatial grid interval, and the velocity, whereas previous methods use the same upper limit of the frequency for different velocities; (2) the method considers different angles of the wave propagation, whereas the previous methods use a fixed angle; and (3) a regularization model by perturbation is proposed to ensure FD coefficients being solved via a well-posed linear system.

With the dispersion analysis of our new method and comparison with previous FD coefficients determination methods using staggered-grid schemes, we conclude that our method is feasible for seismic wave modeling. As a result, our method can be a substitute for the traditional staggered-grid FD coefficients determination methods, whereas these methods are essential in forward seismic wave modeling and reverse-time migration.

## ACKNOWLEDGMENTS

We deeply thank the reviewers and the associate editors for their very helpful comments that improved the manuscript. This work is supported by the National Natural Science Foundation of China under grant nos. 41325016 and 11271349 and R&D of Key Instruments and Technologies for Deep Resources Prospecting (the National R&D Projects for Key Scientific Instruments) under grant no. ZDYZ2012-1-02-04.