Trapezoid-Grid Finite-Difference Time-Domain Method for 3D Seismic Wavefield Modeling Using CPML Absorbing Boundary Condition

Abstract

The large computational memory requirement is an important issue in 3D large-scale wave modeling, especially for GPU calculation. Based on the observation that wave propagation velocity tends to gradually increase with depth, we propose a 3D trapezoid-grid finite-difference time-domain (FDTD) method to achieve the reduction of memory usage without a significant increase of computational time or a decrease of modeling accuracy. It adopts the size-increasing trapezoid-grid mesh to fit the increasing trend of seismic wave velocity in depth, which can significantly reduce the oversampling in the high-velocity region. The trapezoid coordinate transformation is used to alleviate the difficulty of processing ununiform grids. We derive the 3D acoustic equation in the new trapezoid coordinate system and adopt the corresponding trapezoid-grid convolutional perfectly matched layer (CPML) absorbing boundary condition to eliminate the artificial boundary reflection. Stability analysis is given to generate stable modeling results. Numerical tests on the 3D homogenous model verify the effectiveness of our method and the trapezoid-grid CPML absorbing boundary condition, while numerical tests on the SEG/EAGE overthrust model indicate that for comparable computational time and accuracy, our method can achieve about 50% reduction on memory usage compared with those on the uniform-grid FDTD method.

1 Introduction

Reverse time migration (RTM) (Baysal et al., 1983; Xuan et al., 2014; Qu et al., 2015; Xu et al., 2021a; Du et al., 2021) and full-waveform inversion (FWI) (Tarantola, 1984; Virieux and Operto, 2009; Cai and Zhang, 2015; Xia et al., 2017; Jia et al., 2019) play a fundamental role in geophysical exploration. Since forward modeling of the wave equation consumes most computational time in the RTM and FWI processes (Jing et al., 2019; Xu et al., 2021b; Liu et al., 2021), how to achieve the improvement of efficiency and the reduction of memory usage without a significant decrease of accuracy for 3D large-scale modeling is a key problem of seismic modeling. The finite-difference (FD) method has been regarded as one of the most popular wave modeling methods for its easy implementation and high-computational efficiency (Antunes et al., 2014; Abreu et al., 2015; Xu and Gao, 2018; Robertsson and Blanch, 2020). However, the numerical dispersion of the traditional FD method leads to the use of fine grids or high-order operators (Dablain, 1986; Liu and Sen, 2011b), which inevitably affects the efficiency of simulation.

The conventional FD method literally adopts a weighted summation of neighboring grid points’ values to estimate the derivative for a designated grid point (Zhou et al., 2021), where the grid size (h) is fixed and the FD coefficients are calculated by Taylor expansion. In this way, the approximation error ϵ can be expressed as (Liu and Sen, 2011a; Wu et al., 2019b) follows:where 2M is the length of the FD operator, h is the spatial interval, f is the frequency, and v is the seismic wave velocity. Considering that λ = v/f is the wavelength and G = λ/h is the number of grid points per wavelength (NPPW), we can rewrite Eq. 1 as follows:

Eq. 2 indicates that the modeling accuracy of the conventional FD method is proportional to G and M. Because the seismic wave velocity is varying in different positions, the wavelength is short in low-velocity regions and long in high-velocity regions (Liu, 2020). Therefore, a part of computing resources is wasted in the high-velocity regions for the fixed spatial interval and FD order. With respect to this problem, there are two kinds of techniques corresponding to the different understanding of Eq. 2. The first one is the variable-operator FD method (Liu and Sen, 2011a), which adopts the long and short FD stencils in the low- and high-velocity region, respectively. For the scheme in Liu and Sen (2011a), the variable-length FD stencils are designed by approaching the dispersion relation in the time–space domain, and Liu (2020) subsequently optimizes their FD coefficients.

The second one is the variable-grid FD method, which adopts different gird sizes in different regions and can efficiently reduce the oversampling in the high-velocity region. The key problem of the variable-grid FD method is the processing of the transition area between the fine grids and the coarse grids. The variable-grid FD method based on interpolation (Hayashi and Burns, 2005; Pasalic and McGarry, 2010) is the easiest one, in which the lacking information in the transition area is completed by interpolation. However, the resulting artificial reflection in the transition area and the possible instability make it inefficient for high-accuracy seismic wave simulation. Another variable-grid FD method adopts irregular FD coefficients to process the transition region (Huang and Dong, 2009; Liu et al., 2014), which can significantly avoid the artificial reflection and improve the stability. The disadvantage of this type of variable-grid method is the additional computing cost brought by calculating irregular FD coefficients.

The trapezoid-grid FD method (Chen and Xu, 2012; Gao et al., 2018; Wu et al., 2018, Wu et al., 2019a, Wu et al., 2019b) is one of the practical variable-grid methods. It uses the trapezoid-grid mesh to fit the trend of velocity increasing with depth, which can effectively reduce the number of required grid points. Meanwhile, the use of trapezoid coordinate transformation can avoid the difficulty of processing ununiform grids in the physical Cartesian coordinate system. On the other hand, the significant reduction of memory requirement of trapezoid-grid FDTD can improve the easy implementation of GPU calculation (Fujii et al., 2013; Li et al., 2016). The existing research on trapezoid-grid FDTD methods mainly focuses on 2D wavefield modeling. Therefore, it is essential to expand trapezoid-grid FDTD from 2D to 3D for realistic seismic exploration research.

In this article, we propose a 3D trapezoid-grid FDTD method for acoustic wave modeling. First, we design the 3D trapezoid coordinate transformation and derive the 3D acoustic equation in the trapezoid coordinate system. Second, to reduce the artificial boundary reflection (Ma et al., 2018, Ma et al., 2019), we apply the corresponding trapezoid-grid convolutional perfectly matched layer (CPML) absorbing boundary condition. Third, stability analysis is given to generate stable modeling results. We then test our proposed method on the 3D homogenous model and the SEG/EAGE overthrust model and compare the efficiency and accuracy of the trapezoid-grid FDTD method with the uniform-grid FDTD method. Finally, conclusions are shown in the last section.

2 Methods

2.1 3D Trapezoid Coordinate System

In this article, the 3D trapezoid coordinate transformation is defined aswhere (x₀, y₀, z₀) is the Cartesian coordinate system, and (x, y, z) is the defined trapezoid coordinate system. In Eq. 3, α and β are central horizontal positions of the 3D trapezoid mesh, and γ is the scaling parameter for lateral coordinates. The velocity-related function g(z) is the sampling function for z₀-axis, which should be first- and second-order continuous for deriving 3D wave equations in the trapezoid coordinate system. The discrete points of g(z) are given by the following recursion:where f₀ is the dominant frequency of the source term, N₀ is the preferred NPPW and is related to the accuracy of the adopted FD scheme, and v_min(g(z)) is the selected minimum velocity at depth g(z) in the physical Cartesian coordinate system and is smoothed by solving a local polynomial fitting problem with the constraint that v_min(g(z)) should not be greater than the model minimum velocity at depth g(z). The central value of each local polynomial corresponds to a value of v_min(g(z)). In particular, we usually set the order of the local polynomial as three. Such v_min(g(z)) can lead to a smooth sampling function g(z) for discontinuous velocity variation while satisfying the required number of points per wavelength in the z₀-direction.

If the grid sizes for the trapezoid-grid FDTD in the trapezoid coordinate system are defined as Δx, Δy, and Δz, then the corresponding grid sizes in the physical Cartesian system can be described asIn our work, γ and g(z) are determined adaptively according to the model velocity. By selecting γ such that Δx₀(z) and Δy₀(z) are always smaller than or equal to Δz₀(z), and a variable-grid mesh adaptive to the velocity model can be achieved in the physical Cartesian coordinate system. Figure 1 shows the schematic of the 3D trapezoid coordinate transform. Figure 1A shows the trapezoid-grid mesh in the Cartesian coordinate system, while Figure 1B shows the corresponding uniform grid mesh in the transformed trapezoid coordinate system. In particular, the two gray regions in Figure 1 represent the same physical region.

FIGURE 1

2.2 3D Acoustic Equation with CPML Absorbing Boundary Condition in the Trapezoid Coordinate System

According to the theory of Pasalic and McGarry (2010), the time-domain-discretization form of the 3D isotropic acoustic equation with the CPML absorbing boundary condition in the Cartesian coordinate system can be described aswhere u^j = u(x₀, y₀, z₀, t_j) represents the scalar wavefield at the jth time step in the Cartesian coordinate system; v is the velocity; Δt is the time interval; f(t) is the source term; is the position of source; ln, where R denotes the designated theoretical boundary reflection coefficient, v^max is the maximum velocity of the model, is the thickness of CPML absorbing boundary along the τ₀ direction, and denotes the distance to the inner area in the τ₀ direction; and α^max = πf₀.

In order to derive the acoustic equation in the trapezoid coordinate system, we first need to transform the derivatives in the Cartesian coordinate system into the derivatives in the trapezoid coordinate system. Based on the definition of the trapezoid coordinate system in Eq. 3 and the derivation rule of the composite function, the relationships of first- and second-order derivatives in the two coordinate systems can be given as

Since the mixed spatial derivatives in Eq. 7f are hard to discrete directly with the FD method, to transform the mixed spatial derivatives into non-mixed spatial derivatives, we define three rotation transformations in the trapezoid coordinate system aswhere and are axes along diagonals in the (x, z) planes, and are axes along diagonals in the (y, z) planes, and and are axes along diagonals in the (x, y) planes. θ₁ is the angle between x and axes, θ₂ is the angle between y and axes, and θ₃ is the angle between x and axes. A schematic of the coordinate transform in Eq. 8a is shown in Figure 2.

FIGURE 2

By using Eqs. 8a–c, the mixed spatial derivatives in Eq. 7f can be transformed asFor simplicity, we usually use equal grid sizes in the trapezoid coordinate system (Δx = Δy = Δz = Δ), which means . By substituting Eq. 7 and Eq. 9 into Eq. 6, we get the time-domain-discretization form of the 3D acoustic equation with the CPML absorbing boundary condition in the trapezoid coordinate system aswhere u^j = u(x, y, z, t_j) represents the scalar wavefield at the jth time step in the trapezoid coordinate system; (x^s, y^s, z^s) is the position of the source in the trapezoid coordinate system; ln, where L_τ is the thickness of CPML absorbing boundary along the τ direction, denotes the distance to the inner area in the τ direction; and , τ ∈ {x, y, z}.

A schematic of the grid discretization in the 3D trapezoid-grid CPML area is shown in Figure 3. In this work, 30 and 20 absorbing boundary layers are usually used for the trapezoid-grid CPML area in the horizontal and vertical directions, respectively.

FIGURE 3

2.3 Stability Analysis

Stability condition is usually required for the FD scheme to give a stable time step. From Eq. 10a, we use a local frozen coefficients technique in each discrete point and can get the full-discretization form of the 3D trapezoid coordinate system acoustic equation without the CPML boundary condition and source function:where is the wavefield at (x_m, y_n, z_l, t_j), x_m = x⁰ + (m − 1)Δx, y_n = y⁰ + (n − 1)Δy, z_l = z⁰ + (l − 1)Δz, t_j = t⁰ + (j − 1)Δt, N_x, N_y, N_z, N_xy, N_xz, N_yz are half-of-spatial FD orders, η^x, η^y, η^z, η^xy, η^xz, η^yz are corresponding FD coefficients of the second-order derivative, and c^x, c^y, c^z, c^xy, c^xz, c^yz are corresponding FD coefficients of the first-order derivative.

To derive the stability condition, we use the plane wave solution that is defined aswhere is the amplitude of the plane wave, i is the imaginary unit, ω is the angular frequency, and k_x, k_y, k_z are wavenumbers in the x-, y- and z-directions, respectively. Similar to stability analysis of Kosloff and Baysal (1982), by substituting Eq. 12 into Eq. 11 and only considering the maximum wavenumber, the stability condition of the 3D acoustic equation in the trapezoid coordinate system can be expressed aswhere mod is the function for the getting remainder, and represents the maximum value of the objective function at those discrete points (x_m, y_n, z_l).

3 Numerical Results

In the following numerical examples, Eq. 10 is discretized by the eighth-order FD in the space, and conventional Taylor-expansion–based high-order FD coefficients (Dablain, 1986) are adopted.

3.1 Homogenous Model

First, we use a 3D homogenous model with a constant velocity of 2000 m/s to verify the effectiveness of our trapezoid-grid FDTD method and corresponding CPML absorbing boundary condition. A Ricker wavelet with a dominant frequency of 20 Hz is located at the center of the model as the source. The FD time step is taken as 1.6 ms. The scaling parameter γ is set as 2.78 × 10^–4, the sampling function g(z) = z, and the lateral grid sizes in the Cartesian coordinate system increase from 7.5 to 10 from top to bottom. Figure 4A shows the snapshot obtained by the trapezoid-grid FDTD with CPML at 0.45 s, while Figure 4B shows the corresponding snapshot without CPML. Figure 5 shows the comparison between the recorded seismograms computed by the uniform-grid FDTD and the trapezoid-grid FDTD at (x₀, y₀, z₀) = (600, 600, 0 m). The comparison between Figures 4A,B demonstrates that trapezoid-grid CPML can effectively reduce boundary reflections, while Figure 5 demonstrates the accuracy of the trapezoid-grid FDTD method for the homogenous model.

FIGURE 4

FIGURE 5

3.2 Overthrust Model

Then, we apply our method to the SEG/EAGE overthrust model (Figure 6A), which is based on the real overthrusts of South America. Figures 6B,C show the modeling area of the SEG/EAGE overthrust model in the trapezoid coordinate system and the Cartesian coordinate system, respectively. A Ricker wavelet with the dominate frequency of 4.2 Hz is located at (10 km, 10 km, 0.5 km) as the source. The grid sizes for the uniform-grid FDTD are 50 m × 50 m × 50 m, which means the minimum NPPW in each direction is close to 10. We therefore set the minimum NPPW in x₀-, y₀-, and z₀-direction as 10 for the trapezoid-grid FDTD method, and get the scaling parameter as γ = 2.07 × 10^–4. Figure 7A shows the vertical sampling function g(z) used for this model, and the vertical grid sizes in the Cartesian coordinate system increase from 51.8 m in the shallow region to a maximum value of 142.6 m in the deep region, as shown in Figure 7B. Based on the stability analysis, the time step for the trapezoid-grid FDTD and the uniform-grid FDTD are 3.697 and 3.585 ms, respectively. Receivers are located on the surface along y₀ = 10 km. Figure 8A shows the snapshot at 2.5 s computed by the uniform-grid FDTD, and Figure 8B is the snapshot at 2.5 s in the trapezoid coordinate system computed by our trapezoid-grid method. Using coordinate transformation and cubic spline interpolation, we can get the corresponding snapshot in the Cartesian coordinate system, as shown in Figure 8C. Figures 8A,C show good agreement. To give more detailed comparisons, single-trace seismograms at (7.5, 10, 0 km), (10, 10, 0 km), and (12.5, 10, 0 km) for both the uniform-grid (black solid line) and the trapezoid-grid (red dash line) FDTD are shown in Figure 9. Figure 9 also shows good agreement between the uniform-gird FDTD and the trapezoid-grid FDTD. On our computing platform (Intel(R) Xeon(R) Sliver 4216 CPU @ 2.10GHz, 256GB of memory, and C++ codes), using 16-threads computation and similar code optimization techniques, the running time for the trapezoid-grid FDTD and the uniform-grid FDTD is calculated as 2203 s and 2925 s, respectively, which shows a calculation efficiency improvement of 24.7%. The memories for the trapezoid-grid FDTD and the uniform-grid FDTD are about 336 and 1213 MB, respectively, which shows a memory reduction of 72.3%. Considering that the simulation area of our trapezoid-grid method is almost 60% of that of the uniform-grid method, for the common simulation area, we can achieve about 50% reduction on memory usage.

FIGURE 6

FIGURE 7

FIGURE 8

FIGURE 9

4 Conclusion

In this article, we propose a 3D trapezoid-grid FDTD seismic wave modeling method based on the increasing trend of seismic wave velocity with depth. The trapezoid-grid mesh in the physical Cartesian system can effectively reduce the oversampling in the high-velocity region compared with the uniform-grid method, and the design of 3D trapezoid coordinate transform greatly avoids the difficulty of processing an irregular grid. We derive the 3D acoustic equation in the trapezoid coordinate system. The corresponding CPML boundary condition is also given to decrease artificial boundary reflection. To obtain a stable and efficient wave modeling result, we combine the plane wave theory and frozen coefficients technique and provide an effective stability condition for the 3D trapezoid-grid FDTD method. The discretization of the 3D acoustic equation in the trapezoid coordinate system is completed by the eighth order and second order finite-difference method in the space and time domain, respectively. The 3D homogenous model is given to verify the effectiveness of trapezoid-grid FDTD and the performance of the CPML boundary. Numerical tests on the SEG/EAGE overthrust model indicate the accuracy and the significant memory reduction of our method compared with uniform-grid FDTD. The key idea of our method is the combination of the trapezoid coordinate transformation and the FD stencils. Such idea can be generalized to many other wave equations such as elastic equation (Zhan et al., 2017) and Maxwell’s equations (Zhan et al., 2021). Besides, our method is actually dealing with the regular grids in the trapezoid coordinate system, which means that we can combine other methods to treat the irregular surface (Li et al., 2020) or curved interfaces (Zhan et al., 2020).

References

• 1

  AbreuR.StichD.MoralesJ. (2015). The Complex-Step-Finite-Difference Method. Geophys. J. Int.202, 72–93. 10.1093/gji/ggv125

  • CrossRef
  • Google Scholar

• 2

  AntunesA. J. M.Leal-ToledoR. C. P.FilhoO. T. d. S.ToledoE. M. (2014). Finite Difference Method for Solving Acoustic Wave Equation Using Locally Adjustable Time-Steps. Proced. Computer Sci.29, 627–636. 10.1016/j.procs.2014.05.056

  • CrossRef
  • Google Scholar

• 3

  BaysalE.KosloffD. D.SherwoodJ. W. C. (1983). Reverse Time Migration. Geophysics48, 1514–1524. 10.1190/1.1441434

  • CrossRef
  • Google Scholar

• 4

  CaiJ.ZhangJ. (2015). “Acoustic Full Waveform Inversion with Physical Model Data,” in 2015 Workshop: Depth Model Building: Full-Waveform Inversion. Beijing, China: Society of Exploration Geophysicists (SEG), 146–149. 10.1190/FWI2015-036

  • CrossRef
  • Google Scholar

• 5

  ChenF.XuS. (2012). “Pyramid-shaped Grid for Elastic Wave Propagation,” in SEG Technical Program Expanded Abstracts 2012 (Las Vegas: Society of Exploration Geophysicists (SEG)), 1–5. 10.1190/segam2012-0890.1

  • CrossRef
  • Google Scholar

• 6

  DablainM. A. (1986). The Application of High‐order Differencing to the Scalar Wave Equation. GEOPHYSICS51, 54–66. 10.1190/1.1442040

  • CrossRef
  • Google Scholar

• 7

  DuQ.ZhangX.ZhangS.ZhangF.FuL.-Y. (2021). The Pseudo-laplace Filter for Vector-Based Elastic Reverse Time Migration. Front. Earth Sci.9, 538. 10.3389/feart.2021.687835

  • CrossRef
  • Google Scholar

• 8

  FujiiY.AzumiT.NishioN.KatoS.EdahiroM. (2013). “Data Transfer Matters for Gpu Computing” in Proceeding of the 2013 International Conference on Parallel and Distributed Systems, Seoul, Korea (South), 15-18 Dec. 2013, (IEEE), 275–282. 10.1109/ICPADS.2013.47

  • CrossRef
  • Google Scholar

• 9

  GaoJ.XuW.WuB.LiB.ZhaoH. (2018). Trapezoid Grid Finite Difference Seismic Wavefield Simulation with Uniform Depth Sampling Interval. Chin. J. Geophys.61, 3285–3296. 10.6038/cjg2018M0313

  • CrossRef
  • Google Scholar

• 10

  HayashiK.BurnsD. R. (1999). “Variable Grid Finite‐difference Modeling Including Surface Topography,” in SEG Technical Program Expanded Abstracts 1999 (Houston: Society of Exploration Geophysicists (SEG)), 528–531. 10.1190/1.1821071

  • CrossRef
  • Google Scholar

• 11

  Hong-QiaoX.Xiao-YiW.Chen-YuanW.Jiang-JieZ. (2021a). Sparse Constrained Least-Squares Reverse Time Migration Based on Kirchhoff Approximation. Front. Earth Sci.9, 770. 10.3389/feart.2021.731697

  • CrossRef
  • Google Scholar

• 12

  HuangC.DongL.-G. (2009). Staggered-Grid High-Order Finite-Difference Method in Elastic Wave Simulation with Variable Grids and Local Time-Steps. Chin. J. Geophys.52, 1324–1333. 10.1002/cjg2.1457

  • CrossRef
  • Google Scholar

• 13

  JiaJ.WuB.PengJ.GaoJ. (2019). Recursive Linearization Method for Inverse Medium Scattering Problems with Complex Mixture Gaussian Error Learning. Inverse Probl.35, 075003. 10.1088/1361-6420/ab08f2

  • CrossRef
  • Google Scholar

• 14

  JingH.YangG.WangJ. (2019). An Optimized Time-Space-Domain Finite Difference Method with Piecewise Constant Interpolation Coefficients for Scalar Wave Propagation. J. Geophys. Eng.16, 309–324. 10.1093/jge/gxz008

  • CrossRef
  • Google Scholar

• 15

  KosloffD. D.BaysalE. (1982). Forward Modeling by a Fourier Method. GEOPHYSICS47, 1402–1412. 10.1190/1.1441288

  • CrossRef
  • Google Scholar

• 16

  LiJ.TsengH.-W.LinC.PapakonstantinouY.SwansonS. (2016). HippogriffDB. Proc. VLDB Endow.9, 1647–1658. 10.14778/3007328.3007331

  • CrossRef
  • Google Scholar

• 17

  LiX.YaoG.NiuF.WuD. (2020). An Immersed Boundary Method with Iterative Symmetric Interpolation for Irregular Surface Topography in Seismic Wavefield Modelling. J. Geophys. Eng.17, 643–660. 10.1093/jge/gxaa019

  • CrossRef
  • Google Scholar

• 18

  LiuS.YanZ.ZhuW.HanB.GuH.HuS. (2021). An Illumination‐compensated Gaussian Beam Migration for Enhancing Subsalt Imaging. Geophys. Prospecting69, 1433–1440. 10.1111/1365-2478.13117

  • CrossRef
  • Google Scholar

• 19

  LiuX.YinX.WuG. (2014). Finite-difference Modeling with Variable Grid-Size and Adaptive Time-step in Porous media. Earthq Sci.27, 169–178. 10.1007/s11589-013-0055-7

  • CrossRef
  • Google Scholar

• 20

  LiuY. (2020). Acoustic and Elastic Finite-Difference Modeling by Optimal Variable-Length Spatial Operators. GEOPHYSICS85, T57–T70. 10.1190/geo2019-0145.1

  • CrossRef
  • Google Scholar

• 21

  LiuY.SenM. K. (2011a). Finite-difference Modeling with Adaptive Variable-Length Spatial Operators. GEOPHYSICS76, T79–T89. 10.1190/1.3587223

  • CrossRef
  • Google Scholar

• 22

  LiuY.SenM. K. (2011b). Scalar Wave Equation Modeling with Time-Space Domain Dispersion-Relation-Based Staggered-Grid Finite-Difference Schemes. Bull. Seismological Soc. America101, 141–159. 10.1785/0120100041

  • CrossRef
  • Google Scholar

• 23

  MaX.YangD.HeX.HuangX.SongJ. (2019). Nonsplit Complex-Frequency-Shifted Perfectly Matched Layer Combined with Symplectic Methods for Solving Second-Order Seismic Wave Equations - Part 2: Wavefield Simulations. GEOPHYSICS84, T167–T179. 10.1190/geo2018-0349.1

  • CrossRef
  • Google Scholar

• 24

  MaX.YangD.HuangX.ZhouY. (2018). Nonsplit Complex-Frequency Shifted Perfectly Matched Layer Combined with Symplectic Methods for Solving Second-Order Seismic Wave Equations - Part 1: Method. GEOPHYSICS83, T301–T311. 10.1190/geo2017-0603.1

  • CrossRef
  • Google Scholar

• 25

  PasalicD.McGarryR. (2010). “A Discontinuous Mesh Finite Difference Scheme for Acoustic Wave Equations,” in SEG Technical Program Expanded Abstracts 2010 (Denver: Society of Exploration Geophysicists (SEG)), 2940–2944. 10.1190/1.3513457

  • CrossRef
  • Google Scholar

• 26

  QuY.LiZ.HuangJ.DengW.LiJ. (2015). “The Least-Squares Reverse Time Migration for Viscoacoustic Medium Based on a Stable Reverse-Time Propagator,” in SEG Technical Program Expanded Abstracts 2015(New Orleans: Society of Exploration Geophysicists (SEG)), 3977–3980. 10.1190/segam2015-5835196.1

  • CrossRef
  • Google Scholar

• 27

  RobertssonJ. O. A.BlanchJ. O. (2020). “Numerical Methods, Finite Difference,” in Encyclopedia of Solid Earth Geophysics. Editor GuptaH. K. (Cham: Springer International Publishing), 1–9. 10.1007/978-3-030-10475-7_135-1

  • CrossRef
  • Google Scholar

• 28

  TarantolaA. (1984). Inversion of Seismic Reflection Data in the Acoustic Approximation. Geophysics49, 1259–1266. 10.1190/1.1441754

  • CrossRef
  • Google Scholar

• 29

  VirieuxJ.OpertoS. (2009). An Overview of Full-Waveform Inversion in Exploration Geophysics. Geophysics74, WCC1–WCC26. 10.1190/1.3238367

  • CrossRef
  • Google Scholar

• 30

  WuB.LiB.YangH.JiaJ. (2019a). “Trapezoid Grid Finite Difference for Acoustic Wave Modeling,” in SEG 2018 Workshop: SEG Seismic Imaging Workshop. Beijing, China: Society of Exploration Geophysicists (SEG), 52–55. 12–14 November 2018. 10.1190/SEIM2018-13.1

  • CrossRef
  • Google Scholar

• 31

  WuB.XuW.JiaJ.LiB.YangH.ZhaoH.et al (2018). “Convolutional Perfect-Matched Layer Boundary for Trapezoid Grid Finite-Difference Seismic Modeling,” in SEG Technical Program Expanded Abstracts 2018 (Anaheim: Society of Exploration Geophysicists (SEG)), 3989–3993. 10.1190/segam2018-2995754.1

  • CrossRef
  • Google Scholar

• 32

  WuB.XuW.LiB.JiaJ. (2019b). Trapezoid Coordinate Finite Difference Modeling of Acoustic Wave Propagation Using the CPML Boundary Condition. J. Appl. Geophys.168, 101–106. 10.1016/j.jappgeo.2019.06.006

  • CrossRef
  • Google Scholar

• 33

  XiaD.-m.XieC.SongP.TanJ.LiJ.-s.ZhaoB. (2017). “A Time Domain Full Waveform Inversion Method Based on Well-Constrained Regularization,” in SEG 2017 Workshop: Full-Waveform Inversion and beyond. Beijing, China: Society of Exploration Geophysicists (SEG), 136–139. 20-22 November 2017. 10.1190/FWI2017-036

  • CrossRef
  • Google Scholar

• 34

  XuW.GaoJ. (2018). Adaptive 9-point Frequency-Domain Finite Difference Scheme for Wavefield Modeling of 2D Acoustic Wave Equation. J. Geophys. Eng.15, 1432–1445. 10.1088/1742-2140/aab015

  • CrossRef
  • Google Scholar

• 35

  XuW.ZhongY.WuB.GaoJ.LiuQ. H. (2021b). Adaptive Complex Frequency with V-Cycle Gmres for Preconditioning 3D Helmholtz Equation. GEOPHYSICS86, T349–T359. 10.1190/geo2020-0901.1

  • CrossRef
  • Google Scholar

• 36

  XuanK.YingS.XuebaoG.ShizhuL. (2014). “Carbonate Reservoir Pre-stack Reverse-Time Migration and Gpu/cpu Heterogeneous Computing Research,” in Beijing 2014 International Geophysical Conference & Exposition. Beijing, China: Society of Exploration Geophysicists (SEG), 21–24. April 2014. 1101–1104. 10.1190/IGCBeijing2014-279

  • CrossRef
  • Google Scholar

• 37

  ZhanQ.SunQ.RenQ.FangY.WangH.LiuQ. H. (2017). A Discontinuous Galerkin Method for Simulating the Effects of Arbitrary Discrete Fractures on Elastic Wave Propagation. Geophys. J. Int.210, 1219–1230. 10.1093/gji/ggx233

  • CrossRef
  • Google Scholar

• 38

  ZhanQ.WangY.FangY.RenQ.YangS.YinW.-Y.et al (2021). An Adaptive High-Order Transient Algorithm to Solve Large-Scale Anisotropic Maxwell's Equations. IEEE Trans. Antennas Propagat., 1. 10.1109/tap.2021.3111639

  • CrossRef
  • Google Scholar

• 39

  ZhanQ.ZhuangM.MaoY.LiuQ. H. (2020). Unified Riemann Solution for Multi-Physics Coupling: Anisotropic Poroelastic/elastic/fluid Interfaces. J. Comput. Phys.402, 108961. 10.1016/j.jcp.2019.108961

  • CrossRef
  • Google Scholar

• 40

  ZhouH.LiuY.WangJ. (2021). Optimizing Orthogonal-Octahedron Finite-Difference Scheme for 3D Acoustic Wave Modeling by Combination of taylor-series Expansion and Remez Exchange Method. Exploration Geophys.52, 335–355. 10.1080/08123985.2020.1826890

  • CrossRef
  • Google Scholar