Next Article in Journal
Two-View Structure-from-Motion with Multiple Feature Detector Operators
Next Article in Special Issue
Strong-Scattering Multiparameter Reconstruction Based on Elastic Direct Envelope Inversion and Full-Waveform Inversion with Anisotropic Total Variation Constraint
Previous Article in Journal
Observations of Ionospheric Clutter at Near Equatorial High Frequency Radar Stations
Previous Article in Special Issue
Automatic Horizon Picking Using Multiple Seismic Attributes and Markov Decision Process
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation

School of Earth Science, China University of Petroleum (East China), Qingdao 266000, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(3), 604; https://doi.org/10.3390/rs15030604
Submission received: 30 November 2022 / Revised: 10 January 2023 / Accepted: 10 January 2023 / Published: 19 January 2023
(This article belongs to the Special Issue Geophysical Data Processing in Remote Sensing Imagery)

Abstract

:
High-precision finite difference (FD) wavefield simulation is one of the key steps for the successful implementation of full-waveform inversion and reverse time migration. Most explicit FD schemes for solving seismic wave equations are not compact, which leads to difficulty and low efficiency in boundary condition treatment. Firstly, we review a family of tridiagonal compact FD (CFD) schemes of various orders and derive the corresponding optimization schemes by minimizing the error between the true and numerical wavenumber. Then, the optimized CFD (OCFD) schemes and a second-order central FD scheme are used to approximate the spatial and temporal derivatives of the 2D acoustic wave equation, respectively. The accuracy curves display that the CFD schemes are superior to the central FD schemes of the same order, and the OCFD schemes outperform the CFD schemes in certain wavenumber ranges. The dispersion analysis and a homogeneous model test indicate that increasing the upper limit of the integral function helps to reduce the spatial error but is not conducive to ensuring temporal accuracy. Furthermore, we examine the accuracy of the OCFD schemes in the wavefield modeling of complex structures using a Marmousi model. The results demonstrate that the OCFD4 schemes are capable of providing a more accurate wavefield than the CFD4 scheme when the upper limit of the integral function is 0.5π and 0.75π.

Graphical Abstract

1. Introduction

Seismic forward modeling plays an important role in seismic data processing and interpretation and is the basis of migration and full-waveform inversion [1,2,3,4,5]. Acoustic wave equation has attracted wide interest in 2D and 3D seismic exploration due to its advantages, such as low computational costs and storage requirements. To date, finite difference (FD) methods [6,7,8,9], finite-element methods [10,11], spectral-element methods [12,13,14,15], boundary integral methods [16,17,18,19], pseudo-spectral methods [20,21,22], and finite volume methods [23] have been used to calculate the propagation of seismic waves in various media. Among them, FD methods are widely used in exploration seismology for their simple implementation and high efficiency.
Conventional FD methods can be simply classified into two categories: explicit and implicit (compact) schemes. Explicit FD schemes are based on recursion relations and are popular for numerically approximating wave equations due to their lower computational costs and simplicity [24,25]. However, the computational accuracy of lower-order explicit FD methods is lower, while the stability conditions of higher-order explicit FD schemes are more stringent [26,27]. The compact finite difference (CFD) schemes use a linear combination of function values to represent the derivative of the function, which can effectively improve accuracy and stability. Early in the 1970s, several compact Hermitian FD schemes were proposed to solve partial differential equations in the field of fluid mechanics [28,29,30]. In [31], the author introduced a family of CFD schemes with spectral-like resolution and went through a careful analysis of them. Mahesh [32] proposed a more general version of the standard CFD schemes summarized by Lele [31]. Later, the compact spectral scheme was used to approximate first-order and second-order derivatives in the simulation of three-dimensional turbulent flows [33]. Combined CFD schemes were developed to solve the convection–diffusion equation [34], Navier–Stokes equation [35], and advection equation [36] because they are more accurate and compact than conventional CFD schemes. The CFD schemes mentioned above are implemented on regular grids. Geodheer and Potters [37] and Shukla and Zhong [38] reported the application of the CFD schemes on non-equidistant grids. Up to now, various CFD schemes have been widely used in aerodynamics, hydrodynamics, and electromagnetics.
The CFD schemes also show great potential in solving seismic wave equations in exploration seismology. Yang [39] applied the CFD schemes to model the elastic and acoustic wave propagation in a 2D TI medium. Du et al. [40] combined the CFD scheme with a staggered-grid technique, developed a compact staggered-grid FD method, and applied it to simulate elastic wave propagation in a VTI medium. Kosloff [41] proposed a general CFD scheme that can calculate the derivative value at the current point using the function values at any number of points and used it to solve the elastic and acoustic wave equation. Subsequently, Chu and Stoffa [42] and Liu [43] applied the CFD schemes to solve seismic wave propagation in the frequency domain. Liao [44] discussed the accuracy, stability, and dispersion of a Padé approximation-based CFD scheme for the numerical simulation of a 3D acoustic wave equation. In [45], the authors proposed a new CFD scheme to solve a 3D acoustic wave equation, which has fourth-order accuracy in time and space and is simple to implement.
Similar to explicit FD schemes, the CFD schemes also suffer from grid dispersion, especially when large time and spatial steps are involved. Increasing the order of a specific CFD scheme can suppress the dispersion to some extent. However, higher-order CFD schemes require more computing and storage costs. To decrease the dispersion error at a large wavenumber range, optimization-based strategies are usually adopted to obtain modified FD coefficients. Kim and Lee [46] applied the dispersion–relation–preserving (DRP) method proposed by Tam and Webb [47] to minimize the dispersive errors in the wavenumber domain and obtained the optimum coefficients of tridiagonal and pentadiagonal CFD schemes. Liu et al. [48] optimized the pentadiagonal compact scheme using a sequential quadratic programming method and demonstrated its increased performance. Yu [49] proposed an optimized DRP-combined CFD scheme to solve the advection equation. Based on the least-squares method, Venutelli [50] developed two optimized fourth-order CFD schemes and presented classical applications for 1D and 2D nonlinear shallow water equations.
In this paper, we will investigate the validity of a family of tridiagonal CFD schemes and its optimization schemes for solving the 2D acoustic wave equation in an isotropic medium. The rest of this paper is as follows. In Section 2, the new CFD schemes with optimized coefficients are derived, which is followed by the dispersion and stability analysis. In Section 3, two numerical examples are used to demonstrate the feasibility and effectiveness of the proposed method, which is followed by a discussion and a conclusion in Section 4 and Section 5, respectively.

2. Theory and Methods

2.1. Tridiagonal Compact Finite-Difference Schemes for 2D Acoustic Wave Equation

Consider the 2D acoustic wave equation in an isotropic medium:
1 v 2 2 u t 2 = 2 u x 2 + 2 u z 2 + s ( t ) = Δ u + s ( t )
where u ( x , y , z ) is the scalar acoustic wavefield, v ( x , y , z ) is the velocity field, s ( t ) denotes the source function, and Δ denotes the Laplace operator.
The key to solving Equation (1) is to approximate Δ u using high-order FD schemes. For simplicity, let f ( x ) be a function of one variable, x i be the i-th grid point in the x-direction, and f i be the value of f ( x i ) . In [31], a multi-parameter family of pentadiagonal CFD schemes was proposed to approximate the second derivative of f ( x ) :
β f i 2 + α f i 1 + f i + α f i + 1 + β f i + 2 = m = 1 M a m f i + m 2 f i + f i m m 2 h x 2
where M { 1 , 2 , 3 } , h x denotes the grid interval in the x-direction, and α , β and a m are difference coefficients. By matching the Taylor series coefficients of different orders, we obtain the relations between α and a m as follows:
a 1 + a 2 + a 3 = 1 + 2 α + 2 β ( second   order )
a 1 + 2 2 a 2 + 3 2 a 3 = 4 ! 2 ! ( α + 2 2 β ) ( fourth   order )
a 1 + 2 4 a 2 + 3 4 a 3 = 6 ! 4 ! ( α + 2 4 β ) ( sixth   order )
a 1 + 2 6 a 2 + 3 6 a 3 = 8 ! 6 ! ( α + 2 6 β ) ( eighth   order )
a 1 + 2 8 a 2 + 3 8 a 3 = 10 ! 8 ! ( α + 2 8 β ) ( tenth   order )
According to Equations (3)–(7), only the tenth-order compact scheme has unique FD coefficients since there are five coefficients to be determined. Other lower-order compact schemes are not uniquely determined because of insufficient constraint conditions.
For β = 0 , a family of tridiagonal schemes is generated:
α f i 1 + f i + α f i + 1 = m = 1 M a m f i + m 2 f i + f i m m 2 h x 2
when M = 1 and Equations (3) and (4) are chosen as constraints, the classical Padé scheme with fourth-order accuracy (CFD4) is recovered. If M = 2 , we can obtain a sixth-order compact (CFD6) scheme constrained by Equations (3)–(5). For M = 3 , an eighth-order compact (CFD8) scheme constrained by Equations (3)–(6) is obtained, which is the highest-order scheme that Equation (8) can achieve. The coefficients of these schemes are tabulated in Table 1. If a further choice of α = 0 is made, conventional central FD schemes are obtained. We focus on the ability of Equation (8) to approximate Δ u in this paper.
Equation (8) can be solved using the Thomas (chase) algorithm because the coefficient matrix consisting of α is tridiagonal. The computational effort of the Thomas algorithm for solving linear equations is (5X-4), where X is the number of grid points.

2.2. Optimization of Compact Finite-Difference Coefficients

By relaxing the constraints of each CFD scheme, we can obtain corresponding low-order schemes, which provides the possibility to improve the difference coefficients using various optimization methods. For example, if the CFD4 scheme is only restricted to Equation (3), it will degenerate into a new scheme with second-order accuracy.
Fourier analysis is usually used to measure the accuracy of the FD schemes. Let w denote the product of the scaled true wavenumber k and the sampling interval h . Applying a Fourier transform to Equation (8), we obtain the scaled modified (numerical) wavenumber as follows:
L ( w ) = 2 m = 1 M a m [ 1 cos ( m w ) ] 1 + 2 α cos ( w )
where L ( w ) denotes the product of the scaled numerical wavenumber and the sampling interval h . We define a weighted deviation (integrated error) that aims to minimize w and L ( w ) :
E ( α , a 1 , a 2 , a 3 ) = 0 t ( L 2 ( w ) w 2 ) 2 G 2 ( w ) d w
where E denotes the error between w and L ( w ) , 0 < t π is the upper limit of the integral function, and G ( w ) denotes the weighting function. The weighting function makes Equation (10) integrable analytically and reduces the numerical dispersion in the high-wavenumber ranges by weighting the integrated error. The following weighting function is adopted:
G ( w ) = 1 + 2 α cos ( w )
The optimization condition that makes E a local minimum is as follows:
E α = 0
Equations (3)–(6) and (12) provide a system of linear equations by which the optimization coefficients of the tridiagonal compact schemes can be determined. Table 2 shows the optimization schemes (OCFD) with different upper limits of integral functions.
For comparison, as shown in Figure 1, we plot the difference error curves of these CFD and OCFD schemes in which the square of the numerical wavenumber varies with the true wavenumber. Figure 1a–c suggest that when the upper limit of the integral function is 0.5π or 0.75π, the curves of the OCFD schemes are closer to the exact curve than that of the CFD scheme. When the upper limit of the integral function is π, over-optimization occurs in a certain wavenumber range. In addition, as the upper limit of the integral function decreases, the accuracy curve of the OCFD scheme gradually approaches that of the CFD scheme. In fact, the OCFD scheme will degenerate into the CFD scheme if the upper limit of the integral function infinitely approaches zero. Figure 1d displays a comparison of the CFD schemes and standard central FD schemes. We observe that the accuracy of the CFD schemes is much higher than that of the FD schemes of the same order.

2.3. Dispersion Analysis and Stability Analysis

Numerical dispersion is a significant issue in FD forward modeling because the discretized form of the seismic wave equations is usually dispersive [51,52,53]. The dispersion can be measured by the error of the phase velocity v F D with respect to the true velocity v T R . Let u i , j n denote the numerical solution of the acoustic wavefield at time level n τ and grid point ( x i , z j ) . The numerical approximation of Equation (1) using second-order central finite-difference in time and various-orders compact finite-difference in space is written as:
u i , j n + 1 = 2 u i , j n u i , j n 1 + τ 2 v 2 [ ( u x x ) i , j n + ( u z z ) i , j n ]
where τ denotes the time step. The spatial derivatives u x x and u z z are solved by Equation (8). Based on the plane-wave theory, we let:
u i , j n = exp [ I ( i k x h x + j k z h z ω ˜ n τ ) ]
where θ is the inclination angle between the x-axis and the direction of plane wave propagation, the vector wavenumber k = ( k x , k z ) = ( | k | sin θ , | k | cos θ ) , ω ˜ = v F D | k | = v F D k , and I = 1 .
Without loss of generality, let h = h x = h z . We obtain the dispersion relation of the CFD schemes by substituting Equation (14) into Equation (13):
cos ( ω ˜ τ ) = 1 r 2 { m = 1 M a m [ 1 cos ( m h k x ) ] 2 α cos ( h k x ) + 1 + m = 1 M a m [ 1 cos ( m h k z ) ] 2 α cos ( h k z ) + 1 }
where r = v T R τ h . The phase velocity error is defined as follows:
γ = v F D v T R 1 = k v F D τ k v T R τ 1 = ω ˜ τ k v T R τ h h 1 = ω ˜ τ k r h 1 = arccos ( cos ( ω ˜ τ ) ) k r h 1
If γ equals 0, the CFD scheme is non-dispersive. If γ is negative and far from 0, strong spatial dispersion will occur. When γ is positive, the temporal dispersion of the scheme is large. According to Equations (15) and (16), the phase error γ depends on θ and r . Next, we fix one of them and calculate γ for k h ranging from 0 to π to compare the dispersion curves for the other parameter.
First, we study the effect of θ on the dispersion. In Figure 2, we plot γ vs. k h for different θ , while r is fixed at 0.1. From Figure 2, when r and θ are fixed, the phase velocity error of the higher-order schemes is much smaller than that of the lower-order schemes. From Figure 2a–f or Figure 2g–i, we observe that the spatial dispersion of these schemes decreases as θ increases in the high wavenumber ranges. In addition, the OCFD schemes show great potential in eliminating the spatial dispersion in the high wavenumber ranges as the upper limit of the integral function increases. However, they are not suitable for suppressing the temporal dispersion, especially when the upper limit of the integral function is π. When the upper limit of the integral function is 0.5π, we obtain three relatively good optimization schemes: OCFD4 (0.5π), OCFD6 (0.5π), and OCFD8 (0.5π). The temporal error of these optimization schemes is smaller than that of other optimization schemes.
Next, we fix the other parameter at θ = 0 to study the effect of r . Figure 2a,d,g and Figure 3 show that all the CFD and OCFD schemes are sensitive to r , which indicates that using a small time step may help to reduce the temporal dispersion error when the spatial grid interval is fixed.
Stability analysis is the other key issue in measuring the effect of an FD scheme. Almost all of the FD methods for approximating the seismic wave equations are conditionally stable and subject to constraints on velocity, spatial grid interval, and time step. The stability condition of the CFD schemes mentioned above can be derived by the Fourier method:
v max τ 1 h x 2 + 1 h z 2 2 ( 1 2 α ) m = 1 M a m [ 1 ( 1 ) m ] = C
where v max denotes the maximum velocity of the model, and C is the Courant–Friedrichs–Lewy (CFL) number. Table 3 shows the CFL number of the CFD and OCFD schemes. We observe that the CFL number of the lower-order CFD schemes is larger than that of the higher-order schemes. It should also be noted that the stability condition of the OCFD schemes is slightly stricter than that of the general scheme.

3. Numerical Examples

3.1. Homogeneous Model

We use a simple homogeneous model to evaluate the simulation accuracy of the new schemes shown in Equation (13). The model is discretized into 200 × 200 grids. The velocity of this model is 3000 m/s. A Ricker wavelet with 30 Hz dominant frequency is used for the simulation, which is located at the center of the model. The time step is 1 ms and the spatial grid interval is 20 m. Figure 4 displays the snapshots at 0.5 s simulated by the CFD4 scheme (Figure 4a), OCFD4 schemes (Figure 4b–d), and conventional central FD schemes (Figure 4e,f). In Figure 4a, the spatial grid dispersion is strong because the spatial grid interval is large. In Figure 4b–d, the spatial dispersion is gradually eliminated as the upper limit of the integral function increases (indicated by the red arrows). However, when the upper limit of the integral function is π, the time error generated in acoustic wave propagation is too large to be ignored (indicated by the black arrow).
We extract four traces from the snapshots simulated by the CFD4 and OCFD schemes for further analysis, as shown in Figure 5a. The dashed arrows indicate the spatial numerical dispersion, and the dashed box indicates the temporal dispersion. The single trace waveform of the CFD4 and OCFD4 (π) schemes between 0.4 and 0.8 km is clearly different from that of other optimization schemes, which shows that increasing the upper limit of the integral function contributes less to temporal simulation accuracy. To prove that the simulation accuracy of the proposed scheme is better than that of conventional central FD schemes, we extract single traces at a lateral distance of 2 km from the snapshots shown in Figure 4b,e,f. The result is shown in Figure 5b. From Figure 5b, we conclude that the modeling accuracy of the OCFD4 (0.5π) scheme is not inferior to that of the sixth-order central FD scheme.
The snapshots simulated by the CFD6 and CFD8 schemes and their optimization schemes are shown in Appendix A. The simulation results displayed in Figure A1 and Figure A2 are similar to those shown in Figure 4. Table 4 displays the computing time for the CFD schemes of different orders. The result shows that the computational cost of the CFD4 scheme is much lower than that of the higher-order schemes.

3.2. Marmousi Model

We use a Marmousi model to further illustrate the accuracy and effectiveness of the OCFD4 schemes in modeling acoustic wave propagation. The model is shown in Figure 6, which is discretized with 460 × 250 grids. The spatial grid spacing is 10 m and the time step is 0.5 ms. The source is located at (2300, 10) m, which is a Ricker wavelet with 30 Hz dominant frequency. A perfectly matched layer absorbing boundary [54,55] is used for wavefield extrapolation.
Figure 7 displays the wavefield snapshots at 0.9 s computed by the CFD4 (Figure 7a) and OCFD4 (Figure 7c–e) schemes. The snapshot in Figure 7b is simulated by the CFD6 scheme, which is used as a reference. When the upper limit of the integral function is 0.5π or 0.75π, the wavefield is calculated correctly. As indicated by the red arrows in Figure 7e, the wavefield modeled by the OCFD4 (π) scheme is different from those simulated by other schemes. Figure 8 shows the difference in the simulated wavefield. The error of the residual wavefield in Figure 8a–d relative to the reference wavefield in Figure 7b at L2-norm is 1.65%, 1.06%, 1.25%, and 1.74%, respectively. Obviously, the OCFD4 (0.5π) scheme is the most helpful in improving the simulation accuracy of the Marmousi model. Figure 9 displays the shot records simulated by the CFD4, CFD6, and OCFD4 (π) schemes, respectively. We extract two traces from the records for further comparison, as shown in Figure 10. The result shows that the simulation accuracy of the OCFD4 (π) scheme is better than that of the CFD4 scheme. From Figure 7, Figure 8, Figure 9 and Figure 10, we conclude that the OCFD4 (π) scheme is effective, and it can be used to obtain more accurate simulation results than the general CFD4 scheme.

4. Discussion

Compared to CFD schemes, most explicit FD schemes often used in seismic forward modeling are not compact, which gives rise to difficulty in dealing with boundary conditions. For example, a sixth-order central FD scheme requires a 13-points stencil and three layers of boundary conditions for approximating the second derivative in 2D cases while the CFD4 scheme with similar approximation accuracy only requires a 5-point stencil and one layer of boundary conditions. In 3D cases, the CFD4 scheme requires a 7-points stencil and one layer of boundary conditions, while the sixth-order central FD scheme requires a 19-point stencil and three layers of boundary conditions. Therefore, the CFD schemes require fewer boundary conditions than the explicit FD schemes. This advantage of the CFD schemes also helps to save some storage space in numerical simulation.
We studied three CFD schemes and their corresponding optimization schemes in this paper. In fact, the other two CFD schemes have more than one optimization scheme except for the CFD4 scheme. For example, if the CFD6 scheme is only constrained by Equation (3), a second-order OCFD6 scheme is generated. In this case, two optimization conditions, such as Equation (12), must be solved to make Equation (10) a local minimum. For the second-order OCFD8 scheme, it is necessary to solve a system of equations consisting of three optimization conditions to obtain the optimization coefficients. All possible tridiagonal optimization schemes are summarized in Table A1 of Appendix B. The coefficients for the OCFD6_2, OCFD8_2, and OCFD8_4 schemes are not given in the paper due to the limitation of computational power. In addition, Equation (10) can be regarded as a least square problem. Therefore, some optimization algorithms, such as the minimax approximation method [56], Newton method [57], and conjugate gradient method [58], can be used to obtain these coefficients.
The selection of the upper limit of the integral function has a significant impact on the OCFD schemes. For each optimization scheme, we only discussed three representative cases where the upper limit of the integral is 0.5π, 0.75π, and π, respectively. The coefficients in these cases are easily obtained by solving Equation (12) because some terms cancel out in the calculation. The weighting function is also important for optimization results. Some exponential weighting functions [46,59] can help to obtain better optimization coefficients by weighting the integrated error in the high wave number range. However, it makes it difficult to find the primitive function of Equation (10).

5. Conclusions

We have proposed a high-accuracy forward modeling scheme for simulating acoustic wave propagation. A family of tridiagonal OCFD schemes with 2M-order accuracy in space and the central FD scheme with second-order accuracy in time are used to approximate the spatial and temporal second derivatives, respectively. The accuracy curves show that the CFD schemes are superior to conventional FD schemes of the same order, and the OCFD schemes are better than the CFD schemes in a certain wavenumber range. The dispersion analysis demonstrates that increasing the upper limit of the integral function helps to reduce spatial error but contributes less to improving temporal dispersion. In addition, the wavefield simulations on a homogeneous model confirm the conclusion drawn in the dispersion analysis. The efficiency tests for the homogeneous model show that the lower-order CFD scheme is more effective than the higher-order schemes. The numerical experiments on the Marmousi model suggest that the OCFD4 (0.5π and 0.75π) schemes facilitate improvement in the wavefield modeling of complex structures.

Author Contributions

Conceptualization, L.C. and J.H. (Jianping Huang); formal analysis, W.P.; investigation, C.S. and J.H. (Jiale Han); methodology, L.C.; resources, J.H. (Jianping Huang) and L.-Y.F.; software, C.S. and J.H. (Jiale Han); supervision, J.H. (Jianping Huang); writing—original draft, L.C.; writing—review and editing, J.H. (Jianping Huang). All authors have read and agreed to the published version of the manuscript.

Funding

This research is supported by the Marine S&T Fund of Shandong Province for Pilot National Laboratory for Marine Science and Technology (Qingdao) (grant number 2021QNLM020001), the National Key R&D Program of China (grant number 2019YFC0605503C), the National Outstanding Youth Science Foundation (grant number 41922028), and the Major Scientific and Technological Projects of China National Petroleum Corporation (grant number ZD2019-183-003).

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1 shows the snapshots of the homogeneous model simulated by the CFD6 and its optimization schemes. Figure A2 shows the snapshots of the homogeneous model simulated by the CFD8 and its optimization schemes. The two figures are shown below:
Figure A1. Snapshots of the homogeneous model at 0.5 s. They are simulated by: (a) CFD6 scheme, (b) OCFD6 scheme (0.5π), (c) OCFD6 scheme (0.75π), and (d) OCFD6 scheme (π), respectively. The red arrows indicate the spatial dispersion and the black arrow indicates the temporal dispersion.
Figure A1. Snapshots of the homogeneous model at 0.5 s. They are simulated by: (a) CFD6 scheme, (b) OCFD6 scheme (0.5π), (c) OCFD6 scheme (0.75π), and (d) OCFD6 scheme (π), respectively. The red arrows indicate the spatial dispersion and the black arrow indicates the temporal dispersion.
Remotesensing 15 00604 g0a1aRemotesensing 15 00604 g0a1b
Figure A2. Snapshots of the homogeneous model at 0.5 s. They are simulated by: (a) CFD8 scheme, (b) OCFD8 scheme (0.5π), (c) OCFD8 scheme (0.75π), and (d) OCFD8 scheme (π), respectively. The red arrows indicate spatial dispersion and the black arrow indicates the temporal dispersion.
Figure A2. Snapshots of the homogeneous model at 0.5 s. They are simulated by: (a) CFD8 scheme, (b) OCFD8 scheme (0.5π), (c) OCFD8 scheme (0.75π), and (d) OCFD8 scheme (π), respectively. The red arrows indicate spatial dispersion and the black arrow indicates the temporal dispersion.
Remotesensing 15 00604 g0a2

Appendix B

Table A1 shows all tridiagonal optimization schemes with different orders. OCFD6_2 denotes the second-order OCFD6 scheme, OCFD8_4 denotes the fourth-order OCFD8 scheme, and OCFD8_2 denotes the second-order OCFD8 scheme.
Table A1. All tridiagonal optimization schemes with different orders.
Table A1. All tridiagonal optimization schemes with different orders.
SchemesConstraintsOrderOptimization Conditions
CFD4Equations (3) and (4)fourth/
OCFD4Equation (3)second E / α = 0
CFD6Equations (3)–(5)sixth/
OCFD6Equations (3) and (4)fourth E / α = 0
OCFD6_2Equation (3)second E / α = 0 , E / a 1 = 0
CFD8Equations (3)–(6)eighth/
OCFD8Equations (3)–(5)sixth E / α = 0
OCFD8_4Equations (3) and (4)fourth E / α = 0 , E / a 1 = 0
OCFD8_2Equation (3)second E / α = 0 , E / a 1 = 0 , E / a 2 = 0

References

  1. Li, X.; Yao, G.; Niu, F.; Wu, D.; Liu, N. Waveform inversion of seismic first arrivals acquired on irregular surface. Geophysics 2022, 87, 1MJ-V246. [Google Scholar] [CrossRef]
  2. Yao, G.; Wu, B.; da Silva, N.V.; Debens, H.A.; Wu, D.; Cao, J. Least-squares reverse-time migration with a multiplicative Cauchy constraint. Geophysics 2022, 87, 1MJ-V246. [Google Scholar] [CrossRef]
  3. Han, Y.; Wu, B.; Yao, G.; Ma, X.; Wu, D. Eliminate time dispersion of seismic wavefield simulation with semi-supervised deep learning. Energies 2022, 15, 7701. [Google Scholar] [CrossRef]
  4. Baysal, E.; Kosloff, D.D.; Sherwood, J.W. Reverse time migration. Geophysics 1983, 48, 1514–1524. [Google Scholar] [CrossRef] [Green Version]
  5. Abubakar, A.; Pan, G.; Li, M.; Zhang, L.; Habashy, T.M.; Berg, P.M.V.D. Three-dimensional seismic full-waveform inversion using the finite-difference contrast source inversion method. Geophys. Prospect. 2011, 59, 874–888. [Google Scholar] [CrossRef]
  6. Alford, R.M.; Kelly, K.R.; Boore, D.M. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics 1974, 39, 834–842. [Google Scholar] [CrossRef] [Green Version]
  7. Kelly, K.R.; Ward, R.W.; Treitel, S.; Alford, R.M. Synthetic seismograms: A finite difference approach. Geophysics 1976, 41, 2–27. [Google Scholar] [CrossRef]
  8. Virieux, J. SH-wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics 1984, 49, 1933–1942. [Google Scholar] [CrossRef]
  9. Virieux, J. P-SV wave propagation in heterogeneous media velocity-stress finite-difference method. Geophysics 1984, 51, 889–901. [Google Scholar] [CrossRef]
  10. Marfurt, K.J. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics 1984, 49, 533–549. [Google Scholar] [CrossRef]
  11. Kay, I.; Krebes, E.S. Applying finite-element analysis to the memory variable formulation of wave propagation in anelastic media. Geophysics 1999, 64, 300–307. [Google Scholar] [CrossRef]
  12. Komatitsch, D.; Vilotte, J.P. The spectral-element method: An efficient tool to simulate the seismic response of 2-D and 3-D geological structures. Bull. Seismol. Soc. Am. 1998, 88, 368–392. [Google Scholar] [CrossRef]
  13. Komatitsch, D.; Tromp, J. Introduction to the spectral-element method for 3-D seismic wave propagation. Geophys. J. Int. 1999, 139, 806–822. [Google Scholar] [CrossRef] [Green Version]
  14. Komatitsch, D.; Tromp, J. Spectral-element simulations of global seismic wave propagation—Part I: Validation. Geophys. J. Int. 2002, 149, 390–412. [Google Scholar] [CrossRef]
  15. Komatitsch, D.; Tromp, J. Spectral-element simulations of global seismic wave propagation—Part II: 3-D models, oceans, rotation, and self-gravitation. Geophys. J. Int. 2002, 150, 303–318. [Google Scholar] [CrossRef]
  16. Sánchez-Sesma, F.J.; Esquivel, J. Ground motion on alluvial valleys under incident plane SH waves. Bull. Seismol. Soc. Am. 1979, 69, 1107–1120. [Google Scholar] [CrossRef]
  17. Gaffet, S.; Bouchon, M. Source location and valley shape effects on the P-SV displacement field using a boundary integral equation-discrete wavenumber representation method. Geophys. J. Int. 1991, 106, 341–355. [Google Scholar] [CrossRef] [Green Version]
  18. Sánchez-Sesma, F.J.; Campillo, M. Diffraction of P, SV, and Raileigh waves by topographic features: A boundary integral formulation. Bull. Seismol. Soc. Am. 1991, 81, 2234–2253. [Google Scholar] [CrossRef]
  19. Bouchon, M.; Schultz, C.A.; Toksoz, M.N. A fast implementation of boundary integral equation methods to calculate the propagation of seismic waves in laterally varying layered medium. Bull. Seismol. Soc. Am. 1995, 85, 1679–1687. [Google Scholar] [CrossRef]
  20. Fornberg, B. The pseudospectral method: Comparisons with finite differences for the elastic wave equation. Geophysics 1987, 52, 483–501. [Google Scholar] [CrossRef]
  21. Fornberg, B. The pseudospectral method: Accurate representation of interfaces in elastic wave calculations. Geophysics 1988, 53, 625–637. [Google Scholar] [CrossRef]
  22. Faccioli, E.; Maggio, F.; Paolucci, R.; Quarteroni, A. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method. J. Seismol. 1997, 1, 237–251. [Google Scholar] [CrossRef]
  23. Dormy, E.; Tarantola, A. Numerical simulation of elastic wave propagation using a finite volume method. J. Geophys. Res. 1995, 100, 2123–2133. [Google Scholar] [CrossRef]
  24. Robertsson, J.O.A. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography. Geophysics 1996, 61, 1921–1934. [Google Scholar] [CrossRef]
  25. Das, S.; Liao, W.; Gupta, A. An efficient fourth-order low dispersive finite difference scheme for a 2-D acoustic wave equation. J. Comput. Appl. Math. 2014, 258, 151–167. [Google Scholar] [CrossRef]
  26. Ren, Y.; Huang, J.; Yong, P.; Liu, M.; Cui, C.; Yang, M. Optimized staggered-grid finite-difference operators using window functions. Appl. Geophys. 2018, 15, 253–260. [Google Scholar] [CrossRef]
  27. Yang, D.H.; Wang, N.; Chen, S.; Song, G.J. An explicit method based on the implicit Runge–Kutta algorithm for solving the wave equations. Bull. Seismol. Soc. Am. 2009, 99, 3340–3354. [Google Scholar] [CrossRef]
  28. Orszag, S.A.; Israeli, M. Numerical simulation of viscous incompressible flows. Annu. Rev. Fluid Mech. 1974, 6, 281–318. [Google Scholar] [CrossRef]
  29. Hirsh, R.S. Higher order accurate difference solutions of fluid mechanics problems by a compact differencing technique. J. Comput. Phys. 1975, 19, 90–109. [Google Scholar] [CrossRef]
  30. Adam, Y. A Hermitian finite difference method for the solution of parabolic equations. Comput. Math. Appl. 1975, 1, 393–406. [Google Scholar] [CrossRef]
  31. Lele, S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 1992, 103, 16–42. [Google Scholar] [CrossRef]
  32. Mahesh, K. A family of high order finite difference schemes with good spectral resolution. J. Comput. Phys. 1998, 145, 332–358. [Google Scholar] [CrossRef] [Green Version]
  33. Lee, C.; Seo, Y. A new compact spectral scheme for turbulence simulations. J. Comput. Phys. 2002, 183, 438–469. [Google Scholar] [CrossRef]
  34. Chu, P.C.; Fan, C. A three-point combined compact difference scheme. J. Comput. Phys. 1998, 140, 370–399. [Google Scholar] [CrossRef] [Green Version]
  35. Sengupta, T.K.; Lakshmanan, V.; Vijay, V.V.S.N. A new combined stable and dispersion relation preserving compact scheme for non-periodic problems. J. Comput. Phys. 2009, 228, 3048–3071. [Google Scholar] [CrossRef]
  36. Mohebbi, A.; Abbaszadeh, M. Compact finite difference scheme for the solution of time fractional advection-dispersion equation. Numer Algorithms. 2013, 63, 431–452. [Google Scholar] [CrossRef]
  37. Geodheer, W.J.; Potters, J.H.H.M. A compact finite difference scheme on a non-equidistant mesh. J. Comput. Phys. 1985, 61, 269–279. [Google Scholar] [CrossRef]
  38. Shukla, R.K.; Zhong, X. Derivation of high-order compact finite difference schemes for non-uniform grid using polynomial interpolation. J. Comput. Phys. 2005, 204, 404–429. [Google Scholar] [CrossRef]
  39. Yang, D.; Wang, S.; Zhang, Z.; Teng, J. N-times absorbing boundary conditions for compact finite-difference modeling of acoustic and elastic wave propagation in the 2D TI medium. Bull. Seismol. Soc. Am. 2003, 93, 2389–2401. [Google Scholar] [CrossRef]
  40. Du, Q.; Li, B.; Hou, B. Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme. Appl. Geophys. 2009, 6, 42–49. [Google Scholar] [CrossRef]
  41. Kosloff, D.; Pestana, R.C.; Tal-Ezer, H. Acoustic and elastic numerical wave simulations by recursive spatial derivative operators. Geophysics 2010, 75, T167–T174. [Google Scholar] [CrossRef]
  42. Chu, C.; Stoffa, P.L. Frequency domain modeling using implicit spatial finite difference operators. In SEG Technical Program Expanded Abstracts 2010; Society of Exploration Geophysicists: Houston, TX, USA, 2010; pp. 3076–3080. [Google Scholar] [CrossRef]
  43. Liu, Y. An optimal 5-point scheme for frequency-domain scalar wave equation. J. Appl. Geophys. 2014, 108, 19–24. [Google Scholar] [CrossRef]
  44. Liao, W. On the dispersion, stability and accuracy of a compact higher-order finite difference scheme for 3D acoustic wave equation. J. Comput. Appl. Math. 2014, 270, 571–583. [Google Scholar] [CrossRef]
  45. Li, K.; Liao, W. An efficient and high accuracy finite-difference scheme for the acoustic wave equation in 3D heterogeneous media. J. Comput. Sci. 2020, 40, 101063. [Google Scholar] [CrossRef] [Green Version]
  46. Kim, J.W.; Lee, D.J. Optimized compact finite difference schemes with maximum resolution. AIAA 1996, 34, 887–893. [Google Scholar] [CrossRef] [Green Version]
  47. Tam, C.K.W.; Webb, J.C. Dispersion–relation–preserving schemes for computational aeroacoustics. AIAA 1992, 1, 92-02-033. [Google Scholar]
  48. Liu, Z.; Huang, Q.; Zhao, Z.; Yuan, J. Optimized compact finite difference schemes with high accuracy and maximum resolution. Int. J. Aeroacoustics 2008, 7, 123–146. [Google Scholar] [CrossRef]
  49. Yu, C.H.; Wang, D.; He, Z.; Pähtz, T. An optimized dispersion-relation-preserving combined compact difference scheme to solve advection equations. J. Comput. Phys. 2015, 300, 92–115. [Google Scholar] [CrossRef] [Green Version]
  50. Venutelli, M. New optimized fourth-order compact finite difference schemes for wave propagation phenomena. Appl. Numer. Math. 2014, 87, 53–73. [Google Scholar] [CrossRef]
  51. Huang, J.; Yang, J.; Liao, W.; Wang, X.; Li, Z. Common-shot Fresnel beam migration based on wave-field approximation in effective vicinity under complex topographic conditions. Geophys. Prospect. 2016, 64, 554–570. [Google Scholar] [CrossRef]
  52. Zhang, Y.; Fu, L.-Y.; Zhang, L.X.; Wei, W.; Guan, X. Finite difference modeling of ultrasonic propagation (coda waves) in digital porous cores with un-split convolutional PML and rotated staggered grid. J. Appl. Geophys. 2014, 104, 75–89. [Google Scholar] [CrossRef]
  53. Hou, W.; Fu, L.-Y.; Carcione, J.M.; Wang, Z.W.; Wei, J. Simulation of thermoelastic waves based on the Lord-Shulman theory. Geophysics 2021, 86, T155–T164. [Google Scholar] [CrossRef]
  54. Zhang, L.X.; Fu, L.-Y.; Pei, Z.L. Finite difference modeling of Biot’s poroelastic equations with unsplit convolutional PML and rotated staggered grid. Chin. J. Geophys. 2010, 53, 2470–2483. [Google Scholar] [CrossRef]
  55. Yang, H.; Fu, L.-Y.; Fu, B.-Y.; Du, Q. Poro-acoustoelasticity finite-difference simulation of elastic wave propagation in prestressed porous media. Geophysics 2022, 87, T329–T345. [Google Scholar] [CrossRef]
  56. Yang, L.; Yan, H.Y.; Liu, H. Optimal staggered-grid finite-difference schemes based on the minimax approximation method with the Remez algorithm. Geophysics 2017, 82, T27–T42. [Google Scholar] [CrossRef]
  57. Yong, P.; Huang, J.P.; Li, Z.C.; Liao, W.Y.; Qu, L.P.; Li, Q.Y.; Liu, P.J. Optimized equivalent staggered-grid FD method for elastic wave modelling based on plane wave solutions. Geophys. J. Int. 2017, 208, 1157–1172. [Google Scholar] [CrossRef]
  58. Jin, K.; Huang, J.P.; Zou, Q.; Wang, Z.; Tong, S.; Liu, B.; Hu, Z. Optimization of staggered grid finite-difference coefficients based on conjugate gradient method. J. Seism. Explor. 2021, 31, 33–52. [Google Scholar]
  59. Haras, Z.; Ta’asan, S. Finite difference schemes for long–time integration. J. Comput. Phys. 1994, 114, 265–279. [Google Scholar] [CrossRef]
Figure 1. Accuracy curves: (a) comparison of the CFD4 and OCFD4 schemes, (b) comparison of the CFD6 and OCFD6 schemes, (c) comparison of the CFD8 and OCFD8 schemes, (d) comparison of the central FD and CFD schemes.
Figure 1. Accuracy curves: (a) comparison of the CFD4 and OCFD4 schemes, (b) comparison of the CFD6 and OCFD6 schemes, (c) comparison of the CFD8 and OCFD8 schemes, (d) comparison of the central FD and CFD schemes.
Remotesensing 15 00604 g001
Figure 2. Dispersion curves of various θ : (ac) the CFD4 scheme and its optimization schemes, (df) the CFD6 scheme and its optimization schemes, (gi) the CFD8 scheme and its optimization schemes.
Figure 2. Dispersion curves of various θ : (ac) the CFD4 scheme and its optimization schemes, (df) the CFD6 scheme and its optimization schemes, (gi) the CFD8 scheme and its optimization schemes.
Remotesensing 15 00604 g002
Figure 3. Dispersion curves of various r . (a,b) the CFD4 scheme and its optimization schemes, (c,d) the CFD6 scheme and its optimization schemes, (e,f) the CFD8 scheme and its optimization schemes.
Figure 3. Dispersion curves of various r . (a,b) the CFD4 scheme and its optimization schemes, (c,d) the CFD6 scheme and its optimization schemes, (e,f) the CFD8 scheme and its optimization schemes.
Remotesensing 15 00604 g003
Figure 4. Snapshots of the homogeneous model at 0.5 s. They are simulated by: (a) CFD4 scheme, (b) OCFD4 scheme (0.5π), (c) OCFD4 scheme (0.75π), (d) OCFD4 scheme (π), (e) fourth-order central FD scheme, and (f) sixth-order central FD scheme, respectively. The red arrows indicate the spatial dispersion and the black arrow indicates the temporal dispersion.
Figure 4. Snapshots of the homogeneous model at 0.5 s. They are simulated by: (a) CFD4 scheme, (b) OCFD4 scheme (0.5π), (c) OCFD4 scheme (0.75π), (d) OCFD4 scheme (π), (e) fourth-order central FD scheme, and (f) sixth-order central FD scheme, respectively. The red arrows indicate the spatial dispersion and the black arrow indicates the temporal dispersion.
Remotesensing 15 00604 g004aRemotesensing 15 00604 g004b
Figure 5. Single traces of the snapshots at lateral distance of 2 km. (a) Comparison of the CFD4 and OCFD4 schemes and (b) comparison of the OCFD4 (0.5π) and central FD schemes. The dashed box indicates the temporal dispersion.
Figure 5. Single traces of the snapshots at lateral distance of 2 km. (a) Comparison of the CFD4 and OCFD4 schemes and (b) comparison of the OCFD4 (0.5π) and central FD schemes. The dashed box indicates the temporal dispersion.
Remotesensing 15 00604 g005aRemotesensing 15 00604 g005b
Figure 6. Marmousi model.
Figure 6. Marmousi model.
Remotesensing 15 00604 g006
Figure 7. Snapshots of the Marmousi model at 0.5 s. They are simulated by: (a) CFD4 scheme, (b) CFD6 scheme, (c) OCFD4 scheme (0.5π), (d) OCFD4 scheme (0.75π), and (e) OCFD4 scheme (π), respectively. The red arrows indicate the dispersion error.
Figure 7. Snapshots of the Marmousi model at 0.5 s. They are simulated by: (a) CFD4 scheme, (b) CFD6 scheme, (c) OCFD4 scheme (0.5π), (d) OCFD4 scheme (0.75π), and (e) OCFD4 scheme (π), respectively. The red arrows indicate the dispersion error.
Remotesensing 15 00604 g007
Figure 8. The difference of the simulated wavefield between: (a) Figure 7a,b, (b) Figure 7b,c, (c) Figure 7b,d, and (d) Figure 7b,e.
Figure 8. The difference of the simulated wavefield between: (a) Figure 7a,b, (b) Figure 7b,c, (c) Figure 7b,d, and (d) Figure 7b,e.
Remotesensing 15 00604 g008
Figure 9. Shot records simulated by the: (a) CFD4 scheme, (b) CFD6 scheme, (c) OCFD4 (0.5π).
Figure 9. Shot records simulated by the: (a) CFD4 scheme, (b) CFD6 scheme, (c) OCFD4 (0.5π).
Remotesensing 15 00604 g009
Figure 10. Single traces exacted from the shot records in Figure 9 at the distance of (a) 1.5 km and (b) 3.75 km.
Figure 10. Single traces exacted from the shot records in Figure 9 at the distance of (a) 1.5 km and (b) 3.75 km.
Remotesensing 15 00604 g010
Table 1. Coefficients of the tridiagonal compact schemes.
Table 1. Coefficients of the tridiagonal compact schemes.
SchemesConstraintsOrder α a 1 a 2 a 3
CFD4Equations (3) and (4)fourth0.11.200
CFD6Equations (3)–(5)sixth0.1818181.0909090.0618180
CFD8Equations (3)–(6)eighth0.2368420.9671050.536842−0.030263
Table 2. Coefficients of the tridiagonal compact optimization schemes.
Table 2. Coefficients of the tridiagonal compact optimization schemes.
SchemesConstraintsOrderIntegral
Limit
α a 1 a 2 a 3
OCFD4Equation
(3)
secondπ0.1660541.33210900
0.75π0.1315111.26302100
0.5π0.1125311.22506300
OCFD6Equations
(3) and (4)
fourthπ0.2773270.9635640.5910900
0.75π0.2243041.0342620.4143460
0.5π0.1980531.0692620.3268440
OCFD8Equations
(3)–(5)
sixthπ0.3325450.7517750.996214−0.082900
0.75π0.2774860.8756560.731933−0.052617
0.5π0.2519030.9332170.609136−0.038547
Table 3. CFL number of the CFD and OCFD schemes.
Table 3. CFL number of the CFD and OCFD schemes.
SchemesCFL NumberSchemesCFL NumberSchemesCFL Number
CFD40.817CFD60.764CFD80.750
OCFD4 (0.5π)0.795OCFD6 (0.5π)0.752OCFD8 (0.5π)0.745
OCFD4 (0.75π)0.764OCFD6 (0.75π)0.730OCFD8 (0.75π)0.735
OCFD4 (π)0.708OCFD6 (π)0.680OCFD8 (π)0.708
Table 4. Computing time for the CFD schemes of different orders.
Table 4. Computing time for the CFD schemes of different orders.
SchemesRecording Time (s)Computing Time (s)Recording Time (s)Computing Time (s)Recording Time (s)Computing Time (s)
CFD4536.71073.220145.6
CFD6570.910142.020283.7
CFD85108.510217.420433.1
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, L.; Huang, J.; Fu, L.-Y.; Peng, W.; Song, C.; Han, J. A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation. Remote Sens. 2023, 15, 604. https://doi.org/10.3390/rs15030604

AMA Style

Chen L, Huang J, Fu L-Y, Peng W, Song C, Han J. A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation. Remote Sensing. 2023; 15(3):604. https://doi.org/10.3390/rs15030604

Chicago/Turabian Style

Chen, Liang, Jianping Huang, Li-Yun Fu, Weiting Peng, Cheng Song, and Jiale Han. 2023. "A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation" Remote Sensing 15, no. 3: 604. https://doi.org/10.3390/rs15030604

APA Style

Chen, L., Huang, J., Fu, L. -Y., Peng, W., Song, C., & Han, J. (2023). A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation. Remote Sensing, 15(3), 604. https://doi.org/10.3390/rs15030604

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