Next Article in Journal
Numerical Method Using Homotopic Iterative Functions Based on the via Point for the Joint-Space Trajectory Generation
Previous Article in Journal
The Effect of Dead-Time and Damping Ratio on the Relative Performance of MPC and PID on Second Order Systems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Efficient Finite-Difference Stencil with High-Order Temporal Accuracy for Scalar Wave Modeling

1
School of Information and Communication Engineering, University of Electronic Science and Technology of China, No.4, Section 2, North Jianshe Road, Chengdu 610054, China
2
Laboratory of Imaging Detection and Intelligent Perception, University of Electronic Science and Technology of China, Chengdu 610054, China
3
BGP Inc., CNPC, Zhuozhou 072751, China
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2023, 13(2), 1140; https://doi.org/10.3390/app13021140
Submission received: 12 December 2022 / Revised: 1 January 2023 / Accepted: 5 January 2023 / Published: 14 January 2023

Abstract

:
Solving a scalar wave equation by the finite-difference (FD) method is a key step for advanced seismic imaging, in which the numerical accuracy is significantly affected by the FD stencil. High-order spatial and temporal approximations of the FD stencil can effectively improve the numerical accuracy and mitigate dispersion error. However, the huge costs of high-order stenciling in computation and storage hinder the application of large-scale modeling. In this paper, we propose a new efficient FD stencil with high-order temporal accuracy for numerical seismic modeling. The new stencil has a radial shape, including a standard cross-stencil and a rotated cross-stencil with a ( π / 4 ) degree, and it can reach sixth-order accuracy in the time approximation. Compared with the well-known temporal high-order cross-rhombus stencil, the new stencil involves fewer grid nodes and thus has higher computational efficiency, especially in high-order cases. Dispersion and stability analyses show that the new stencil has great improvements in mitigating the dispersion error and stability problem compared with the conventional methods. Numerical accuracy and execution time analyses show that the new stencil is an economical and feasible method for large-scale modeling.

1. Introduction

Finite-difference methods are the most popular tool for understanding complex wave phenomena due to their straightforward implementation and relatively small computational cost [1,2,3,4,5,6,7]. Solving scalar wave equations through FD methods is a basic step in seismic depth migration [8] and velocity modeling [9,10,11]. An FD stencil defines the grid nodes related to the FD operator, such as the widely used cross-stencil involving a series of horizontal and vertical grid nodes. High-order FD stencils for approximating temporal derivatives have practical advantages in mitigating dispersion errors and stability problems.
To obtain an FD stencil with high-order time approximation, the authors of [1] presented a Lax–Wendroff scheme in which the high-order FD operators of the temporal derivatives are transformed into the FD operators of the mixed spatial derivatives. The authors of [12,13] further developed this method with fourth- and sixth-order temporal accuracies and derived a stability formula for the scalar wave equation. Afterward, in [14], the authors presented a new FD stencil with a rhombus shape. The rhombus stencil is also similar to the Lax–Wendroff method, and it can reach arbitrary even-order accuracy for both the temporal and spatial accuracies. To improve the computational efficiency, the authors of [15] presented a mixed stencil with a small rhombus stencil and a large cross stencil. They demonstrated that this cross-rhombus stencil can reach (2M)th-order spatial and (2N)th-order temporal accuracies with the time–space domain FD coefficients. The authors of [16,17] further developed the cross-rhombus stencil in a staggered grid and gave the analytical formula for solving FD coefficients with fourth- and sixth-order temporal accuracies. Then, the authors of [18] presented a method for solving FD coefficients with arbitrary even-order accuracies in a staggered grid and investigated the optimal FD coefficients using a combination of the Taylor series expansion and least squares method.
A standard cross-rhombus stencil with a high-order time approximation will involve more grid nodes, thus increasing the cost exponentially. For the problem of low computational efficiency in high-order cases, we propose a new temporal high-order FD stencil. The new stencil has a radial shape, including a standard cross-stencil and a rotated cross-stencil with a ( π / 4 ) degree, in which the rotated cross stencil is similar to the Lax–Wendroff approach to improve the temporal accuracy. Compared with the standard cross-rhombus stencil method, the new stencil involves fewer grid nodes, and this is significant for reducing computational costs in large-scale seismic modeling. We prove that the new stencil can reach sixth-order temporal accuracy with the time–space domain FD coefficients. Although the new stencil cannot guarantee higher accuracy than the sixth-order time approximation, it is a compromise scheme which sacrifices a little accuracy in exchange for a large increase in computational efficiency. Then, we use the Taylor-series expansion (TE) method [19] to obtain the FD coefficients of the proposed FD stencil. As an alternative, we present the dispersion relationship-preserving (DRP) method [6,7,20,21] to obtain the FD coefficients satisfying more elevated wavenumbers in the dispersion relationship. The DRP-based FD coefficients can effectively mitigate the dispersion errors from the high-wavenumber components. To verify the feasibility of the proposed method, we designed a series of experiments to analyze the performance of the new stencil. The results indicate that the new stencil is a significant improvement in mitigating dispersion error and computational efficiency compared with the conventional methods.

2. Method

2.1. Review of the Staggered-Grid FD Scheme with the Cross-Rhombus Stencil

The 2D constant density scalar wave equation for space ( x , z ) and time t can be expressed as follows:
1 v 2 2 p t 2 = 2 p x 2 + 2 p z 2 ,
where p ( x , z , t ) is the scalar wavefield and v is the velocity. The Taylor expansion of Equation (1) with respect to time [1] yields
p 1 2 p 0 + p 1 Δ t 2 v 2 2 p 0 x 2 + 2 p 0 z 2 + j = 0 N w j 2 N p 0 x 2 j z 2 N 2 j .
Here, p l = p ( x , z , t + l Δ t ) is the discretized scalar wavefield, Δ t is the time step, and w j represents the weights of the spatial mixed derivatives. Then, Equation (2) can be approximated by the FD operator with a cross-rhombus stencil such that
a 0 p 0 , 0 0 + m = 1 M a m p m , 0 0 + p m , 0 0 + p 0 , m 0 + p 0 , m 0 c r o s s s t e n c i l + m = 1 N 1 n = 1 N m b m , n p m , n 0 + p m , n 0 + p m , n 0 + p m , n 0 r h o m b u s s t e n c i l r 2 p 1 2 p 0 + p 1 ,
Here, p m , n l = p ( x + m h , z + n h , t + l Δ t ) , h represents the grid spacing and r = v Δ t / h is the Courant number, where the parameters M and N define the size of the cross- and rhombus stencils, respectively. The cross-rhombus stencil can achieve (2M)th-order spatial and (2N)th-order temporal accuracies when the FD coefficients are derived from the time–space domain dispersion relation [15].

2.2. A New Efficient FD Stencil with High-Order Spatial and Temporal Accuracies

To further improve the computational efficiency, we simplify the cross-rhombus stencil to a new FD stencil with a radial shape. As shown in Figure 1, the new stencil involves grid nodes only with ( π n / 4 ) degrees, where n is an arbitrary integer. Compared with the standard cross-rhombus stencil, the new stencil has less grid nodes, especially in the high-order cases. We can prove that the stencil still has a higher-order temporal approximation, and the details of the proof are shown in Appendix A.
According to the definition of the new stencil (Figure 1), we obtained the new FD scheme:
a 0 p 0 , 0 0 + m = 1 M a m p m , 0 0 + p m , 0 0 + p 0 , m 0 + p 0 , m 0 c r o s s s t e n c i l + n = 1 N 1 b n p n , n 0 + p n , n 0 + p n , n 0 + p n , n 0 ) r o t a t e d c r o s s s t e n c i l w i t h a ( π / 4 ) d e g r e e r 2 p 1 2 p 0 + p 1 .
The new radiation stencil is a combination of a standard cross-stencil and a rotated cross-stencil with ( π / 4 ) degrees (Figure 1). Let a m represent the FD coefficients of the standard cross-stencil and b n represent the FD coefficients of the rotated cross-stencil. In the following parts, we present two methods to obtain the FD coefficients of the proposed stencil and analyze the approximate accuracy theoretically.

2.3. Determining FD Coefficients of the New Stencil through Taylor-Series Expansion

Assuming the plane wave propagates in the grid, we let
p i , j l = e i k x ( x + i h ) + k z ( z + j h ) ω ( t + l Δ t ) .
Here, k x = k cos ( θ ) , k z = sin ( θ ) , and θ is the propagation direction of the plane wave, while ω represents the angular frequency. By substituting Equation (5) into Equation (4), we obtain
a 0 + 2 m = 1 M a m cos ( m k x h ) + cos ( m k z h ) + 4 n = 1 N 1 b n cos ( n k x h ) cos ( n k z h ) r 2 2 + 2 cos ( ω Δ t ) .
Equation (6) represents the time–space domain dispersion relation of our new FD scheme. Let v = ω / k . The Taylor-series expansions of the cosine functions cos ( m k x h ) , cos ( m k z h ) and cos ( ω Δ t ) give
a 0 + 2 m = 1 M a m j = 0 ( 1 ) j m 2 j ( k x 2 j + k z 2 j ) h 2 j 2 j ! + 4 n = 1 N 1 b n ξ = 0 ζ = 0 1 ξ + ζ n 2 ( ξ + ζ ) 2 ξ ! 2 ζ ! k x 2 ξ k z 2 ζ h 2 ξ + ζ 2 j = 1 1 j r 2 j 2 k h 2 j 2 j ! .
When comparing the coefficients of h 0 ( j = 0 ), we obtain
a 0 + 4 m = 1 M a m + 4 n = 1 N 1 b n = 0 .
Let
k 2 j = k x 2 + k z 2 j = ξ = 0 j j ! ξ ! j ξ ! k x 2 ξ k z 2 j 2 ξ .
By substituting Equation (9) into Equation (7) and comparing the weights of h 2 j ( j > 0 ), we obtain
m = 1 M a m ( 1 ) j m 2 j ( k x 2 j + k z 2 j ) 2 j ! + 2 n = 1 N 1 b n ξ = 0 j 1 j n 2 j 2 ξ ! 2 j 2 ξ ! k x 2 ξ k z 2 j 2 ξ ξ = 0 j ( 1 ) j j ! r 2 j 2 2 j ! ξ ! ( j ξ ) ! k x 2 ξ k z 2 j 2 ξ .
When comparing the weights of k x 2 ξ k z 2 j 2 ξ , we obtain
m = 1 M m 2 j a m + 2 n = 1 N 1 n 2 j b n = r 2 j 2 ( ξ = 0 o r ξ = j ) 2 n = 1 N 1 n 2 j 2 j 2 ξ ! 2 ξ ! b n = j ! r 2 j 2 2 j ! ξ ! j ξ ! ( ξ = 1 , 2 , , j 1 ) .
We only need M + N 1 independent equations to solve the FD coefficients a m and b n . Through redundancy analysis, we simplify Equation (11) as
m = 1 M m 2 j a m + 2 n = 1 N 1 n 2 j b n = r 2 j 2 ( j = 1 , 2 , , M ) ; n = 1 N 1 n 2 j b n = j ! 2 ξ ! 2 j 2 ξ ! 2 2 j ! ξ ! j ξ ! r 2 j 2 ( j = 2 , 3 , , N ) ; a n d ξ = i n t j / 2 ) .
where i n t represents a function to find the integer part of the real value. In Appendix A, we prove that the FD coefficients obtained by Equation (12) can reach at least sixth-order temporal accuracy. Aside from that, we give the solution of Equation (12) in an analytic form when N = 2 , and the details are shown in Appendix B. Next, we analyze the relationship between the new stencil and some existing methods:
(1) When M > 0 and N > 0 (for example, M = 3 and N = 3 ), this is a general situation, and Equation (12) can be expressed as
1 4 9 2 8 18 1 16 81 2 32 162 1 16 81 2 128 1458 0 0 0 1 16 81 0 0 0 1 64 729 0 0 0 1 256 6561 a 1 a 2 a 3 b 1 b 2 b 3 = 1 r 2 r 4 r 2 / 6 r 4 / 10 3 r 6 / 10 .
The FD coefficients a m and b n can be determined by solving the above linear equation, and a 0 can be determined by Equation (8).
(2) When N = 0 , there is no rotated cross-stencil. Equation (12) can be simplified as follows:
m = 1 M m 2 j a m = r 2 j 2 ( j = 1 , 2 , , M ) .
The new FD scheme is simplified to the conventional time–space domain FD scheme with a cross-stencil [19].
(3) When N = 1 , Equation (12) becomes the fourth-order temporal accuracy FD scheme with a cross-rhombus stencil [14,15].
(4) When r = 0 , the FD coefficients b n = 0 , Equation (12) can be expressed as
m = 1 M m 2 a m = 1 ( j = 1 ) m = 1 M m 2 j a m = 0 ( j = 2 , , M ) .
In this case, the new FD scheme is simplified to the space domain FD scheme with a cross-stencil.

2.4. Determining the FD Coefficients of the New Stencil with the Dispersion Relationship-Preserving Method

TE-based FD coefficients are prone to dispersion errors for the high-wavenumber components [17,22]. As an alternative, we present the dispersion relationship-preserving method to obtain FD coefficients satisfying more elevated wavenumbers in the dispersion relationship [6,7,20,21]. The DRP-based FD coefficients can effectively mitigate the dispersion errors from the high-wavenumber components. Following our previous work [7], we define a new function ψ m , β , θ to represent the weights of the FD coefficients a m in the dispersion relation (Equation (6)). Then, ψ m , β , θ can be defined as
ψ m , β , θ = 2 cos ( m β cos ( θ ) ) + cos ( m β sin ( θ ) ) .
Here, β = k h , and θ represents the propagation angle. Similarly, we define another function φ n , β , θ to represent the weights of b n . The function φ n , β , θ can be defined as
φ n , β , θ = 4 cos ( n β cos ( θ ) ) cos ( n β sin ( θ ) ) .
Then, we can extend ψ m , β , θ to matrix A ( θ ) involving a series of β and a fixed angle θ . Matrix A ( θ ) is
A ( θ ) = 1 ψ 1 , β 1 , θ ψ 2 , β 1 , θ ψ M , β 1 , θ 1 ψ 1 , β 2 , θ ψ 2 , β 2 , θ ψ M , β 2 , θ 1 ψ 1 , β ξ , θ ψ 2 , β ξ , θ ψ M , β ξ , θ ,
where β i = β m a x / ξ i and β m a x can be obtained by the optimization algorithm [23]. We also extend the function φ n , β , θ to the matrix
B ( θ ) = φ 1 , β 1 , θ φ 2 , β 1 , θ φ N 1 , β 1 , θ φ 1 , β 2 , θ φ 2 , β 2 , θ φ N 1 , β 2 , θ φ 1 , β ξ , θ φ 2 , β ξ , θ φ N 1 , β ξ , θ .
Additionally, the right-hand side of Equation (6) is expanded to
D ( θ ) = r 2 2 + 2 cos ( β 1 r ) r 2 2 + 2 cos ( β 2 r ) r 2 2 + 2 cos ( β ξ r ) .
Thus, the time–space domain dispersion relation of our new method involving ξ wavenumbers and ζ angles can be expressed as follows:
A ( θ 1 ) B ( θ 1 ) A ( θ 2 ) B ( θ 2 ) A ( θ ζ ) B ( θ ζ ) a 0 a 1 b N 1 = D ( θ 1 ) D ( θ 2 ) D ( θ ζ ) .
Then, numerically solving this over-determined system can yield the DRP-based FD coefficients.

3. Dispersion Analysis

We denote the phase velocity v F D = ω / | k | and then yield the phase velocity ratio
δ ( β , θ ) = v FD v = 1 r β arccos ( 1 + 0.5 r 2 q ( β , θ ) ) .
Here, q is
q = a 0 , 0 + 2 m = 1 M a m cos ( m β cos ( θ ) ) + cos ( m β sin ( θ ) ) + 4 n = 1 N 1 b n cos ( n β cos ( θ ) ) cos ( n β sin ( θ ) ) ,
In addition, β = k h . We analyze and compare the ratios δ of different FD methods, and the abbreviations of these FD methods are listed in Table 1. Figure 2 shows the ratios δ varying with parameters β and θ .
It is clear that the ratios δ of the TE-C method were far larger than one for the high-wavenumber region (Figure 2a), resulting in a serious temporal dispersion error. The corresponding ratios of the new method (Figure 2b) were significantly improved, indicating that the TE-R method can effectively mitigate the temporal dispersion error for the high-wavenumber components. Figure 2c shows that the DRP-R method had the smallest dispersion error compared with the other methods.
The TE-CR method will multiply the computational cost by the increase in the temporal order N. Then, we analyzed the ratios δ of the TE-CR and TE-R methods with different orders N, and the results are shown in Figure 3. It can be seen that the performance of the TE-R method was slightly less than that of the TE-CR method, but the cost of the TE-CR method increased exponentially in the high-order cases. Although the new stencil sacrificed a little approximation accuracy, it had great improvement in computational efficiency.

4. Stability Analysis

The function cos ( ω Δ t ) in Equation (6) must lie in the interval [ 1 , 1 ] , where
1 cos ( ω Δ t ) = 1 + 0.5 r 2 q 1 .
We mostly use the Nyquist wavenumber to analyze the stability [19,24] such that
k x h = k z h = π .
By substituting Equation (25) into Equation (24), we obtain the 2D stability condition of the proposed FD scheme as follows:
r m = 1 M a m ( 1 ) m 1 + 1 + n = 1 N b n ( 1 ) 2 n 1 + 1 1 / 2 .
Then, the 2D stability factor s can be defined as
s = m = 1 M a m ( 1 ) m 1 + 1 + n = 1 N 1 b n ( 1 ) 2 n 1 + 1 1 / 2 ,
where the stability factor s is related to the FD coefficients, and these FD coefficients are determined by the Courant number r. Therefore, the FD scheme must satisfy the following stability condition:
r s .
In the following, we analyze the stability factors s varying with the Courant number r, and the results are shown in Figure 4. It can be seen that the TE-C method had the smallest stability factor, indicating that it is prone to instability. The stability factors of the TE-R and DRP-R methods were slightly less than the TE-CR method but much larger than that of the TE-C method. Figure 4b shows the maximum factor s satisfying r s for different orders M. The maximum stability factors of the TE-R and DRP-R methods were obviously larger than that of the TE-C method, which showed an improvement in stability.

5. Numerical Experiments

5.1. Seismic Modeling of a 2D Homogeneous Velocity Model

In this section, we use a 2D homogeneous velocity model to examine our new stencil. The 2D homogeneous model had 512 × 512 grid nodes with a velocity υ = 1500 (m/s). We set M = 8 , N = 8 and a time step Δ t = 0.0015 (s) for the FD methods. A Ricker wavelet with 40 (Hz) was injected as an explosion source. Figure 5 shows the snapshots of different methods with grid sampling h = 4 (m) and h = 6 (m).
It can be seen that the TE-C method had visible temporal dispersion errors (red arrows) with the grid sampling h = 4 (m) and had both serious spatial (white arrows) and temporal dispersion errors with the grid sampling h = 6 (m). The corresponding temporal dispersion errors were reduced in the TE-R method, indicating that the new stencil had better temporal accuracy. The spatial dispersion errors were significantly reduced in the DRP-R method. Thus, the DRP-based FD coefficients could effectively mitigate the spatial dispersion with a large grid spacing.

5.2. Numerical Accuracy and CPU Execution Time Analyses

An increase in the temporal order N will significantly increase the computational cost. In this part, we set a series of temporal orders N to analyze the relative errors and CPU execution times (CPU: AMD 5600h) of the cross-rhombus stencil and our new stencil. Here, we adopted h = 6 (m), Δ t = 0.0015 (s) and υ = 1500 (m/s) for all FD methods. We set up three receivers at positions ( 768 , 1536 ) (m), ( 1536 , 768 ) (m) and ( 768 , 768 ) (m), which represent the simulations of different propagation angles. We calculated their total relative errors with the reference method [7], and the total relative error was defined as
E = i = 1 3 t i r i 2 ,
where t i represents the records of the three receivers and r i represents the records of the references.
Table 2 lists the total relative errors and CPU execution times of the TE-R and TE-CR methods with different orders N. It can be seen that the relative errors of the TE-CR method were slightly less than that of the TE-R method (Table 2). However, the CPU execution times of the TE-CR method were much larger than that of the TE-R method. Table 2 demonstrates that the numerical accuracy of the new method was slightly lower than that of the standard cross-rhombus stencil, but its computational efficiency was greatly improved, especially in the high-order cases.

5.3. Seismic Modeling of the 2D Inhomogeneous Velocity Model

We used a modified 2D Sigsbee velocity model to test the new stencil in the complex velocity model. The modified Sigsbee model (Figure 6) has 500 × 1200 grid nodes, with the velocities varying from 1550 to 4500 (m/s). We set M = 8 , N = 8 , grid spacing h = 8.5 (m) and the time step Δ t = 0.0009 (s) for all FD schemes. The source with the peak frequency f p = 30 (Hz) (Ricker wavelet) was injected at the center of the surface. For convenience, we ignored the influence of the discrete grid representation on the strong material interfaces [25,26].
Figure 7 shows the snapshots and seismic records of different FD methods. The TE-C method had obvious dispersion errors in the low-velocity layers, while the corresponding dispersion errors were reduced in the TE-R and DRP-R methods. Figure 8 shows the seismic records at position ( 2125 , 0 ) (m), at which the reference trace r used a high-order temporal and spatial FD scheme ( M = 12 , N = 12 ) with the DRP-based coefficients [7]. Traces a-r b-r and c-r represent the relative errors of the TE-C, TE-R and DRP-R methods compared with the reference method, respectively. The dispersion errors of the TE-C method were serious for the first arrival waves and reflection waves. The corresponding errors were significantly alleviated in the TE-R and DRP-R methods. Table 3 lists the relative errors of all FD methods. It is clear that the relative errors of the TE-R and DRP-R methods were far less than that of the TE-C method. The numerical experiments with the Sigsbee model demonstrated that the new stencil also had better performance in the complex velocity model for mitigating dispersion errors.

6. Conclusions

In this work, we proposed a new FD stencil with high-order temporal and spatial accuracies. The new stencil can reach at least sixth-order temporal accuracy and involves fewer grid nodes than that of the well-known cross-rhombus stencil. We presented two methods to solve the FD coefficients of the new stencil by Taylor-series expansion and the dispersion relationship-preserving method, in which the DRP-based FD coefficients can effectively mitigate the dispersion errors from the high-wavenumber components. Dispersion, stability analyses and numerical experiments demonstrated that the new stencil had better accuracy and stability than those of the conventional cross-stencil. Numerical accuracy and execution time analyses proved that the accuracy of the new stencil was slightly lower than that of the standard cross-rhombus stencil in the high-order cases, but its computational efficiency was greatly improved. Thus, the new stencil is a more economical and feasible method for large-scale modeling.

Author Contributions

Conceptualization, G.C. and Y.L.; methodology, G.C.; validation, Z.P.; formal analysis, G.C.; investigation, G.C.; resources, Y.L.; writing—original draft preparation, G.C.; writing—review and editing, Z.P. and Y.L.; visualization, G.C.; supervision, Z.P.; project administration, Z.P.; funding acquisition, Z.P. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Natural Science Foundation of the Sichuan Province of China (Grant No. 2022NSFSC40574) and partially supported by the National Natural Science Foundation of China (Grant No. 61775030 and Grant No. 61571096).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This study benefited from the reproducible codes and documentation provided in Madagascar open-source software [27]. The authors appreciate the owner of the copyrights.

Conflicts of Interest

The authors declare no conflict of interest.The funders had no role in the design of the study or in the collection or analyses.

Appendix A

In this part, we analyze the accuracy of our new stencil with the time–space domain TE-based FD coefficients (Equation (12)). First, we define the absolute error ε according to Equation (6) as follows:
ε = a 0 + 2 m = 1 M a m cos ( m k x h ) + cos ( m k z h ) + 4 n = 1 N 1 b n cos ( n k x h ) cos ( n k z h ) r 2 2 + 2 cos ( ω Δ t ) .
By substituting Equations (7) and (10) into Equation (A1), we obtain
ε = a 0 + 2 m = 1 M a m j = 0 ( 1 ) j m 2 j ( k x 2 j + k z 2 j ) h 2 j 2 j ! + 4 n = 1 N 1 b n ξ = 0 ζ = 0 1 ξ + ζ n 2 ( ξ + ζ ) 2 ξ ! 2 ζ ! k x 2 ξ k z 2 ζ h 2 ξ + ζ 2 r 2 j = 1 1 j r 2 j k h 2 j 2 j ! = a 0 + 2 m = 1 M a m j = 0 ( 1 ) j m 2 j ( k x 2 j + k z 2 j ) h 2 j 2 j ! + 4 n = 1 N 1 b n j = 1 ( 1 ) j n 2 j ( k x 2 j + k z 2 j ) h 2 j 2 j ! + 4 n = 1 N 1 b n ξ = 1 ζ = 1 1 ξ + ζ n 2 ( ξ + ζ ) 2 ξ ! 2 ζ ! k x 2 ξ k z 2 ζ h 2 ξ + ζ 2 r 2 j = 1 1 j r 2 j k h 2 j 2 j ! = a 0 + 4 m = 1 M a m + 4 n = 1 N 1 b n + 2 j = 1 ( 1 ) j m = 1 M a m m 2 j ( k x 2 j + k z 2 j ) 2 j ! + 2 n = 1 N 1 b n n 2 j ( k x 2 j + k z 2 j ) 2 j ! h 2 j + 4 n = 1 N 1 b n ξ = 1 ζ = 1 1 ξ + ζ n 2 ( ξ + ζ ) 2 ξ ! 2 ζ ! k x 2 ξ k z 2 ζ h 2 ξ + ζ 2 r 2 j = 1 1 j r 2 j ( k x 2 j + k z 2 j ) h 2 j 2 j ! 2 r 2 j = 1 ξ = 1 j 1 ( 1 ) j r 2 j k x 2 ξ k z 2 j 2 ξ j ! h 2 j ( 2 j ) ! ( ξ ) ! ( j ξ ) ! 2 j = M + 1 ( 1 ) j m = 1 M a m m 2 j ( k x 2 j + k z 2 j ) 2 j ! + 2 n = 1 N 1 b n n 2 j ( k x 2 j + k z 2 j ) 2 j ! h 2 j + 4 n = 1 N 1 b n ξ = 1 , ζ = 1 , ξ + ζ > 3 , 1 ξ + ζ n 2 ( ξ + ζ ) 2 ξ ! 2 ζ ! k x 2 ξ k z 2 ζ h 2 ξ + ζ 2 r 2 j = M + 1 1 j v 2 j Δ t 2 j ( k x 2 j + k z 2 j ) 2 j ! 2 r 2 j = 4 ξ = 1 , ξ i n t ( j / 2 ) j 1 ( 1 ) j v 2 j Δ t 2 j k x 2 ξ k z 2 j 2 ξ j ! ( 2 j ) ! ( ξ ) ! ( j ξ ) ! o ( h 2 M ) + o ( h 6 ) o ( Δ t 2 M ) o ( Δ t 6 ) o ( h 6 ) + o ( Δ t 6 ) .
where the minimum power of h and Δ t is six, and thus the modeling accuracy is of the sixth order.

Appendix B

When N = 2 , Equation (12) gives
b 1 = r 2 6 .
Then, we rewrite Equation (12) as the following Vandermonde matrix form:
1 0 2 0 M 0 1 2 2 2 M 2 1 2 M 2 2 2 M 2 M 2 M 2 1 2 a 1 + r 2 / 3 2 2 a 2 M 2 a m = 1 r 2 r 2 M 2 .
By solving this linear system, we obtain the analytic FD coefficients:
b 1 = r 2 6 ; a m = ( 1 ) m + 1 1 n M , n m | n 2 r 2 n 2 m 2 | r 2 3 ( m = 1 ) ; a m = ( 1 ) m + 1 m 2 1 n M , n m | n 2 r 2 n 2 m 2 | ( m = 2 , , M ) .

References

  1. Dablain, M. The application of high-order differencing to the scalar wave equation. Geophysics 1986, 51, 54–66. [Google Scholar] [CrossRef]
  2. van Vossen, R.; Robertsson, J.O.; Chapman, C.H. Finite-difference modeling of wave propagation in a fluid-solid configuration. Geophysics 2002, 67, 618–624. [Google Scholar] [CrossRef]
  3. Moczo, P.; Kristek, J.; Halada, L. 3D fourth-order staggered-grid finite-difference schemes: Stability and grid dispersion. Bull. Seismol. Soc. Am. 2000, 90, 587–603. [Google Scholar] [CrossRef]
  4. Etgen, J.T.; O’ Brien, M.J. Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial. Geophysics 2007, 72, SM223–SM230. [Google Scholar] [CrossRef] [Green Version]
  5. Moczo, P.; Kristek, J.; Gális, M. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures; Cambridge University Press: Cambridge, UK, 2014. [Google Scholar]
  6. Etemadsaeed, L.; Moczo, P.; Kristek, J.; Ansari, A.; Kristekova, M. A no-cost improved velocity-stress staggered-grid finite-difference scheme for modelling seismic wave propagation. Geophys. J. Int. 2016, 207, 481–511. [Google Scholar] [CrossRef] [Green Version]
  7. Chen, G.; Wang, Y.; Wang, Z.; Zhang, S. Dispersion-relationship-preserving seismic modelling using the cross-rhombus stencil with the finite-difference coefficients solved by an over-determined linear system. Geophys. Prospect. 2020, 68, 1771–1792. [Google Scholar] [CrossRef]
  8. Baysal, E.; Dan, D.K.; Sherwood, J. Reverse-Time Migration. Geophysics 1983, 48, 1514–1524. [Google Scholar] [CrossRef] [Green Version]
  9. Virieux, J.; Operto, S. An overview of full-waveform inversion in exploration geophysics. Geophysics 2009, 74, WCC1–WCC26. [Google Scholar] [CrossRef]
  10. Chen, G.; Wang, Z. Robust full-waveform inversion based on particle swarm optimization. In SEG Technical Program Expanded Abstracts 2017; Society of Exploration Geophysicists: Houston, TX, USA, 2017; pp. 1302–1306. [Google Scholar] [CrossRef]
  11. Zhang, C.; Fan, L.; Chen, G.; Li, J. AVO-Friendly Velocity Analysis Based on the High-Resolution PCA-Weighted Semblance. Appl. Sci. 2022, 12, 6098. [Google Scholar] [CrossRef]
  12. Chen, J.B. High-order time discretizations in seismic modeling. Geophysics 2007, 72, SM115–SM122. [Google Scholar] [CrossRef]
  13. Chen, J.B. A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation. Geophysics 2011, 76, T37–T42. [Google Scholar] [CrossRef] [Green Version]
  14. Liu, Y.; Sen, M.K. Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation. J. Comput. Phys. 2013, 232, 327–345. [Google Scholar] [CrossRef]
  15. Wang, E.; Liu, Y.; Sen, M.K. Effective finite-difference modelling methods with 2-D acoustic wave equation using a combination of cross and rhombus stencils. Geophys. J. Int. 2016, 206, 1933–1958. [Google Scholar] [CrossRef]
  16. Tan, S.; Huang, L. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation. Geophys. J. Int. 2014, 197, 1250–1267. [Google Scholar] [CrossRef]
  17. Tan, S.; Huang, L. A staggered-grid finite-difference scheme optimized in the time–space domain for modeling scalar-wave propagation in geophysical problems. J. Comput. Phys. 2014, 276, 613–634. [Google Scholar] [CrossRef] [Green Version]
  18. Ren, Z.; Li, Z.C. Temporal high-order staggered-grid finite-difference schemes for elastic wave propagationTemporal high-order SFD schemes. Geophysics 2017, 82, T207–T224. [Google Scholar] [CrossRef]
  19. Liu, Y.; Sen, M.K. A new time–Space domain high-order finite-difference method for the acoustic wave equation. J. Comput. Phys. 2009, 228, 8779–8806. [Google Scholar] [CrossRef]
  20. Ye, F.; Chu, C. Dispersion-relation-preserving finite difference operators: Derivation and application. In SEG Technical Program Expanded Abstracts 2005; Society of Exploration Geophysicists: Houston, TX, USA, 2005; pp. 1783–1786. [Google Scholar]
  21. Liang, W.; Wang, Y.; Yang, C. Determining finite difference weights for the acoustic wave equation by a new dispersion-relationship-preserving method. Geophys. Prospect. 2015, 63, 11–22. [Google Scholar] [CrossRef]
  22. Liu, Y. Globally optimal finite-difference schemes based on least squares. Geophysics 2013, 78, T113–T132. [Google Scholar] [CrossRef]
  23. Chen, G.; Peng, Z.; Li, Y. A framework for automatically choosing the optimal parameters of finite-difference scheme in the acoustic wave modeling. Comput. Geosci. 2022, 159, 104948. [Google Scholar] [CrossRef]
  24. Liu, Y.; Sen, M.K. Scalar Wave Equation Modeling with Time-Space Domain Dispersion-Relation-Based Staggered-Grid Finite-Difference Schemes. Bull. Seismol. Soc. Am. 2011, 101, 141–159. [Google Scholar] [CrossRef]
  25. Moczo, P.; Gregor, D.; Kristek, J.; de la Puente, J. A discrete representation of material heterogeneity for the finite-difference modelling of seismic wave propagation in a poroelastic medium. Geophys. J. Int. 2019, 216, 1072–1099. [Google Scholar] [CrossRef] [Green Version]
  26. Gregor, D.; Moczo, P.; Kristek, J.; Mesgouez, A.; Lefeuve-Mesgouez, G.; Kristekova, M. Subcell-resolution finite-difference modelling of seismic waves in Biot and JKD poroelastic media. Geophys. J. Int. 2020, 224, 760–794. [Google Scholar] [CrossRef]
  27. Fomel, S.; Sava, P.; Vlad, I.; Liu, Y.; Bashkardin, V. Madagascar: Open-source software project for multidimensional data analysis and reproducible computational experiments. J. Open Res. Softw. 2013, 1, e8. [Google Scholar]
Figure 1. Schematic diagram of the new stencil, where the new radiation stencil is composed of a standard cross-stencil and a rotated cross-stencil of ( π / 4 ) degrees.
Figure 1. Schematic diagram of the new stencil, where the new radiation stencil is composed of a standard cross-stencil and a rotated cross-stencil of ( π / 4 ) degrees.
Applsci 13 01140 g001
Figure 2. Dispersion curves varying with the parameters β = k h and propagation angle θ . Here, Δ t = 0.0015 (s), υ = 1500 (m/s), h = 6 (m), M = 8 , and N = 8 : (a) TE-C method, (b) TE-R method and (c) DRP-R method.
Figure 2. Dispersion curves varying with the parameters β = k h and propagation angle θ . Here, Δ t = 0.0015 (s), υ = 1500 (m/s), h = 6 (m), M = 8 , and N = 8 : (a) TE-C method, (b) TE-R method and (c) DRP-R method.
Applsci 13 01140 g002
Figure 3. Dispersion curves of TE-CR and TE-R methods with different orders N. Here, υ = 1500 (m/s), h = 6 (m), M = 8 , and propagation angle α = π / 16 : (a) TE-CR method and (b) TE-R method.
Figure 3. Dispersion curves of TE-CR and TE-R methods with different orders N. Here, υ = 1500 (m/s), h = 6 (m), M = 8 , and propagation angle α = π / 16 : (a) TE-CR method and (b) TE-R method.
Applsci 13 01140 g003
Figure 4. Stability curves of different FD methods: (a) stability factors varying with the Courant numbers r and (b) maximum stability factors satisfying r s with different orders M.
Figure 4. Stability curves of different FD methods: (a) stability factors varying with the Courant numbers r and (b) maximum stability factors satisfying r s with different orders M.
Applsci 13 01140 g004
Figure 5. Numerical experiments of the homogeneous model with grid sampling h = 4 (m) and h = 6 (m). The left panel shows the snapshots with the grid sampling h = 4 (m), and the right panel represents the grid sampling h = 6 (m): (a,b) TE-C method, (c,d) TE-R method and (e,f) DRP-R method.
Figure 5. Numerical experiments of the homogeneous model with grid sampling h = 4 (m) and h = 6 (m). The left panel shows the snapshots with the grid sampling h = 4 (m), and the right panel represents the grid sampling h = 6 (m): (a,b) TE-C method, (c,d) TE-R method and (e,f) DRP-R method.
Applsci 13 01140 g005aApplsci 13 01140 g005b
Figure 6. The 2D Sigsbee model had 500 × 1200 grid nodes with a grid spacing h = 8.5 (m) and variation in the velocities from 1550 to 4500 (m/s).
Figure 6. The 2D Sigsbee model had 500 × 1200 grid nodes with a grid spacing h = 8.5 (m) and variation in the velocities from 1550 to 4500 (m/s).
Applsci 13 01140 g006
Figure 7. Snapshots and corresponding seismic records of the 2D Sigsbee velocity model. Here, M = 8 , N = M , Δ t = 0.0009 (s) and h = 8.5 (m). (a,b) The TE-C method. (c,d) The TE-R method. (e,f) The DRP-R method.
Figure 7. Snapshots and corresponding seismic records of the 2D Sigsbee velocity model. Here, M = 8 , N = M , Δ t = 0.0009 (s) and h = 8.5 (m). (a,b) The TE-C method. (c,d) The TE-R method. (e,f) The DRP-R method.
Applsci 13 01140 g007
Figure 8. Seismogram at position ( 0 , 2125 ) (m) of different FD methods. Trace a represents the record of the TE-C method. Trace b represents the TE-R method. Trace c represents the DRP-R method. Trace r represents the reference trace. Trace a-r represents the relative error between trace a and reference trace r. Trace b-r represents the relative error of trace b. Trace c-r represents the relative error of trace c.
Figure 8. Seismogram at position ( 0 , 2125 ) (m) of different FD methods. Trace a represents the record of the TE-C method. Trace b represents the TE-R method. Trace c represents the DRP-R method. Trace r represents the reference trace. Trace a-r represents the relative error between trace a and reference trace r. Trace b-r represents the relative error of trace b. Trace c-r represents the relative error of trace c.
Applsci 13 01140 g008
Table 1. Abbreviation table of different FD methods used for dispersion analyses.
Table 1. Abbreviation table of different FD methods used for dispersion analyses.
AbbreviationsFD CoefficientsFD Stencils
TE-CTE-based FD coefficientsCross-stencil
TE-CRTE-based FD coefficientsCross-rhombus stencil
TE-RTE-based FD coefficientsRadiation stencil
DRP-RDRP-based FD coefficientsRadiation stencil
Table 2. Relative errors and CPU execution times of the TE-CR and TE-R methods with different orders N.
Table 2. Relative errors and CPU execution times of the TE-CR and TE-R methods with different orders N.
CasesMethodsMNTotal Relative Errors ( 10 1 Pa)Execution Times (s)
1TE-C1242.6383241.22905
2TE-CR1262.0253480.51018
3TE-R1263.8692653.72843
4TE-CR1281.48282109.32694
5TE-R1282.5116561.30650
6TE-CR12100.85564156.39934
7TE-R12101.5266365.83632
8TE-CR12120.24282204.95636
9TE-R12120.7244670.14191
Table 3. Relative errors of different FD methods in the 2D Sigsbee velocity model.
Table 3. Relative errors of different FD methods in the 2D Sigsbee velocity model.
OrdersTotal Relative Errors (Pa)
MNTE-CTE-CRTE-RDRP-R
887.9747632.2032472.3482450.715832
842.3903452.4998540.846557
6612.5221964.3802394.4737681.091585
634.5023924.6285711.154466
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Chen, G.; Peng, Z.; Li, Y. An Efficient Finite-Difference Stencil with High-Order Temporal Accuracy for Scalar Wave Modeling. Appl. Sci. 2023, 13, 1140. https://doi.org/10.3390/app13021140

AMA Style

Chen G, Peng Z, Li Y. An Efficient Finite-Difference Stencil with High-Order Temporal Accuracy for Scalar Wave Modeling. Applied Sciences. 2023; 13(2):1140. https://doi.org/10.3390/app13021140

Chicago/Turabian Style

Chen, Guiting, Zhenming Peng, and Yalin Li. 2023. "An Efficient Finite-Difference Stencil with High-Order Temporal Accuracy for Scalar Wave Modeling" Applied Sciences 13, no. 2: 1140. https://doi.org/10.3390/app13021140

APA Style

Chen, G., Peng, Z., & Li, Y. (2023). An Efficient Finite-Difference Stencil with High-Order Temporal Accuracy for Scalar Wave Modeling. Applied Sciences, 13(2), 1140. https://doi.org/10.3390/app13021140

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop