Zhenzhu Fang
Zhundong Second Mine of State Grid Energy Xinjiang Zhundong Coal and Electricity Co., Ltd., Changji, China.
DOI: 10.4236/ojg.2025.159027   PDF    HTML   XML   78 Downloads   421 Views   Citations

Abstract

To overcome the problem of calculation errors in the Born approximation when the forward accumulation effect is strong in VTI media, this article combines the De Wolf approximation method with VTI media to propose a numerical simulation method for scattering waves in VTI media based on the De Wolf approximation. Scattering equations are divided into forward and backward scattering equations, and the Born series is renormalized to improve its convergence in VTI media. During the calculation process, multiple interlayer waves are ignored, and only multiple forward and one backward scattering signals are retained. To accelerate computational efficiency, thin slab approximation and screen approximation methods are employed for rapid implementation. I use the concave model and complex model for algorithm testing, and the research results demonstrate that the De Wolf approximation can effectively mitigate the issue of calculation errors in the Born approximation, as well as reveal the propagation law of seismic waves in VTI media. Due to the influence of the anisotropic parameters of the medium, there may be phase deviations in the forward modeling records of seismic waves using the thin slab approximation method and the screen approximation method; because the thin slab approximation method and the screen approximation method only account for one reflection, there are no signals of multiple waves in the seismic records, and the comparative results verify the effectiveness of the algorithm.

Keywords

De Wolf Approximation, VTI Medium, Lippmann-Schwinger Equation, Scattering Wave Numerical Simulation, Screen Approximation Method

Share and Cite:

1. Introduction

VTI medium (transversely isotropic media with a vertical axis of symmetry) is an approximate representation of the anisotropic characteristics of actual Earth media (underground sedimentary layers are mostly layered and exhibit lateral isotropy), which can affect the travel time information of seismic waves [1] [2]. When the geological structure of the exploration target is complex, faults are developed, the dip angle of the strata is steep, or the lithology changes laterally, and heterogeneous geological bodies of different scales coexist, an extremely complex seismic wavefield with multiple wave groups interfering with each other will be formed [3]. In this case, using only reflected and diffracted wave information cannot achieve precise imaging of complex areas. Research has shown that scattered waves also carry geometric and physical information related to complex structures and complex rock types [4]. Therefore, the numerical simulation of scattering waves in VTI media has research value.

The numerical simulation methods of scattered wave fields can be divided into two categories: deterministic methods and statistical methods. The former is a numerical simulation of special velocity structures such as terraces and sedimentary basins through numerical calculations, while the latter is based on studying the statistical characteristics of the medium, treating velocity or terrain as random variables rather than studying the determined position of each scatter in the medium. At present, the widely used numerical simulation methods mainly include the finite difference method [5] [6], finite element method [7], volume integral equation method [8], F-K domain integration method [9], and other methods. The finite difference method can analyze the distribution of scattering fields and the characteristics of wave field transport without any limitation on the scale of velocity changes. However, the finite difference method also has unavoidable limitations. The simulated wavefield contains non-scattering information, and the scattered wave energy is very small. Its detailed information is submerged by a single location [10]. In addition, this method is also limited by computer memory and computational accuracy. Compared with the finite difference method, the F-K domain integration method does not include a first-order field in the simulated wavefield, and common direct waves do not appear in the profile. Moreover, the accuracy of this algorithm is higher than that of the finite difference method. However, its computational efficiency is low, and when the number of calculation points is the same, its computational workload is much larger than that of the finite difference method [11]. The advantage of the finite element method is that it can simulate infinite complex bodies with finite and interrelated elements. No matter how complex the geometry is, it can be simplified with corresponding elements to model, analyze, and calculate the results [12]. Simplifying complex engineering problems that seem difficult to approach is the greatest advantage. The finite element method is expressed in matrix form and has high programmability. For linear elastic problems, convergent solutions can be obtained when the actual structural displacement field function is continuous and smooth [13]. In theory, it is always possible to obtain sufficiently approximate simulations for any complex structure through the method of subdividing units. However, this method has a huge computational load and is difficult to apply in numerical simulations of scattered waves in large models. The L-S equation belongs to the body integral equation, which can effectively describe the scattering process of seismic waves in non-uniform media (primary scattering and multiple scattering). Its equation has semi-analytical characteristics and is a powerful tool for studying the propagation of scattered waves [14]-[16].

If we use the iterative series method to obtain the Born series, the convergence speed of the Born series is slow in strongly disturbed media. The Born approximation is a first-order approximation. However, the Born approximation is a weak scattering approximation and is not suitable for long-distance propagation at high frequencies. To solve this problem, some experts introduced the concept of renormalized scattering series, which improves the convergence of the Born series in strongly perturbed media by renormalizing it, i.e., obtaining more accurate numerical solutions with fewer iterations. Regarding this issue, Mr. Wu et al. developed the De Wolf approximation (the MFSB approximation or local Born approximation), among which the thin slab approximation and the screen approximation are two typical implementations of the De Wolf approximation [17] [18]. Compared to the full-wave finite difference and finite element methods, the advantage of this unidirectional propagation method is its fast calculation speed, which can save a lot of memory during seismic wave forward simulation [19]. This method is widely applied to the propagation of seismic wave scattering wavefields in non-uniform media [20]-[26]. This method is accurate and efficient for modeling seismic waves, especially for long-distance propagation of high-frequency seismic waves. In the 1990s, Wu used this method to simulate the propagation of high-frequency elastic Lg waves excited by nuclear tests, which are difficult to achieve using traditional full-wave finite difference methods [27]. Sun and Sun [2] applied the De Wolf approximation to one-way wave migration imaging in VTI media. However, the numerical simulation application of the De Wolf approximation for scattered waves in VTI media has not been reported yet.

Based on this, this article combines the De Wolf approximation with VTI media to derive a numerical simulation method for scattering waves in VTI media based on the De Wolf approximation. The scattering equations are divided into forward and backward scattering equations, and the Born series is renormalized to improve its convergence in VTI media. During the calculation process, multiple interlayer waves are ignored, and only multiple forward and one backward scattering signals are retained. To accelerate computational efficiency, thin slab approximation and screen approximation methods are employed for rapid implementation. This article uses the concave model and complex model for algorithm testing to further verify the effectiveness of the algorithm.

2. Principle

2.1. Lippmann-Schwinger Equation for VTI Media

The dispersion equation under two-dimensional VTI medium conditions:

${\omega }^4-\left(1+2\epsilon \right)v_{qp}^2{\omega }^2{k}_x^2-v_{qp}^2{\omega }^2{k}_z^2+2\left(\epsilon -\delta \right)v_{qp}^4{k}_x^2{k}_z^2=0$ . (1)

where $\omega$ is the circular frequency, $v_{qp}$ is the velocity parameter, $\epsilon$ and $\delta$ are the anisotropy parameters, and $k_x$ and $k_z$ are divided into transverse and vertical wavenumbers. After simplification, I obtain:

$k_x^2+k_z^2=\frac{{\omega }^2/v_{qp}^2-\chi k_x^2+\eta k_z^2}{\eta }$ , (2)

where $\chi =1+2\epsilon$ , $\eta =1+2\left(\epsilon -\delta \right)\frac{v_{qp}^2}{{\omega }^2}{k}_x^2$ .

Transform Formula (2) into the frequency space domain and add spatial coordinates:

$\begin{array}{l}\frac{{\partial }^2}{\partial x^2}p\left(x,z;\omega \right)+\frac{{\partial }^2}{\partial z^2}p\left(x,z;\omega \right)\\ =-\frac{\left[\frac{{\omega }^2}{v_{qp}^2\left(x,z\right)}+\chi \left(x,z\right)\frac{{\partial }^2}{\partial x^2}+\eta \left(x,z\right)\frac{{\partial }^2}{\partial x^2}\right]p\left(x,z;\omega \right)}{\eta \left(x,z\right)},\end{array}$ (3)

where $\chi \left(x,z\right)=1+2\epsilon \left(x,z\right)$ , $\eta \left(x,z\right)=1+2\left[\epsilon \left(x,z\right)-\delta \left(x,z\right)\right]\frac{v_{qp}^2\left(x,z\right)}{{\omega }^2}\frac{{\partial }^2}{\partial x^2}$ , $p\left(x,z;\omega \right)$ is the frequency space domain wavefield.

I set $v_0\left(z\right)$ as the background velocity, ${\epsilon }_0\left(z\right)$ and ${\delta }_0\left(z\right)$ as background anisotropy parameters, which are only functions of depth $z$ and are constant within each layer. Let:

$F_0\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)=\frac{\left[\frac{{\omega }^2}{v_0^2\left(z\right)}+{\chi }_0\left(z\right)\frac{{\partial }^2}{\partial x^2}+{\eta }_0\left(x,z\right)\frac{{\partial }^2}{\partial x^2}\right]}{{\eta }_0\left(x,z\right)}$ , (4)

where ${\eta }_0\left(x,z\right)=1+2\left[{\epsilon }_0\left(x,z\right)-{\delta }_0\left(x,z\right)\right]\frac{v_0^2\left(x,z\right)}{{\omega }^2}\frac{{\partial }^2}{\partial x^2}$ , I add (4) to Equation (3) and further simplify it to obtain:

$\begin{array}{l}\frac{{\partial }^2}{\partial x^2}p\left(x,z;\omega \right)+\frac{{\partial }^2}{\partial z^2}p\left(x,z;\omega \right)+F_0\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)p\left(x,z;\omega \right)\\ =-\left(F\left(x,z,v_{qp},\epsilon ,\delta \right)-F_0\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)\right)p\left(x,z;\omega \right).\end{array}$ (5)

Equation (5) is the non-homogeneous Helmholtz equation for qP waves in VTI media, which can be rewritten as:

$\left({\nabla }^2+F_0^2\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)\right)p\left(x,z;\omega \right)=-O\left(x,z,v_{qp},\epsilon ,\delta \right)p\left(x,z;\omega \right)$ , (6)

where $O\left(x,z,v_{qp},\epsilon ,\delta \right)=F\left(x,z,v_{qp},\epsilon ,\delta \right)-F_0\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)$ is the scattering source term (slowness perturbation and anisotropy parameter perturbation), ${\nabla }^2=\partial /\partial x^2+\partial /\partial z^2$ is the Laplacian operator.

Decompose the wavefield $p\left(x,z;\omega \right)$ of the qP wave in Formula (6) into two parts: the background velocity and background anisotropic reference medium wavefield $p_0\left(x,z;\omega \right)$ , and the anisotropic perturbation and velocity perturbation scattering wavefield $p_F\left(x,z;\omega \right)$ . Then, rewrite Equation (6) as:

$\begin{array}{l}\left({\nabla }^2+F_0^2\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)\right)\left(p_0\left(x,z;\omega \right)+p_F\left(x,z;\omega \right)\right)\\ =-O\left(x,z,v_{qp},\epsilon ,\delta \right)p\left(x,z;\omega \right).\end{array}$ (7)

The wavefield in an anisotropic constant velocity background medium satisfies the homogeneous Helmholtz equation:

$\left({\nabla }^2+F_0^2\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)\right)p_0\left(x,z;\omega \right)=0$ . (8)

If I substitute Formula (8) into (7), the equation satisfied by the scattered wavefield can be written as:

$\left({\nabla }^2+F_0^2\left(x,z,v_0,{\epsilon }_0,{\delta }_0\right)\right)p_F\left(x,z;\omega \right)=-O\left(x,z,v_{qp},\epsilon ,\delta \right)p\left(x,z;\omega \right)$ . (9)

Equation (9) can be solved using the Green’s function method, resulting in:

$p_F\left(x,z;\omega \right)={\displaystyle {\int }_{\Omega }G\left(x,z;x^{\prime },z^{\prime };\omega \right)O\left(x,z,v_{qp},\epsilon ,\delta \right)p\left(x,z;\omega \right)\text{d}{x}^{\prime }\text{d}{z}^{\prime }}$ . (10)

Equation (10) is the formula for calculating the scattering field. The total wavefield of qP waves is:

$p\left(x,z;\omega \right)=p_0\left(x,z;\omega \right)+{\displaystyle {\int }_{\Omega }G\left(x,z;x^{\prime },z^{\prime };\omega \right)O\left(x,z,v_{qp},\epsilon ,\delta \right)p\left(x,z;\omega \right)\text{d}{x}^{\prime }\text{d}{z}^{\prime }}$ . (11)

Equation (11) is the Lippmann-Schwinger (L-S) equation for VTI media.

2.2. De Wolf Approximation in VTI Media

If the scattering in VTI is relatively weak, the background wavefield $p_0\left(x,z;\omega \right)$ is generally used instead of the total wavefield $p\left(x,z;\omega \right)$ . To improve the adaptability of the Born approximation in complex media, the De Wolf approximation is introduced. Firstly, simplify Equation (11):

$p=p_0+G_{0}Op$ , (12)

where O is a diagonal operator in the spatial domain, $G_0$ is a non-diagonal integral operator. If the reference medium is uniform, then $G_0$ is the integral with Green’s function $G\left(x,z;x^{\prime },z^{\prime };\omega \right)$ as the kernel.

The approximate expression for De Wolf (Wu, 1994) [17] is:

$p=p_f+G_f{O}_b{p}_f$ , (13)

where $G_f={\displaystyle \sum _{m=0}^M{\left[G_0{O}_f\right]}^m{G}_0}$ , $p_f={\displaystyle \sum _{m=0}^M{\left[G_0{O}_f\right]}^m{p}_0}$ .

2.3. Multiple Forward and One Backward Continuation Operator in VTI Media

The De Wolf approximation is a method based on multiple forward scattering and single backward scattering approximations, which ignores internal reverberation and partitions the model. $z_0$ is the entrance of the first thin plate, $z^{\prime }$ is the exit of the first thin slab and also the entrance of the second thin plate; $z^{″}$ is the outlet of the second thin slab. The position of the outgoing wave field is $\left(x^{*},z^{*}\right)$ , and the position of the incident wave is $\left(x^{\prime },z^{\prime }\right)$ . For the forward scattering field (positive direction of the Z-axis), the outgoing field is located at $\left(x^{*}=x^{″},z^{*}=z^{″}\right)$ .

By using local Born approximations in each layer of the thin slab, the scattered wavefield can be written as:

$p_s\left(x^{*},z^{*};\omega \right)=k_0^2\left(z;\omega \right){\displaystyle {\iint }_V{\tilde{G}}_0\left(x^{*},z^{*};r^{\prime };\omega \right)O\left(r^{\prime },v_{qp},\epsilon ,\delta \right)p_0\left(r^{\prime };\omega \right)\text{d}{r}^{\prime }}$, (14)

where $r^{\prime }=\left(x^{\prime },z^{\prime }\right)$ \u3002 I use the Fourier transform to transform Equation (14):

$p_s\left(k_x^{*},z^{*};\omega \right)=k_0^2\left(z;\omega \right){\displaystyle {\int }_{z^{\prime }}^{z^{*}}\text{d}z{\displaystyle \int \text{d}x{\tilde{G}}_0\left(k_x^{*},z^{*};r^{\prime };\omega \right)}}\times O\left(r^{\prime },v_{qp},\epsilon ,\delta \right)p_0\left(r^{\prime };\omega \right)$. (15)

The entrance wavefield $p_0\left(r^{\prime };\omega \right)$ of a thin slab can be written as (Wu et al., 2007) [18]:

$p_0\left(r^{\prime };\omega \right)=\frac{1}{2\pi }{\displaystyle \iint \text{d}{k^{\prime }}_x{p}_0\left({k^{\prime }}_x,z^{\prime };\omega \right){\text{e}}^{\text{i}{k}_{z0}\left(z^{\prime }\right)\left(z-z^{\prime }\right)}{\text{e}}^{\text{i}{k^{\prime }}_x{x}^{\prime }}}$ . (16)

I substitute Equation (16) into Equation (15):

$\begin{array}{c}{p}_F\left(k_x^{*},z^{*};\omega \right)={\displaystyle {\iiint }_{\Omega }G\left(k_x^{*},z^{*};r^{\prime };\omega \right)O\left(r^{\prime },v_{qp},\epsilon ,\delta \right)}\\ \times \frac{1}{2\pi }{\displaystyle \iint \text{d}{k^{\prime }}_{xy}p\left({k^{\prime }}_x,z^{\prime };\omega \right){\text{e}}^{i{k}_{z0}\left(z^{\prime }\right)}{\text{e}}^{\text{i}{k^{\prime }}_{x}x}}.\end{array}$ (17)

For convenience in calculation, the Green function of the local background medium is used here instead of the locally renormalized Green function. The expression for the Green function of the locally uniform background medium is as follows:

$G_0\left(k_x^{*},z^{*};r^{\prime };\omega \right)=\frac{\text{i}}{2k_{z0}\left(z^{\prime }\right)}{\text{e}}^{\text{i}{k}_{z0}\left(z^{\prime }\right)\left(z^{*}-z^{\prime }\right)}{\text{e}}^{-\text{i}{k}_x^{*}{x}^{\prime }}$ , (18)

where $k_{z0}\left(z^{\prime }\right)=\sqrt{F_0^2\left(r^{\prime },v_0,{\epsilon }_0,{\delta }_0\right)-{k^{\prime }}_x^2}$ . $k_{z0}\left(z^{\prime }\right)$ can be written as:

$\begin{array}{c}{k}_{z0}\left(z^{\prime }\right)=\sqrt{\frac{\left[\frac{{\omega }^2}{v_0^2\left(z^{\prime }\right)}-{\chi }_0\left(z^{\prime }\right){k^{\prime }}_x^2+{\eta }_0\left({k^{\prime }}_x,z^{\prime }\right){k^{\prime }}_x^2\right]}{{\eta }_0\left({k^{\prime }}_x,z^{\prime }\right)}-{k^{\prime }}_x^2}\\ =\frac{\omega }{v_0\left(z^{\prime }\right)}\sqrt{\frac{1-\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right)\frac{v_0^2\left(z^{\prime }\right){k^{\prime }}_x^2}{{\omega }^2}}{1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right)\frac{v_0^2\left(z^{\prime }\right){k^{\prime }}_x^2}{{\omega }^2}}},\end{array}$ (19)

where ${\eta }_0\left({k^{\prime }}_x^2,z\right)=1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right)\frac{v_0^2\left(z^{\prime }\right){k^{\prime }}_x^2}{{\omega }^2}$ .

$O\left(r^{\prime },v_{qp},\epsilon ,\delta \right)=F^2\left(r^{\prime },v_{qp},\epsilon ,\delta \right)-F_0^2\left(r^{\prime },v_0,{\epsilon }_0,{\delta }_0\right)$ can be written as:

$\begin{array}{l}O\left(r^{\prime },v_{qp},\epsilon ,\delta \right)=\frac{\left[\frac{{\omega }^2}{v_{qp}^2\left(r^{\prime }\right)}-\chi \left(r^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}+\eta \left(x^{\prime },z^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}\right]}{\eta \left(x^{\prime },y^{\prime },z^{\prime }\right)}\\ -\frac{\left[\frac{{\omega }^2}{v_0^2\left(z^{\prime }\right)}-{\chi }_0\left(z^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}+{\eta }_0\left(z^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}\right]}{{\eta }_0\left(z^{\prime }\right)}.\end{array}$ (20)

I further simplify Equation (20) as follows:

$O\left(r^{\prime },v_{qp},\epsilon ,\delta \right)=\frac{A{\eta }_0\left(z^{\prime }\right)-A_0\eta \left(x^{\prime },z^{\prime }\right)}{\eta \left(x^{\prime },z^{\prime }\right){\eta }_0\left(z^{\prime }\right)}$ , (21)

where $A=\frac{{\omega }^2}{v_{qp}^2\left(r^{\prime }\right)}-\chi \left(r^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}+\eta \left(r^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}$ $A_0=\frac{{\omega }^2}{v_0^2\left(z^{\prime }\right)}-{\chi }_0\left(z^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}+{\eta }_0\left(z^{\prime }\right)\frac{{\partial }^2}{\partial {x^{\prime }}^2}$ , to simplify the calculation, the denominator term is replaced with ${\eta }_0\left(z^{\prime }\right)$ for $\eta \left(x^{\prime },z^{\prime }\right)$ . To avoid the calculation of term $\frac{{\partial }^2}{\partial {x^{\prime }}^2}$ , it is transformed into the frequency wavenumber domain, making term $O\left({k^{\prime }}_x^2,z^{\prime },v_{qp},\epsilon ,\delta \right)$ a dual domain calculation.

Equation (21) can be simplified as follows:

$O\left({k^{\prime }}_x^2,z^{\prime },v_{qp},\epsilon ,\delta \right)=a_s\frac{\omega }{v_0\left(z^{\prime }\right)}\Delta s+a_{\epsilon }\left[\epsilon \left(r^{\prime }\right)-{\epsilon }_0\left(z^{\prime }\right)\right]+a_{\delta }\left[\delta \left(r^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right]$ , (22)

where:

$\{\begin{array}{l}{a}_s=k_0\left(z^{\prime }\right)\frac{\left[1-4\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2+2\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right)\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^4\right]}{{\left[1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2\right]}^2}\frac{\omega }{v_0\left(z^{\prime }\right)}\\ a_{\epsilon }=-\frac{2k_0^2\left(z^{\prime }\right)\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right){\lambda }^4}{{\left[1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2\right]}^2}\\ a_{\delta }=-\frac{2k_0^2\left(z^{\prime }\right)\left[1-\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right){\lambda }^2\right]{\lambda }^2}{{\left[1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2\right]}^2}\end{array}$ , (23)

where $\Delta s=v_0^2\left(z^{\prime }\right)/v_{qp}^2\left(r^{\prime }\right)-1$ , $\lambda ={k^{\prime }}_x/k_0\left(z^{\prime }\right)$ , $k_0\left(z^{\prime }\right)=\omega /v_0\left(z^{\prime }\right)$ .

By substituting Equations (18) and (22) into Equation (17), Equation (17) can be simplified as follows:

$\begin{array}{c}{p}_F\left(k_x^{*},z^{*};\omega \right)=\frac{\text{i}}{4\pi k_{z0}\left(z^{\prime }\right)}{\text{e}}^{\text{i}{k}_{z_0}\left(z^{\prime }\right)\left(z^{*}-z^{\prime }\right)}{\displaystyle \int \text{d}{k^{\prime }}_x}\times {\displaystyle {\int }_{z^{\prime }}^{z^{*}}\text{d}z{\displaystyle \int \text{d}{x}^{\prime }}}O\left({k^{\prime }}_x,z^{\prime },v_{qp},\epsilon ,\delta \right)\\ \times p_0\left({k^{\prime }}_x,z^{\prime };\omega \right){\text{e}}^{-\text{i}\left(k_{z_0}\left(z^{*}\right)-k_{z_0}\left(z^{\prime }\right)\right)z^{\prime }}{\text{e}}^{-\text{i}\left(k_x^{*}-{k^{\prime }}_x\right)x^{\prime }},\end{array}$ (24)

where $\Delta z=z^{*}-z^{\prime }$ , Equation (24) is an approximate formula for double domain thin slabs.

To improve computational efficiency, the thin slab is compressed into a screen, which compresses the three-dimensional thin slab into a two-dimensional screen to accelerate computation speed. The transverse wave numbers $k_x$ at the inlet and outlet of the thin slab are much smaller than the background wave numbers $k_0\left(z^{\prime }\right)$ inside the thin slab (Wu et al., 2007) [18].

I use Taylor expansion $k_{z0}\left(z^{\prime }\right)$ term to obtain:

$k_{z0}\left(z^{\prime }\right)=\frac{\omega }{v_0\left(z^{\prime }\right)}\sqrt{1-\frac{\left(1+2{\delta }_0\left(z^{\prime }\right)\right)v_0^2\left(z^{\prime }\right)\frac{{k^{\prime }}_x^2}{{\omega }^2}}{1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right)v_0^2\left(z^{\prime }\right)\frac{{k^{\prime }}_x^2}{{\omega }^2}}}\approx \frac{\omega }{v_0\left(z^{\prime }\right)}\left(1-0.5A^{\prime }\right)$ , (25)

where $A^{\prime }=\frac{\left(1+2{\delta }_0\left(z^{\prime }\right)\right)v_0^2\left(z^{\prime }\right){k^{\prime }}_x^2/{\omega }^2}{1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right)v_0^2\left(z^{\prime }\right){k^{\prime }}_x^2/{\omega }^2}$ .

Forward scattering $k_z\left(z;\omega \right)-{k^{\prime }}_z\left(z^{\prime };\omega \right)$ can be simplified as:

$k_z\left(z;\omega \right)-{k^{\prime }}_z\left(z^{\prime };\omega \right)\approx \left(k_0\left(z;\omega \right)-k_0\left(z^{\prime };\omega \right)\right)-\left(A^{\prime }-A\right)\approx 0$ . (26)

One backscattering $-k_z\left(z;\omega \right)-{k^{\prime }}_z\left(z^{\prime };\omega \right)$ can be simplified as:

$-k_z\left(z;\omega \right)-{k^{\prime }}_z\left(z^{\prime };\omega \right)\approx -\left(k_0\left(z;\omega \right)+k_0\left(z^{\prime };\omega \right)\right)-\left(A^{\prime }+A\right)\approx 2k_0\left(z;\omega \right)$ , (27)

when $\left(x^{*},z^{*}\right)=\left(x^{″},z^{″}\right)$ , $\left(x,z\right)=\left(x^{\prime },z^{\prime }\right)$ , multiple forward scatterings and one backward scattering can be approximated as:

$\begin{array}{c}{p}_F\left({k^{″}}_x,z^{″};\omega \right)\approx \frac{\text{1}}{2k_{z0}\left(z\right)}{\text{e}}^{\text{i}{k}_z\left(z^{\prime };\omega \right)z^{″}-\text{i}{k}_z\left(z^{\prime };\omega \right)z^{\prime }}\\ \times {\displaystyle \iint {\text{e}}^{-\text{i}{k^{″}}_x{x}^{\prime }}}\left[\text{i}O\left({k^{\prime }}_x,z^{\prime },v_{qp},\epsilon ,\delta \right)\Delta z\right]\text{d}{x}^{\prime }\\ \times \frac{1}{2\pi }{\displaystyle \iint \text{d}{k^{\prime }}_x{p}_0\left({k^{\prime }}_x,z^{\prime };\omega \right){\text{e}}^{\text{i}{k^{\prime }}_x{x}^{\prime }}},\end{array}$ (28)

$\begin{array}{c}{p}_b\left({k^{\prime }}_x,z^{\prime };\omega \right)\approx \frac{1}{2k_z\left(z^{\prime };\omega \right)}\text{sinc}\left(\Delta z\right){\text{e}}^{\text{i}{k}_0\left(z^{\prime }\right)\Delta z}{\text{e}}^{\text{i}{k}_z\left(z^{\prime };\omega \right)z^{″}-\text{i}{k}_z\left(z^{\prime };\omega \right)z^{\prime }}\\ \times {\displaystyle \iint {\text{e}}^{-\text{i}{k^{″}}_x{x}^{\prime }}}\left[\text{i}O\left({k^{\prime }}_x,z^{\prime },v_{qp},\epsilon ,\delta \right)\Delta z\right]\text{d}{x}^{\prime }\\ \times \frac{1}{2\pi }{\displaystyle \iint \text{d}{k^{\prime }}_x{p}_0\left({k^{\prime }}_x,z^{\prime };\omega \right){\text{e}}^{\text{i}{k^{\prime }}_x{x}^{\prime }}},\end{array}$ (29)

where $\text{sinc}\left(\Delta z\right)=\sin \left(\Delta z\right)/\Delta z$ .

Performing the small-angle approximation on the multiple forward and one backward operator mentioned above, specifically by using $F_s\left(x^{\prime },z^{\prime }\right)=\frac{1}{2v_0\left(z^{\prime }\right)}\left(\frac{v_0^2\left(z^{\prime }\right)}{v^2\left(x^{\prime },z^{\prime }\right)}-1\right)\approx \frac{1}{v\left(x^{\prime },z^{\prime }\right)}-\frac{1}{v_0\left(z^{\prime }\right)}$ and adding or subtracting the perturbation term $\text{i}\omega F_S\left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)$ in the expression. Therefore, the multiple forward scattering field is:

$\begin{array}{l}{p}_F\left({k^{″}}_x,z^{″};\omega \right)={\text{e}}^{\text{i}{k}_{z_0}\left(z^{\prime }\right)\Delta z}\{FT\left[\text{i}\omega F_s\left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\\ +F_{s1}\left({k^{\prime }}_x,z^{\prime }\right)FT\left[\text{i}{F}_s\left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\\ +F_{\epsilon }\left({k^{\prime }}_x,z^{\prime }\right)FT\left[\text{i}\Delta \epsilon \left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\\ +F_{\delta }\left({k^{\prime }}_x,z^{\prime }\right)FT\left[\text{i}\Delta \delta \left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\}.\end{array}$ (30)

The backscattering field is:

$\begin{array}{l}{p}_b\left({k^{\prime }}_x,z^{\prime };\omega \right)=\text{sinc}\left(\Delta z\right){\text{e}}^{\text{i}{k}_0\left(z^{\prime }\right)\Delta z}{\text{e}}^{\text{i}{k}_{z_0}\left(z^{\prime }\right)\Delta z}\\ \times \{FT\left[\text{i}\omega F_s\left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\\ +F_{s1}\left({k^{\prime }}_x,z^{\prime }\right)FT\left[\text{i}{F}_s\left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\\ +F_{\epsilon }\left({k^{\prime }}_x,z^{\prime }\right)FT\left[\text{i}\Delta \epsilon \left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\\ +F_{\delta }\left({k^{\prime }}_x,z^{\prime }\right)FT\left[\text{i}\Delta \delta \left(x^{\prime },z^{\prime }\right)\Delta z{p}_0\left(x^{\prime },z^{\prime };\omega \right)\right]\}.\end{array}$ (31)

Among them, $FT$ represents the Fourier transform, $F{T}^{-1}$ represents the inverse Fourier transform, and the expressions for $F_{s1}\left({k^{\prime }}_x,z^{\prime }\right)$ , $F_{\epsilon 1}\left({k^{\prime }}_x,z^{\prime }\right)$ , and $F_{\delta 1}\left({k^{\prime }}_x,z^{\prime }\right)$ are as follows:

$\{\begin{array}{l}{F}_{s1}\left({k^{\prime }}_x,{k^{\prime }}_y,z^{\prime }\right)=\frac{k_0\left(z^{\prime }\right)}{k_{z0}\left(z^{\prime }\right)}\frac{\left[1-4\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2+2\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right)\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^4\right]}{{\left[1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2\right]}^2}\omega -\omega \\ F_{\epsilon 1}\left({k^{\prime }}_x,{k^{\prime }}_y,z^{\prime }\right)=-\frac{k_0^2\left(z^{\prime }\right)}{k_{z0}\left(z^{\prime }\right)}\frac{\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right){\lambda }^4}{{\left[1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2\right]}^2}\\ F_{\delta 1}\left({k^{\prime }}_x,{k^{\prime }}_y,z^{\prime }\right)=-\frac{k_0^2\left(z^{\prime }\right)}{k_{z0}\left(z^{\prime }\right)}\frac{\left[1-\left(1+2{\epsilon }_0\left(z^{\prime }\right)\right){\lambda }^2\right]{\lambda }^2}{{\left[1-2\left({\epsilon }_0\left(z^{\prime }\right)-{\delta }_0\left(z^{\prime }\right)\right){\lambda }^2\right]}^2}\end{array}$ . (32)

I set $\left(x^{*},z^{*}\right)=\left(x^{″},z^{″}\right)$ , and the entire forward scattering field $p\left(r^{″};\omega \right)$ can be expressed as:

$p\left(r^{″};\omega \right)=p_0\left(r^{″};\omega \right)+p_s\left(r^{″};\omega \right)+p_{\epsilon }\left(r^{″};\omega \right)+p_{\delta }\left(r^{″};\omega \right)$ , (33)

where

$\{\begin{array}{l}{p}_0\left(r^{″};\omega \right)=F{T}^{-1}\left\{{\text{e}}^{\text{i}{k}_{z0}\left(z^{\prime }\right)\Delta z}FT\left[{\text{e}}^{\text{i}\omega F_s\left(r^{\prime }\right)\Delta z}{p}_0\left(r^{\prime };\omega \right)\right]\right\}\\ p_s\left(r^{″};\omega \right)=F{T}^{-1}\left\{{\text{e}}^{\text{i}{k}_{z0}\left(z^{\prime }\right)\Delta z}{F}_{s1}\left({k^{\prime }}_x,z^{\prime }\right)FT\left[{\text{e}}^{\text{i}\omega F_s\left(r^{\prime }\right)\Delta z}\text{i}\omega \Delta s\left(r^{\prime }\right)\Delta z{p}_0\left(r^{\prime };\omega \right)\right]\right\}\\ p_{\epsilon }\left(r^{″};\omega \right)=F{T}^{-1}\left\{{\text{e}}^{\text{i}{k}_{z0}\left(z^{\prime }\right)\Delta z}{F}_{\epsilon 1}\left({k^{\prime }}_x,z^{\prime }\right)FT\left[{\text{e}}^{\text{i}\omega F_s\left(r^{\prime }\right)\Delta z}\text{i}\Delta \epsilon \left(r^{\prime }\right)\Delta z{p}_0\left(r^{\prime };\omega \right)\right]\right\}\\ p_{\delta }\left(r^{″};\omega \right)=F{T}^{-1}\left\{{\text{e}}^{\text{i}{k}_{z0}\left(z^{\prime }\right)\Delta z}{F}_{\delta 1}\left({k^{\prime }}_x,z^{\prime }\right)FT\left[{\text{e}}^{\text{i}\omega F_s\left(r^{\prime }\right)\Delta z}\text{i}\Delta \delta \left(r^{\prime }\right)\Delta z{p}_0\left(r^{\prime };\omega \right)\right]\right\}\end{array}$ , (34)

where $\Delta z=z^{″}-z^{\prime }$ ，$z^{\prime }$ is the depth position of the incident wavefield, and $z^{″}$ is the depth position of the outgoing wavefield. Equation (34) is a generalized screen approximation extension operator for VTI media based on the De Wolf approximation.

However, due to the existence of term $1/k_{z0}$ , the operator exhibits spatial aliasing effects, making terms $v_0^2\left(z\right)k_x^2/{\omega }^2={\sin }^2\theta$ and $k_{z0}$ equivalent to:

$\begin{array}{c}{k}_{z0}=\frac{\omega }{v_0\left(z\right)}\sqrt{\frac{1-\left(1+2{\epsilon }_0\left(z\right)\right){\sin }^2\theta }{1-2\left({\epsilon }_0(z)-{\delta }_0\left(z\right)\right){\sin }^2\theta }}\\ =\frac{\omega }{v_0\left(z\right)}\sqrt{1-\frac{\left(1+2{\delta }_0\left(z\right)\right){\sin }^2\theta }{1-2\left({\epsilon }_0(z)-{\delta }_0\left(z\right)\right){\sin }^2\theta }}\\ =\frac{\omega }{v_0\left(z\right)}\sqrt{1-\frac{1+2{\delta }_0\left(z\right)}{1/{\sin }^2\theta -2\left({\epsilon }_0(z)-{\delta }_0\left(z\right)\right)}}.\end{array}$ (35)

As the angle approaches 90 degrees, $\sin \theta$ gradually approaches 1, then:

$\frac{1+2{\delta }_0\left(z\right)}{\frac{1}{{\sin }^2\theta }-2\left({\epsilon }_0\left(z\right)-{\delta }_0\left(z\right)\right)}=\frac{1+2{\delta }_0\left(z\right)}{\left(1+2{\delta }_0\left(z\right)\right)-2{\epsilon }_0\left(z\right)}$ , (36)

when ${\epsilon }_0\left(z\right)\ge 0$ , $\frac{1+2{\delta }_0\left(z\right)}{\left(1+2{\delta }_0\left(z\right)\right)-2{\epsilon }_0\left(z\right)}\ge 1$ . Therefore, when the propagation angle approaches 90 degrees, the $1/k_{z0}$ term will still exhibit an unstable phenomenon, namely the spatial aliasing phenomenon. Referring to the isotropic medium processing method, $M^{\prime }$ is introduced here to change the filtering characteristics of the operator on the dip angle of the formation. The improved $1/k_{z0}$ is replaced by $M$ :

$\begin{array}{c}M=\frac{1}{k_{z0}}\cdot M^{\prime }=\frac{1}{k_{z0}}\frac{k_{z0}}{\sqrt{b^2{k}_x^2+k_{z0}^2}}\\ =\frac{1}{\sqrt{b^2{k}_x^2+k_{z0}^2}}=\frac{1}{\sqrt{{\omega }^2/v_0^2\left(z\right)+\left(b^2-c\right)k_x^2}},\end{array}$ (37)

where $c=\sqrt{\frac{{\omega }^2}{v_0^2\left(z\right)}\frac{1+2{\delta }_0\left(z\right)}{1/{\sin }^2\theta -2\left({\epsilon }_0\left(z\right)-{\delta }_0\left(z\right)\right)}}$ . When $b=c$ occurs, spatial aliasing disappears, but it is often not processed cleanly during actual calculations. When $b>c$ occurs, it can further suppress the spatial aliasing on both wings, but also remove some effective signals. When $b<c$ occurs, the spatial aliasing effect is more severe. The longitudinal amplitude attenuation factor is:

$l\left(z\right)=1-\left[1-s_z\right]\sin \left(\frac{\pi }{2}\frac{j}{N-1}\right)$ . (38)

Among them, $l\left(z\right)$ is the correction factor, which takes different values at different depths; N is the grid point in the vertical direction, with specific values of $j=0,1,\cdots ,N-1$ ; $s_z$ is the attenuation coefficient, with slightly different values at different depths. Assume the top layer is 1.0 and the bottom layer is a very small value.

3. Algorithm Procedure

The forward simulation process of the De Wolf approximation (thin slab approximation and generalized screen) in a VTI medium is as follows:

1) Divide the entire underground space into a series of horizontal thin slabs (with a thickness of one grid).

2) In each thin slab (screen), the velocity is represented as background velocity $v_0\left(z\right)$ and $v_{\text{per}}\left(z\right)$ disturbance parameters, and the anisotropic parameters are represented as background anisotropic parameters ${\epsilon }_0\left(z\right)$ and ${\delta }_0\left(z\right)$ , and anisotropic disturbances ${\epsilon }_{\text{per}}\left(z\right)$ and ${\delta }_{\text{per}}\left(z\right)$ , respectively.

3) Place the seismic source wavelet function on the surface and perform a Fourier transform to convert it into the frequency domain.

4) Convert the wavefield at the entrance of the thin slab (screen) to the frequency wavenumber domain, propagate it with background parameters ($v_0\left(z\right)$ , ${\epsilon }_0\left(z\right)$ , ${\delta }_0\left(z\right)$ ), and calculate the forward wavefield $p_0\left(x,z;\omega \right)$ inside the thin slab (screen).

5) Convert all wavefields inside the thin plate (screen) into the frequency space domain and perform spatial disturbance correction. Forward scattered or transmitted waves are mainly controlled by velocity disturbances and anisotropic parameters, while backscattered or reflected waves are affected by impedance disturbances caused by velocity disturbances and anisotropic parameters.

Take the total forward scattering field as the initial condition for the next thin plate, while retaining the primary backward scattering wave field.

7) Repeat steps 4) - 6) until the bottom of the model to obtain all forward and backward wavefields.

8) Using the backscatter field at the bottom as the boundary condition for the first upward thin slab (screen), propagate toward the surface direction.

9) Calculate according to steps 4) - 6), and add the backscattered field saved during the forward calculation process as the total field for the next thin plate (screen). Starting from the bottom of the model’s thin slab (screen), accumulate the scattered field of each thin slab (screen) interface in the opposite vertical direction until reaching the surface, and obtain scattered wave seismic records.

4. Numerical Simulation

4.1. Depression Model

Figure 1(a) is a two-dimensional velocity model. The grid points of the model in the x and z directions are 401 and 301, respectively; the grid spacing is 5.0 m; the VTI (vertically transversely isotropic) medium is a transversely isotropic model characterized by a vertical axis of symmetry. It approximates the sedimentary stratigraphic characteristics of subsurface layers. Compared to other anisotropic models, the propagation characteristics of seismic waves in the VTI medium are more easily described. The anisotropic parameters $\epsilon \left(x,z\right)$ and $\delta \left(x,z\right)$ have the following linear relationship with velocity $v\left(x,z\right)$ (Huang et al., 2017) [27]:

$\{\begin{array}{l}\epsilon \left(x,z\right)=\frac{0.606\left[v\left(x,z\right)-v_{\min }\right]}{v_{\max }}\\ \delta \left(x,z\right)=\frac{0.485\left[v\left(x,z\right)-v_{\min }\right]}{v_{\max }}\end{array}$ . (39)

Among them, $v_{\max }$ and $v_{\min }$ are the maximum and minimum velocities, respectively, and the calculation results are shown in Figure 1(b) and Figure 1(c). Based on the above model, the parameter settings of the observation system are as follows: there are a total of 200 detectors, with a track spacing of 5.0 meters, arranged in the center of the detectors, and the first detector is located at (0.5 km, 0.0). The seismic source wavelet adopts Rayleigh waves with a dominant frequency of 30 Hz, a sampling interval of 1 s, and a sampling length of 1600.

Figure 2 shows the forward modeling results. Figure 2(a) is the generalized screen seismic record of isotropic media (Acoustic-thin slabs, without false frequencies), and Figure 2(b) and Figure 2(c) are the approximate seismic records of VTI media thin slabs (VTI-thin slabs) without and with spatial aliasing, respectively. Figure 2(d) and Figure 2(e) show the approximate seismic records of VTI media screens (VTI-screens) without and with false frequencies, respectively. As shown in the figure, seismic records with spatial aliasing have relatively weak amplitudes of deep reflection waves. The use of amplitude attenuation factors can effectively eliminate the phenomenon of spatial aliasing, and seismic records with spatial aliasing effects have a higher signal-to-noise ratio.

Figure 1. Concave model parameters: (a) Velocity model parameter; (b) $\epsilon$ parameter; (c) $\delta$ parameter.

To verify the effectiveness of the algorithm, I used the traditional visco-acoustic medium finite difference method for comparison. Figure 3(a) shows the finite difference forward simulation results. Compared with the thin slab approximation and screen approximation methods, the seismic records of the three methods are similar. The finite difference method has richer wavefield information, but the phase axis is relatively coarser. Extract one of the waveforms for waveform comparison (as shown in Figure 3(b)). The phase difference between the acoustic thin slab approximation result and the seismic records of the other three methods is relatively large. Ignoring the anisotropy of the medium directly leads to a significant lag in the seismic records of the acoustic thin slab approximation compared to the other three methods, and the wavefield components are relatively simple. The thin slab approximation and screen approximation results are closer to each other. Due to the relatively simple model, the two algorithms’ small-angle situations yield similar results. The above research results verify the effectiveness of the algorithms. Table 1 shows the calculation time of several methods. The calculation time for the VTI medium is higher than that for the acoustic medium, due to the addition of anisotropic parameter calculation. The efficiency of the thin slab approximation calculation under the same conditions is lower than that of the screen approximation method. The reason is that screen approximation is an approximation method for thin slab approximation at small angles, and the calculation efficiency of both methods is higher than that of the traditional finite difference method.

Figure 2. Seismic records: (a) Acoustic-thin slabs (without spatial aliasing); (b) VTI-thin slabs (with spatial aliasing); (c) VTI-thin slabs (without v); (d) VTI-screen (with spatial aliasing); (e) VTI-screen (without spatial aliasing).

Figure 3. Comparison of calculation results of concave model: (a) Finite difference forward modeling results of VTI medium; (b) Comparison of waveforms.

Table 1. Comparison of calculation time by different methods (computer equipment CPU: Intel® Core™ i7-4800MQ CPU @ 2.7 GHz; Memory: 16 GB).

Method	Acoustic-thin slabs (without)	VTI-thin slabs (with)	VTI-thin slabs (without)	VTI-screens (with)	VTI-screens (without)	VTI-FD
Time	27 s	46 s	47 s	40 s	41 s	153 s

4.2. Complicated Models

To further validate the effectiveness of the algorithm, a velocity model is set up as shown in Figure 4(a), with 401 and 301 grid points in the x and z directions, respectively; the grid spacing is 5.0 m; the anisotropy parameters are shown in Figure 4(b) and Figure 4(c). Based on the above model, the parameter settings of the observation system are as follows: there are a total of 200 detectors, with a track spacing of 5.0 meters, arranged in the center of the detectors, and the first detector is located at (0.5 km, 0.0); the seismic source wavelet adopts Rayleigh waves with a dominant frequency of 30 Hz, a sampling interval of 1 s, and a sampling length of 1600.

Figure 4. Complex model parameters: (a) Velocity model parameter; (b) $\epsilon$ parameter; (c) $\delta$ parameter.

I used the thin slab approximation method and the screen approximation method for forward simulation, and Figure 5 shows the seismic forward simulation records. Among them, Figure 5(a) shows the generalized screen seismic record of isotropic media (Acoustic-thin slabs, without spatial aliasing), and Figure 5(b) shows the seismic record results of the global Born approximation of VTI media (VTI-Born, without spatial aliasing); Figure 5(c) and Figure 5(d) show the approximate seismic records of VTI medium thin slabs (VTI-thin slabs) without and without spatial aliasing removal, respectively. Figure 5(e) and Figure 5(f) show the approximate seismic records of VTI media screens (VTI-screens) without and with spatial aliasing, respectively. As shown in the figure, when the cumulative effect of forward scattering is strong, the calculation error of the global Born approximation is relatively large. The thin slab approximation and screen approximation based on the De Wolf approximation can effectively solve this problem. When the amplitude of deep reflection waves in seismic records with spatial aliasing is relatively weak, the use of an amplitude attenuation factor can effectively eliminate the phenomenon of spatial false frequency, and the signal-to-noise ratio of seismic records without the spatial false frequency effect is higher.

Figure 5. Seismic records: (a) Acoustic-thin slabs (without spatial aliasing); (b) VTI-Born (without spatial aliasing); (c) VTI-thin slabs (with spatial aliasing); (d) VTI-thin slabs (without spatial aliasing); (e) VTI-screen (with spatial aliasing); (f) VTI-screen (without spatial aliasing).

To verify the effectiveness of the algorithm, I used the traditional viscous acoustic medium finite difference method for comparison. Figure 6(a) shows the finite difference forward simulation results. Compared with the thin slab approximation and screen approximation methods, the seismic records of the three methods are similar. The finite difference method has richer wavefield information, but the phase axis is relatively coarser. One of the waveforms is extracted for waveform comparison (as shown in Figure 6(b)). The phase difference between the acoustic thin slab approximation results and the seismic records of the other three methods is relatively large. Ignoring the anisotropy of the medium directly leads to a significant lag in the seismic records of the acoustic thin slab approximation compared to the other three methods, and the wavefield components are relatively simple. The results of the thin slab approximation and screen approximation are closer to each other. Due to the relatively simple model, the two algorithms have similar results in small-angle situations. The above research results validate the effectiveness of the algorithm.

Figure 6. Comparative analyses: (a) Finite difference forward modeling results of VTI medium; (b) Comparison of waveforms.

5. Conclusion

The underground sedimentary strata are mostly distributed in layers and exhibit lateral isotropy. VTI media are an approximate representation of the anisotropy characteristics of actual Earth media. Research has shown that scattered waves carry geometric and physical information related to complex structures and rock types. However, the traditional Born approximation is a weak scattering approximation and is not suitable for long-distance propagation. The De Wolf approximation can effectively overcome this problem. To achieve the De Wolf approximation, Wu proposed the thin slab approximation and the screen approximation in the 1990s and provided corresponding implementation algorithms. This method is widely used in the numerical simulation of seismic waves. However, there have been no reports on using the De Wolf approximation to simulate the propagation process of seismic waves in VTI media. Based on this, to overcome the problem of calculation errors in the Born approximation when the forward cumulative effect is strong in VTI media, this paper combines the De Wolf approximation method with VTI media and proposes a numerical simulation method for scattering waves in VTI media based on the De Wolf approximation. The scattering equations are divided into forward and backward scattering equations, and the Born series is renormalized to improve its convergence in VTI media. During the calculation process, multiple interlayer waves are ignored, and only multiple forward and one backward scattering signals are retained. Due to the influence of the anisotropic parameters of the medium, there may be phase deviations in the forward modeling records of seismic waves using the thin slab approximation method and the screen approximation method; because the thin slab approximation method and the screen approximation method only account for one reflection, there are no signals of multiple waves in the seismic records, and the comparative results verify the effectiveness of the algorithm. When the velocity disturbance of the medium along the depth direction is relatively weak (when the non-uniformity is relatively weak), the thin plate can be relatively thicker; when the medium disturbance is relatively strong, thin plates can be thinner; the more thin plates are divided (with smaller thickness), the lower the calculation efficiency and the relatively higher the calculation accuracy; the fewer thin plates are divided (thicker), the higher the calculation efficiency, but the calculation accuracy is relatively low. The key to selecting the background velocity is to ensure that the overall velocity disturbance inside the thin plate is as small as possible.