It is demonstrated in the study that most sedimentary rocks on the crust are showing seismic anisotropy   . Therefore, the study of seismic wave propagation in anisotropic media is an important subject and research trend in the field of seismology and exploration seismology. Numerical simulation technique in the seismic wave field is based on the theory of seismic wave propagation, simulating the propagation process of seismic wave in an underground media. It is one of the important methods to study and understand the seismic anisotropy. At present, the commonly-used numerical simulation methods include ray tracing method  and wave equation method. The wave equation method includes the Pseudo Spectrum Method (PSM)   , Finite Element Method(FEM)  , Boundary Element Method (BEM)  , Spectral Element Method  , finite difference method    and so on. The numerical simulation of the wave equation in the anisotropic medium has been developed rapidly all over the world and gradually improved in the 1980s  -  . It includes methods such as pseudo-spectrum method and finite difference method. The grid is developed from conventional grid to staggered grid, rotating staggered grid, and so on. The applied wave equations are elastic wave equation and quasi-P wave equation.
With the development of numerical simulation technology and the improvement of production practice requirements, the researchers put forward the higher requirements for the precision and computational efficiency of finite difference numerical simulation. The most straightforward way to improve the accuracy of differentials is to increase the number of grid nodes in the differential calculation, which also increases the computational effort and the required memory space, and that is not conducive to boundary processing. The approximate analytic discretization method can solve this contradiction well and construct a high-order interpolation function in space to approximate the spatial partial derivative by using the truncated Taylor formula to obtain an optimal approximation of the original function     . In this method, only three nodes are needed, which can make the differential precision of the partial derivative of the space reach the sixth order precision of the traditional difference method, and can greatly save the memory requirement and the computation amount. It is of great significance for the large 3D scale model forward and inversion modeling. In this paper, the two-dimensional elastic wave equation of VTI medium is simulated by this method.
2. VTI Medium Elastic Wave Equation
Seismologists believe that the anisotropy of fluctuation in the underground rocks is ubiquitous. The characteristics and laws of seismic wave propagation in the underground are related to the propagation direction. The anisotropy is mostly caused by crack-induced anisotropy and cyclical thin interbeds, both of which are hexagonal anisotropic media, also known as Transversely Isotropy (TI). It is a very vital anisotropic media model. Practically about 70% of the sedimentary rocks show TI anisotropy. According to the axis of symmetry of the TI medium, it can be divided into a transversely isotropic medium with a vertical symmetry axis (VTI) and a transversely isotropic medium with a horizontal axis of symmetry (HTI). The elastic matrix of the VTI medium has five independent elastic coefficients, N/m2, and the elastic coefficient matrix can be expressed as the formula (1).
The elastic coefficient matrix of the VTI medium represented by the formula (1) is taken into the formulas (2)-(5) to obtain the elastic wave equation of the 3-D VTI medium.
1) A differential equation of Motion describing the relationship between displacement and stress:
2) A Cauchy equation describing the relationship between strain and displacement:
3) The constitutive equation describing the relationship between stress and strain:
Among them, ρ is the density of the medium, kg/m3; is the displacement vector, m; is the physical density vector, N/m3; ε is the strain tensor, 1; σ is the stress tensor, N/m2; L is the partial derivative operator matrix:
The elastic wave equations of the obtained three-dimensional VTI medium are expressed as formula (6):
Among them, .
If only the two-dimensional elastic waves in the XOZ plane is studied, variables related to y in the above formula and the partial derivative of the variable y are set zero, in order to obtain the simplified elastic wave equations of the two-dimensional VTI medium:
It can be seen from formula (7) that the displacement of particles in the two-dimensional VTI medium is related to the four elastic coefficients (c11, c33, c44, c13), but the four elastic coefficients are only the proportion coefficients that describe the relationship between stress and strain, and have no intuitive physical meaning. It cannot directly reflect the characteristics of seismic wave propagation, especially the seismic wave propagation speed. For the sake of theoretical research and practical application, Thomsen (1986) proposed a set of parameters describing the elastic properties of TI media, defined as follows:
The above five parameters are called Thomsen parameters, and their physical meanings are: 1) vP0 and vS0 respectively denote the phase velocity of the quasi-P and quasi-S waves in the direction of the axis of symmetry (i.e., the vertical direction), m/s. Because the polarization direction of P wave and SV wave particle in the VTI medium is not parallel and perpendicular to the propagation direction of the wave, it is called quasi-P wave (qP wave) and quasi-S wave (qS wave). 2) ε, γ and δ are the three dimensionless factors representing the intensity of the seismic anisotropy of the medium, 1. When they are equal to zero, the medium is isotropic. 3) ε and γ respectively denote the anisotropic intensity of the qP wave and qS wave. The larger they are, the more obvious anisotropy of the qP wave and qS wave is, the greater the difference between the vertical velocity and the horizontal velocity is, and when they are equal to zero, the qP wave and qS wave are anisotropic. In general, the monotonicity of ε and γ is consistent, and they are incremented or equal to zero at the same time. 4) δ indicates the variation of the qP wave’s anisotropic intensity in the main direction of the vertical elasticity, and it is a transitional parameter between vP0 and vS0, which affects the wave velocity in the vicinity of the symmetry axis of the VTI medium.
From the Thomsen parameter, the qP wave and qS wave propagation velocity in the XOZ plane are controlled by the elastic parameters : 1) , while the magnitude of the phase velocity of the qP wave is controlled by in the horizontal and vertical directions. 2) The magnitude of the phase velocity of the qS wave is controlled by in the horizontal and vertical directions, and the phase velocity is equal in both directions. 3) In other horizontal and vertical directions, the phase velocity of the qP wave and qS waves and shape of the seismic wave are affected by .
Thomsen parameters have more definite physical meanings, and by Thomsen parameters, the four elastic coefficients in formula (7) can be expressed as:
3. ONAD Numerical Simulation Method and VTI Medium Elastic Wave Equation Differential Scheme
As a finite difference numerical simulation method, the optimal nearly-analytic discrete method closes the partial derivative in space by using the truncated Taylor formula to construct the high-order interpolation function in time and space to obtain an optimal approximation of the original function. The characteristic of this method is that it includes the different orders partial derivative of the displacement variable when the partial differential equation is solved. The correlation is noted between the original function and the partial derivative of different orders. The method can effectively reduce the loss of the seismic information in the discrete process and improves the accuracy of numerical calculation and the effectiveness of the calculation  .
Following the difference equation of the optimal nearly-analytic discrete method (ONAD), the numerical simulation of the 2-D VTI medium elastic wave equation will be deduced.
Let the solution satisfy the wave Equation (7) to be U, it is defined below:
where u and w represent the magnitude of the x and z displacements, m, ux, uz, wx and wz respectively indicate the gradient of the displacements u and w in the x and z directions, 1.
By Taylor formula, two formulas can be obtained as follows:
By adding and transforming the above two formulas, the following formula can be obtained:
The above equation is the time-layer recursive formula of the ONAD method, and its truncation error is , Where i and j represent the space grid point coordinates of Un−1; n and n + 1 denote the time grid coordinates of U, respectively denoting the time layer n − 1, n, and n + 1. P and the second time derivative of P in time can be transformed into the high-order partial derivative of U in space by the wave Equation (7), for example:
It can be seen from the formulas (14)-(18) that the value of the time layer n+1 can be calculated only by knowing the sum of the displacement U and the value and of U’s first order space partial derivative at the time layer n and n-1. In the formulas (15)-(18), there are higher order derivatives of
the displacement U for the space, for example: , and , etc.
According to the idea of nearly-analytic discrete method, these spatial derivatives are approximated by the displacement and gradient of the surrounding points of the space grid. Due to space constraints, only formula (19) is given here, and the approximate function of the other higher order partial derivatives is given in the literature  .
4. Wave Field Simulation and Feature Analysis
In the previous section, the basic principles of the ONAD method was discussed and the ONAD difference scheme was derived for the elastic wave equation of the 2-D VTI medium. In this section, the uniform model and the layered model are calculated by using the difference formulas (14)-(18). The wave field characteristics of the elastic wave in the VTI medium are analyzed, and the influence of the Thomsen parameter on the propagation law of the seismic wave is further analyzed.
4.1. Uniform VTI Media Model
It can be seen from the formula (9) that the qP wave and qS wave in the two-dimensional XOZ plane are related to the four parameters of vP0, vS0, ε and δ, vP0 and vS0, respectively representing the vertical phase velocity of the qP wave and qS wave. The medium’s anisotropic features are characterized by ε and δ. In order to analyze the influence of ε and δ on seismic wave field characteristics, the uniform medium model is set, the model’s length and depth are both 3000 m, and the vertical velocity of qP wave and qS wave is constant which are 3000 m/s and 1732 m/s respectively. The density of model is 2.29 g/cm3, and the parameters are used in Table 1.
Table 1. Thomsen parameter table for homogeneous media model.
At the center of the model (X = 1500 m, Z = 1500 m), the w component is loaded with the vertical concentrated force source of Ricker wavelet whose dominant frequency is 25 Hz. Taking the time step as 1 ms, the spatial grid is 10 m. Figure 1 is wave field simulation snapshot of displacement u in x direction at 320 ms when the VTI media is calculated using the parameters in Table 1.
As can be seen from the wave field snapshots in Figure 1:
1) Compared with the qS wave, the wavefront variation of the qP wave is not very intense; it is mainly controlled by ε, which can be called the elliptical factor. When ε is negative, the vertical velocity of the qP wave is greater than the horizontal phase velocity, and the wavefront of the qP wave is the upright ellipse with the long axis as the Z axis. When ε is equal to zero, the vertical velocity of the qP wave is equal to the horizontal phase velocity, but by the influence of δ, the wavefront of the qP wave is approximately circular. When ε is positive, the vertical velocity of the qP wave is less than the horizontal phase velocity, and the wavefront of the qP wave is the lateral ellipse of the X axis. When both ε and δ are zero, the medium becomes an isotropic medium, such as the third of the third group, and the wavefront of the qP wave and qS wave is circular.
2) Although the vertical and transverse phase velocities of qS waves are the same in the 2-D VTI medium, they are affected by ε and δ, which causes the wavefront of the qS wave to change severely, showing a variety of shapes, such as circular, square and diamond, and producing the tricorn region.
3) With the change of ε and δ, there will be amplitude singularity in the wave field tricorn zone before the wave front of qS wave, and the position and degree of the tricorn zone are related to the relationship between ε and δ size. The greater the difference between ε and δ is, the more obvious the tricorn zone is, such as the first group of wave field snapshot. The tricorn area has become more and more obvious, while the fifth group is the opposite. Regardless of the size of
Figure 1. Wave field snapshot of U by the elastic wave equation numerical simulation in 2-D homogeneous VTI media model (320 ms). (a) the first group ; (b) the second group ; (c) the third group ; (d) the fourth group ; (e) the fifth group .
the ε and δ, when they are equal, there is no tricorn zone. In addition, it can be seen that the more obvious the qS wave tricorn zone is, the more the wave front of the qS wave deviates from the ellipse, so can be called the non-elliptic factor of the qP wave or qS wave tricorn zone factor.
4) When , the tricorn zone appears between the vertical and horizontal symmetry axis, which is diagonal direction, and the qS wave front is close to the square, such as the first group. When , the tricorn zone appears in the direction of symmetry axis, and the qS wave front is close to the diamond, such as the fifth group. When the tricorn area is not obvious, the wave front is close to the circle.
4.2. Layered VTI Media Model
In order to study the propagation law of elastic wave in layered VTI medium, two layers of horizontal layered media model are set up. The thickness of each layer is 500 m and the model length i is 1000 m.s The model parameters are shown in Table 2. The first layer is the VTI medium, while the second layer is the isotropic medium. The Thomsen parameter of the first layer is changed, and the propagation characteristics of the elastic wave in the layered VTI medium are analyzed by the ground seismic record.
Table 2. Thomsen parameter table for layered media model
In the middle of the ground (X = 500 m, Z = 0 m) to stimulate the equal energy source of the Ricker wavelet whose dominant frequency is 25 Hz.The space step of the vertical and horizontal is 8 m, the time step is 1 ms, and the boundary is PML artificial boundary condition   . Figure 2 shows the wave field simulated snapshots of the displacement u in the X-direction at 320 ms.
The wave field snapshots of Figure 2 show: 1) Only when δ and ε are equal to 0, the wave front of each wave is a circle, and the rest are in the form of an ellipse or a diamond. 2) When Thomsen parameters are different, the wave field snapshots are not the same. With the change of the Thomsen parameter, the wavefront variation of the qP wave are more severe than that of the qS wavefront, and the qS wave emerges from the wavefront sharp corners in some cases, such as, When . They are caused by velocity anisotropy. 3) The U component of the particle displacement polarity reversal is on the Z-direction symmetrical axis at the point of the source. 4) When the wave reaches the
Figure 2. Wavefield snapshot of U by the elastic wave equation numerical simulation in 2D layered VTI media (320 ms). (a) the first group ; (b) the second group ; (c) the third group ; (d) the fourth group ; (e) the fifth group .
(1. Reflected qP wave, 2. Reflected converted qS wave, 3. Reflected converted qP wave, 4. Reflected qS wave, 5. Transmitted qS wave, 6. Transmitted converted qP wave, 7. Transmitted converted qS wave, 8. Transmitted qP wave, 9. Wavefront sharp corners).
layered interface that is the impedance interface, the reflection and transmission will occur, and converted waves will produce, such as reflected converted qP wave and transmitted converted qS wave. Sometimes the various waves interfere with each other, so they are difficult to be distinguished, as shown in Figure 3 that because the time difference is very small, so the reflected converted qP wave and reflected qS wave are stacked together.
(1. direct qP wave, 2. direct qS wave, 3. reflected qP wave, 4. reflected transmitted qS wave, 5. reflected transmitted qP wave, 6. reflected qS wave).
Figure 3. Seismic record of U by the elastic wave equation numerical simulation in 2D layered VTI media (750 ms). (a) the first group ; (b) the second group ; (c) the third group ; (d) the fourth group ; (e) the fifth group .
Figure 3 shows the horizontal component seismic records received on the ground. The length of the record is 0.75 s. The longitudinal and transverse wave field in the record is clear, but the numerical dispersion and artificial boundary reflection almost cannot be seen. The direct wave and other wave and polarity reversal phenomena can be seen in the record. It has great value to study the propagation law of seismic waves in VTI media.
In order to show clearly, the seismic records demonstrate the gain processing. It can be seen from the seismic records in Figure 3 that the direct wave and the reflected wave of the stratigraphic interface show the following characteristics as the changes of Thomsen in the first layer of the medium: 1) Due to the influence of velocity anisotropy, the speed of the direct wave propagating along the horizontal direction is different when the Thomsen parameter is changed, so the time of the direct wave arriving at each detector is changed. For example, the horizontal phase velocity of qP wave is mainly affected by ε, and when ε becomes larger, the horizontal velocity of qP wave increases and the slope of the phase axis decreases. 2) That the U component of the particle displacement polarity reversal on the Z-direction symmetrical axis at the point of the source is the same as the wave field snapshot. 3) With the change of Thomsen parameter, the change of qP wave is not severe than that of the qS wave. With the increase of ε, the reflected qS wave becomes more gentle and the time of arriving the detector changes obviously. It shows that it is more advantageous to use the transverse wave to predict the anisotropy of the formation than the longitudinal wave in the actual seismic exploration.
It is one of the important problems in seismic exploration to study the propagation law and characteristics of seismic waves in anisotropic media. In order to obtain higher simulation accuracy in conventional finite difference numerical simulations, it is necessary to use more nodes. This will not only increase the amount of computation and storage space, but also increase the difficulty of artificial boundary processing.
In this paper, the optimized approximate discretization method was used. Only three nodes were used, and the differential precision of the partial derivative of the space reached the sixth order precision of the traditional difference method. This method can greatly save the memory requirement and the computation. The difference scheme of the fourth order and space sixth order VTI medium wave equation was also established in this paper, and the numerical simulation of the uniform medium and the layered media model was carried out. The simulation results show that the seismic wave field is clear and the applicability of the method is explained.
The propagation law of seismic wave in VTI medium was analyzed in detail by the wave field snapshot and ground seismic record of uniform and layered media model. It is concluded that the ε value is sensitive to the anisotropy of qP wave, the δ value is sensitive to the anisotropy of qS wave, and the negative value δ of the anisotropic characteristics is very obvious. In addition, there will be a trigeminal zone. Meanwhile, the value of Thomsen parameter determines the phase velocity and group velocity of various waves, which also determines the spatial position and geometry of the wavefront. Therefore, in seismic exploration, accurate acquisition of stratigraphic Thomsen parameters will affect the accuracy of seismic data processing and interpretation.