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.

acoustic, finite difference

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.

We consider the 3D acoustic-wave equation:
(1)
where ρ denotes the density, p represents the pressure, K is the bulk modulus such that K=λ+2μ=ρv2 with v being the wave velocity, λ and μ are the Lame’s constants. With the source signature, equation 1 can be written as
(2)
where src(t) is the seismic source. Below, we only consider solving methods for simulation of the wavefields with constant density.

Staggered-grid finite-difference operators

Using the same FD coefficients for the x-, y-, and z-directions and taking the downward direction of z as positive, we have the difference scheme as follows:
(3)
where pm,l,jn=p(x+mh,y+lh,z+jh,t+nτ)⁠; h is the spatial grid interval; t is the time; τ is the time step; m⁠, l⁠, j⁠, and n are the numbers of nodes on each direction; M is the number of FD coefficients; and am (⁠m=1,2,…,M⁠) are the staggered-grid FD coefficients.
The second-order approximation for the time derivative is given by
(4)
where Δt is the time-sampling step in the time direction.
Substituting approximations 3 and 4 into equation 1 results in
(5)
where τ=Δt denotes the time-sampling step.
Using the plane-wave theory, we obtain from equation 5 that (Liu and Sen, 2011)
(6)
where kx=k cos θ cos ϕ⁠, ky=k cos θ sin ϕ⁠, kz=k sin θ⁠, k=kx2+ky2+kz2⁠, r=vτ/h⁠, ω is the angular frequency, θ is the angle of wave propagation measured from the horizontal plane perpendicular to the z-axis, and ϕ is the azimuth of the plane wave.

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

Optimizing regularized staggered-grid finite-difference operators

Assuming that equation 6 is valid for all angles of wave propagation, we obtain
(7)
where kx⁠, ky⁠, and kz are as mentioned above.
Let a be the vector form of the FD coefficients, and denote the left side of equation 7 by
(8)
and the right side of equation 7 by
(9)
Our aim is to minimize the error between the temporal dispersion and the spatial dispersion for a fixed range of wavenumbers; i.e.,
(10)
where k is equally distributed between 0 and K¯⁠, and K¯ is determined by the wave velocity, seismic source, and spatial grid interval; i.e., K¯/Ktotal=fv/(2h) with f being the source frequency and v is the wave velocity.
We observe that direct minimization of the objective function Φ for the FD coefficient a may lead to unstable results. This may happen for an FD operator with long stencils (Liang et al., 2013a; Liang, 2013b; Zhang et al., 2011). Therefore, we resort to regularizing technique to restore stability. The regularization model is established as
(11)
where α>0 is a user-defined regularization parameter and D is a scale operator. In this paper, we choose D as an identity operator. Our new task is solving the minimization problem,
(12)
instead of equation 10.

Regularizing computation of the coefficients

Using the Taylor expansion, we obtain
(13)
where a+ is the initial value, Δa=a−a+ and
(14)
Taking the partial derivative of equation 11 with respect to the increment Δam and noting that am=Δam+a+⁠, and using the fact that the partial derivative should be zero when Jα is minimized, we have
(15)
where α>0 is as mentioned in equation 11,
(16)
l=1,2,…,m⁠, and
(17)
Inserting equations 14, 16, and 17 into equation 15, we obtain the linear algebraic equations for retrieving the FD coefficients in the following form:
(18)
where Aα∈RM×M⁠, x∈RM×1 denotes the FD coefficients, and b∈RM×1⁠, which are in the following forms, respectively:
(19)
in which Cl,: refers to the lth row of the matrix C and cl refers to the lth component of the vector c⁠, l=1,2,…,m⁠. The initial value of a+ can be random numbers. The Von Neumann stability condition for explicit FD schemes is r<1/κ(∑m=1M|am|)-1 (Liu and Sen, 2011), where r is defined as before and κ is the spatial dimension of the medium.
Note that C is symmetric, but C may be ill-conditioned for some FD operators. If we do not consider regularization, direct solution of Ca=c will be unstable. To overcome the ill-conditioned property, we apply the regularization technique as mentioned above. Certainly, we could also apply other regularization techniques to suppress oscillations of solutions (Wang et al., 2011). Thus, instead of minimization problem 12, an a priori constrained optimization problem could be applied as follows:
(20)
and
(21)
where Ω[·] is a function supplying some a priori information to the solution and η is a constant.
Currently, to solve for x⁠, the regularized solution can be obtained by solving
(22)
It is clear that when α=0⁠, the regularized solution reduces to the least-squares solution. In the following, we refer to the regularization method as the new method.

1D dispersion analysis

The dispersion δ is used to measure the dispersion effect of the 1D acoustic-wave equations, which is defined as follows (Liu and Sen, 2011):
(23)
where vFD denotes the numerical phase velocity, v is the true velocity, k is the wavenumber, and h is the spatial grid interval. The dispersion is measured by the distance of δ to 1. Clearly, when δ=1⁠, dispersion disappear. And in equation 23, kh∈[0,π]⁠. Figures 1 and 2 compare dispersion errors of the traditional FD operator (here, we mean the FD operator obtained using the Taylor expansion method in the spatial domain) and Liu and Sen’s method with our new method as M is chosen as 4, 6, and 8, respectively. In Figure 1a and 1b, we choose v=1500  m/s⁠, τ=1  ms⁠, and h=10  m⁠; whereas in Figure 2a and 2b, we choose v=4500  m/s⁠, τ=1  ms⁠, and h=10  m⁠. Comparison of the results indicates that as r increases, the traditional methods will be unstable; whereas for Liu and Sen’s method and our new method, the results are stable. Furthermore, our new method could preserve the dispersion relation within a larger interval of frequencies. As we mentioned before, when α=0⁠, our new method reduces to the least-squares method. In Figure 3, we present the comparison of our new method with the least-squares method. It is evident that the new method yields more stable results than the least-squares method, the reason is that the new method possesses the regularizing property when dealing with ill-posedness, because we tested that the condition number of the matrix C is about O(1010)–O(1017) for different values of M⁠.

2D dispersion analysis

The 2D dispersion δ(θ) is used to measure the dispersion effect, which is defined by Liu and Sen (2011) as
(24)
where q is defined by
(25)
The smaller the absolute value of δ(θ)⁠, the less dispersion of the FD operator will be.

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  m/s and v=4500  m/s⁠, and the same M⁠, h⁠, and τ⁠, 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 α=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

The 3D dispersion δ(θ,ϕ) is used to measure the dispersion effect, which is defined by Liu and Sen (2011) as
(26)
The smaller the absolute value of δ(θ,ϕ)⁠, the less dispersion of the FD operator will be.

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.

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

First, we consider a simple model with velocity 1500  m/s⁠, spatial sampling interval of 20 m, and M=6 for all FD operators. The seismic source is located in the center of the test area and the source function used is as follows:
(27)
where the frequency f0 is 39 Hz and t0 is chosen as 4/f0⁠. The second derivative of the source function is applied in simulating the wavefield. The spectrum of the seismic source is shown in Figure 7 and the dominant frequency is 13 Hz.

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  m/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  Hz⁠, 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  m⁠) in Figure 10a–10c. It is evident that the new method provides the best simulation results with least dispersion.

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.

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.