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 as

where (x0, y0, z0) 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 z0-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:

gz+Δz=gz+vmingzf0N0,(4b)

where f0 is the dominant frequency of the source term, N0 is the preferred NPPW and is related to the accuracy of the adopted FD scheme, and vmin(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 vmin(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 vmin(g(z)). In particular, we usually set the order of the local polynomial as three. Such vmin(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 z0-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 as

In our work, γ and g(z) are determined adaptively according to the model velocity. By selecting γ such that Δx0(z) and Δy0(z) are always smaller than or equal to Δz0(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.

www.frontiersin.org

FIGURE 1. Schematic of the 3D trapezoid coordinate transformation: (A) the trapezoid-grid mesh in the Cartesian coordinate system; (B) the uniform grid mesh in the transformed trapezoid coordinate system. The two gray regions in A and B represent the same physical region.

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 as

1v2uj+12uj+uj1Δt22ujx022ujy022ujz02ψx0jx0ψy0jy0ψz0jz0ζx0jζy0jζz0j=ftjδx0x0sδy0y0sδz0z0s;(6a)
ψτ0j+1=aτ0ψτ0j+bτ0uj+1τ0;(6b)
ζτ0j+1=aτ0ζτ0j+bτ02uj+1τ02+ψτ0j+1τ0;(6c)
aτ0=eστ0+ατ0Δt;(6d)
bτ0=στ0στ0+ατ0aτ01,τ0x0,y0,z0,(6e)

where uj = u(x0, y0, z0, tj) 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;

(x0s,y0s,z0s)

is the position of source;

στ0=3vmax2Lτ0τ̄0Lτ02

ln

1R

, where R denotes the designated theoretical boundary reflection coefficient, vmax is the maximum velocity of the model,

Lτ0

is the thickness of CPML absorbing boundary along the τ0 direction, and

τ̄0

denotes the distance to the inner area in the τ0 direction;

ατ0=αmax1τ̄0Lτ0

and αmax = πf0.

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

x0=11+γgzx;(7a)
y0=11+γgzy;(7b)
z0=γx1+γgzxγy1+γgzy+1gzz;(7c)
2x02=11+γgz22x2;(7d)
2y02=11+γgz22y2;(7e)
2z02=2γ2x1+γgz2x+γ2x21+γgz22x2+2γ2xy1+γgz22xy2γx1+γgzgz2xz+2γ2y1+γgz2y+γ2y21+γgz22y22γy1+γgzgz2yzgzgz3z+1gz22z2.(7f)

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 as

xz=cosθ1sinθ1sinθ1cosθ1x̂ẑ;(8a)
yz=cosθ2sinθ2sinθ2cosθ2ỹz̃;(8b)
xy=cosθ3sinθ3sinθ3cosθ3x̄ȳ,(8c)

where

x̂

and

ẑ

are axes along diagonals in the (x, z) planes,

ỹ

and

z̃

are axes along diagonals in the (y, z) planes, and

x̄

and

ȳ

are axes along diagonals in the (x, y) planes. θ1 is the angle between x and

x̂

axes, θ2 is the angle between y and

ỹ

axes, and θ3 is the angle between x and

x̄

axes. A schematic of the coordinate transform in Eq. 8a is shown in Figure 2.

www.frontiersin.org

FIGURE 2. Schematic of the rotation transformation in the trapezoid coordinate system for transforming mixed spatial derivatives into non-mixed spatial derivatives.

By using Eqs. 8a–c, the mixed spatial derivatives in Eq. 7f can be transformed as

2xz=12sin2θ12x̂22ẑ2;(9a)
2yz=12sin2θ22ỹ22z̃2;(9b)
2xy=12sin2θ32x̄22ȳ2.(9c)

For simplicity, we usually use equal grid sizes in the trapezoid coordinate system (Δx = Δy = Δz = Δ), which means

θ1=θ2=θ3=π4

. 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 as

1v2uj+12uj+uj1Δt2γ2x2+11+γgz22ujx2γ2y2+11+γgz22ujy22γ2x1+γgz2ujx2γ2y1+γgz2ujy+γx1+γgzgz2ujx̂22ujẑ2+γy1+γgzgz2ujỹ22ujz̃2γ2xy1+γgz22ujx̄22ujȳ2+gzgz3ujz1gz22ujz211+γgzψxjx11+γgzψyjy+γx1+γgzψzjz+γy1+γgzψzjy1gzψzjzζxjζyjζzj=ftjxxsyyszzs;(10a)
ψxj+1=axψxj+bx11+γgzuj+1x;(10b)
ψyj+1=ayψyj+by11+γgzuj+1y;(10c)
ψzj+1=azψzj+bzγx1+γgzuj+1xγy1+γgzuj+1y+1gzuj+1z;(10d)
ζxj+1=axζxj+bx11+γgz22uj+1x2+11+γgzψxj+1x;(10e)
ζyj+1=ayζyj+by11+γgz22uj+1y2+11+γgzψyj+1y;(10f)
ζzj+1=azζzj+bzγ2x21+γgz22uj+1x2+γ2y21+γgz22uj+1y2+2γ2x1+γgz2uj+1x+2γ2y1+γgz2uj+1yγx1+γgzgz2uj+1x̂22uj+1ẑ2γy1+γgzgz2uj+1ỹ22uj+1z̃2+γ2xy1+γgz22uj+1x̄22uj+1ȳ2gzgz3uj+1z+1gz22uj+1z2γx1+γgzψzj+1xγy1+γgzψzj+1y+1gzψzj+1z],(10g)

where uj = u(x, y, z, tj) represents the scalar wavefield at the jth time step in the trapezoid coordinate system; (xs, ys, zs) is the position of the source in the trapezoid coordinate system;

στ=3vmax2Lττ̄Lτ2

ln

1R

, where Lτ is the thickness of CPML absorbing boundary along the τ direction,

τ̄

denotes the distance to the inner area in the τ direction; and

ατ=αmax1τ̄Lτ

, τ ∈ {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.

www.frontiersin.org

FIGURE 3. Schematic of the grid discretization in the 3D trapezoid-grid CPML area.

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:

1v2um,n,lj+12um,n,lj+um,n,lj1Δt2γ2xm2+11+γgzl21Δ2p=1Nxηpxump,n,lj2um,n,lj+um+p,n,ljγ2yn2+11+γgzl21Δ2p=1Nyηpyum,np,lj2um,n,lj+um,n+p,lj2γ2xm1+γgzl21Δp=1Nxcpxum+p,n,ljump,n,lj2γ2yn1+γgzl21Δp=1Nycpyum,n+p,ljum,np,lj+γxm1+γgzlgzl12Δ2p=1Nxzηpxzump,n,lpj+um+p,n,l+pjump,n,l+pjum+p,n,lpj+γyn1+γgzlgzl12Δ2p=1Nyzηpyzum,np,lpj+um,n+p,l+pjum,np,l+pjum,n+p,lpjγ2xmyn1+γgzl212Δ2p=1Nxyηpxyump,np,lj+um+p,n+p,ljump,n+p,ljum+p,np,lj+gzlgzl31Δp=1Nzcpzum,n,l+pjum,n,lpj1gzl21Δ2p=1Nzηpzum,n,lpj2um,n,lj+um,n,l+pj=0,(11)

where

um,n,lj

is the wavefield at (xm, yn, zl, tj), xm = x0 + (m − 1)Δx, yn = y0 + (n − 1)Δy, zl = z0 + (l − 1)Δz, tj = t0 + (j − 1)Δt, Nx, Ny, Nz, Nxy, Nxz, Nyz are half-of-spatial FD orders, ηx, ηy, ηz, ηxy, ηxz, ηyz are corresponding FD coefficients of the second-order derivative, and cx, cy, cz, cxy, cxz, cyz are corresponding FD coefficients of the first-order derivative.

To derive the stability condition, we use the plane wave solution that is defined as

ux,y,z,t=u0*eiωtikxxikyyikzz,(12)

where

u0*

is the amplitude of the plane wave, i is the imaginary unit, ω is the angular frequency, and kx, ky, kz 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 as

Δt<Δvmaxmaxm,n,lγ2xm2+11+γgzl2p=1Nxmodp,2ηpx+γ2yn2+11+γgzl2p=1Nymodp,2ηpy+1gzl2p=1Nzmodp,2ηpz12,(13)

where mod is the function for the getting remainder, and

maxm,n,l

represents the maximum value of the objective function at those discrete points (xm, yn, zl).

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 (x0, y0, z0) = (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.

www.frontiersin.org

FIGURE 4. Snapshots obtained by the trapezoid-grid FDTD for the homogenous model (A) with CPML and (B) without CPML.

www.frontiersin.org

FIGURE 5. Comparison of seismograms obtained by the uniform-grid FDTD and the trapezoid-grid FDTD for the homogenous model at (x0, y0, z0) = (600 m, 600 m, 0 m). The black solid line represents the uniform-grid FDTD result, and the red dash line represents the trapezoid-grid FDTD result.

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 x0-, y0-, and z0-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 y0 = 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.

www.frontiersin.org

FIGURE 6. (A) SEG/EAGE overthrust model; (B) SEG/EAGE overthrust model in the trapezoid coordinate system; (C) actual simulating area for the trapezoid-grid FDTD in the Cartesian coordinate system.

www.frontiersin.org

FIGURE 7. (A) Vertical sampling function g(z); (B) variation of Δz0 for the trapezoid-grid FDTD.

www.frontiersin.org

FIGURE 8. Snapshots for the SEG/EAGE overthrust model at 2.5 s: (A) uniform-grid FDTD result; (B) trapezoid-grid FDTD result in the trapezoid coordinate system; (C) trapezoid-grid FDTD result in the Cartesian coordinate system.

www.frontiersin.org

FIGURE 9. Single-trace seismograms for the SEG/EAGE overthrust model at the receiver position of (A) (7.5, 10, 0 km); (B) (10, 10, 0 km); and (C) (12.5, 10, 0 km). The black solid lines represent the uniform-grid FDTD results, and the red dash lines represent the trapezoid-grid FDTD results.

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).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author Contributions

BW contributed to the conception of the study. Writing, compiling, and debugging of the program were completed by WX. WT wrote the first draft of the manuscript, and all authors contributed to the revision.

Funding

We would like to thank the Natural Science Foundation of Shaanxi Province under Grant No. 2020JM-018 and National Natural Science Foundation (41974122) of China for funding this work.

Conflict of Interest

Authors BW, WT, and BL are employed by SINOPEC.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References


Disclaimer:

This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (https://www.frontiersin.org/)