Thin layers in sedimentary rocks lead to seismic anisotropy which makes the wave velocity dependent on the propagation angle. This aspect causes errors in seismic imaging such as mispositioning of migrated events if anisotropy is not accounted for. One of the challenging issues in seismic imaging is the estimation of anisotropy parameters which usually has error due to dependency on several elements such as sparse data acquisition and erroneous data with low signal-to-noise ratio. In this study, an isotropic and anelliptic VTI fast marching eikonal solvers are employed to obtain seismic travel times required for Kirchhoff depth migration algorithm. The algorithm solely uses compressional wave. Another objective is to study the influence of anisotropic errors on the imaging. Comparing the isotropic and VTI travel times demonstrates a considerable lateral difference of wavefronts. After Kirchhoff imaging with true anisotropy, as a reference, and with a model including error, results show that the VTI algorithm with error in anisotropic models produces images with minor mispositioning which is considerable for isotropic one specifically in deeper parts. Furthermore, over- or underestimating anisotropy parameters up to 30 percent are acceptable for imaging and beyond that cause considerable mispositioning.

1. Introduction

It is well-known that hydrocarbon reservoirs and overlying strata are commonly anisotropic [1, 2]. In reality, it is rare to have media with elliptical or weak anisotropy properties. However, anellipticity (deviation of wavefield from ellipse) has been commonly observed in the Earth’s subsurface, and it is a significant characteristic of elastic wave propagation [3, 4].

Another challenging issue in depth imaging is the computation of the travel time taken by a seismic wave from source to receiver. An efficient method to compute travel times is solving the eikonal equation by employing finite differences [5, 6]. Different techniques have been introduced to solve the eikonal equation, such as embedding methods, single-pass methods, sweeping methods, and iterative methods [7]. The main difference of these techniques is in how they cope with the complication of multivalued solutions and in finding solutions in the vicinity of cusps and discontinuities [8]. Anisotropy was initially added to an eikonal solver algorithm by Dellinger [9]. The embedding and iterative methods are both time consuming, particularly in heterogeneous and anisotropic conditions [7]. Fast sweeping methods are originally proposed for isotropic media [10]; however, a modification is executed to handle the anisotropic condition [11]. Single-pass or fast marching method (FMM) is another tool for computing travel times but is not generally applicable for anisotropic medium [5]. This algorithm has since been modified to work for anisotropy [12, 13].

In this study, a prestack depth migration algorithm is developed based on an anelliptic VTI compressional wave equation. Fomel’s anelliptic approximation [14] for both phase and group velocity of P-wave are employed to derive the eikonal equation. The fast marching finite difference approach is used as our eikonal solver since it is fast and stable for travel time computation. In anisotropic study, four anisotropic models are used: a true model which is exactly similar to the model employed for forward modelling, a model with values 30 percent less than the true model, a model with values 40 percent less than the true model, and a model with values 30 percent more than the true model. The calculated travel times are compared and employed in a standard Kirchhoff migration to obtain the image of the subsurface. The Marmousi model, as a complex model, is used to test the algorithm. Finally, we analyze both isotropic and VTI images qualitatively.

2. Methodology

We develop a new algorithm to incorporate anelliptic VTI travel times into prestack depth imaging. Our workflow for PSDM consists of the following: (1) travel time computation and (2) Kirchhoff depth migration. Step (1) provides travel times for the Kirchhoff migration. The algorithm is discussed in detail below.

2.1. VTI Fast Marching Eikonal Solver

The anisotropic wavefield propagation under a high frequency assumption is defined by the eikonal equation:where is the travel time, ,  ,   are the Cartesian coordinates, and is the VTI phase velocity. The FMM solves (1) by considering the fact that the direction of energy propagation follows the group velocity equation. This method is similar to a ray that is perpendicular to wavefronts defined by phase velocity. This ray is called the travel time gradient. A wave equation is needed as a kernel of fast marching algorithm. Fomel [14] enhanced the anelliptic P wave approximation proposed by Muir and Dellinger [15] through replacing the linear approximation with a nonlinear one. By using the shifted hyperbola approximation, he obtains the following equation for P wave phase velocity:where ,  , where are the density-normalized components of the elastic tensor, is the phase angle, and and are the anellipticity coefficient and the elliptical component of the velocity, respectively, defined by where and are Thomsen’s parameters, and . Similarly, for approximating the group velocity, the shifted hyperbola approach is applied on the Muir’s approximation to unlinearize the equation. The new group velocity approximation is where ,  ,  ,   is the group angle, and is the elliptical part:Approximations in (2) and (5) are used for ray tracing in locally homogeneous cells needed in this algorithm. An attentively selected order of travel time evaluation is the main advantage of the FMM. This method is an upwind method which means if a wave propagates from left to right, a difference scheme should be applied for reaching upwind to the left to collect information to construct the solution downwind to the right.

According to Fomel [16], although the algorithm follows a certain procedure, the grid points are divided into three classes, namely, Alive, points which are behind the wavefront and have been already computed; NarrowBand, points on the wavefront awaiting assessment; and FarAway, which remains untouched ahead of the wavefront (Figure 1). In other words, the source positions at the beginning of the evaluation are considered Alive. Given that they are initial points, their travel time is zero. All points that are one grid point away are taken as NarrowBand, and their travel times are computed analytically. All other grid points are marked as FarAway and have an “infinitely large” travel time value [5, 16]. The algorithm includes the following main steps:(1)Find the point with the minimum travel time among the NarrowBand points.(2)Tag the point as Alive and remove it from NarrowBand.(3)Check the neighbours of the minimum point that are not Alive. If any of them are categorized as FarAway, update them as NarrowBand. It means the wavefront is advanced, and the minimum point is behind it.(4)Update travel times for points on the NarrowBand (wavefront) by solving (1) numerically.(5)Repeat the loop until all points are behind the wavefront.For updating, one to three neighbour points are selected and their travel time values need to be less than the current value. After choosing the points, the quadratic equation, should be solved for which is the updated travel time value. is the travel time at the neighbouring points, is the slowness () at the point , and is the grid size in direction.

2.2. Kirchhoff Depth Migration

To image subsurface structures, which is normally recorded as geometrically unfocused, a migration procedure is needed to eliminate the influences on the wave during propagation in the subsurface. Migration technique moves seismic events to their true position, and it weakens diffractions. Since the prestack data have irregular spatial sampling, Kirchhoff migration is often the choice for prestack imaging [17]. Kirchhoff migration is dependent upon the solution of Kirchhoff integral to the wave equation and upon Green’s function theory. Its general form as an integral expression is defined byThe image , given in a two-dimensional space , can be determined from the integral on the data values at the time and weighted by an appropriate element . The time factor is defined as the total computed time while the waves travel from the source to the image point and propagate back to the receiver point at the surface. Overall, for Kirchhoff depth migration, we need two inputs which are seismic data and travel times. An anelliptic VTI fast marching approach is applied on the velocity models to compute travel times, and synthetic data are created via VTI finite difference forward modeling.

3. Numerical Examples

The results of the new VTI PSDM algorithm are demonstrated and discussed in this section. The efficiency of the above-described approach is tested for isotropic and VTI travel time computation on the 2D Marmousi model (Figure 2). The true model was generatedbased on the velocity model which varies between 0.03 and 0.1 (Figure 3(a)). Different percentages of error were introduced to the true model (Figures 3(b)3(d)). The VTI first arrival travel times are computed and compared with the isotropic travel times for sources at  km (Figure 4). These results indicate that the eikonal solver algorithm is able to cover all the desired area in the model even with high complexity such as normal faults and tilted blocks. Stability is another advantage of this method in which it is fast and robust to provide travel times [18]. Furthermore, by comparing the travel times, it is obvious that the anisotropic wavefronts laterally move faster than the isotropic ones; however, the isotropic and anisotropic wavefronts propagate in a same speed vertically. This difference is due to the parameter that affects the wave propagation mainly laterally. Waheed et al. [6] demonstrated that the maximum influence of on the propagation of the wave is in the direction orthogonal to the symmetry axis. The errors introduced to the anisotropy parameter affected the travel times which below 30 percent is not too much but above 30 percent causes a noticeable difference. The gap between isotropic and VTI warfronts increases while the wavefronts propagate away from the source. To know how these variations influence imaging, the Marmousi synthetic data was migrated by using these travel times.

Figure 5 illustrates the images of VTI and isotropic PSDM for Marmousi model as well as the overlay of the images with the reflectivity model which was derived from the velocity model and its calculated acoustic impedance. Comparison of Figures 5(a) and 5(b) demonstrates several differences. Circles, number 1, show misshaping of pinchout and a wrong dip for the reflector in the isotropic image. Also, the spacing between the isotropic interfaces is more than the anisotropic ones.

The resolution is another effect where, in the anisotropic image, layers appear to be more clear and sharp, but, in the isotropic image, they are poorly imaged.

In rectangles, number 2, more fractures are detected by VTI imaging, and different positioning is obvious. Number 3 again indicates changing in the reflector’s dip where the isotropic layer has steep tilt compared to anisotropic layer. Comparing the images with the model confirms that the VTI result is much more accurate, and it is better matched.

Considering Figure 5(d) as a reference, which is imaged by true model and by comparing it with the images from the modified models, it can be concluded that the images difference for those with below 30 percent errors is not remarkable. In contrast, for those above 30 percent errors, whether higher or lower than true value, there is more mispositioning in the deeper sections. Hence, if the estimated anisotropy parameter has error lower than 30 percent, the result of our VTI algorithm is reliable for interpretation and other applications. However, by increasing error in model, mispositioning becomes more apparent that geophysicists need to consider it.

4. Conclusions

A prestack VTI depth imaging algorithm was designed and applied for the isotropic and anisotropic media. A VTI fast marching eikonal solver based on anelliptic was utilized to compute the seismic travel times. The Marmousi velocity model and its seismic dataset are used for imaging. After Kirchhoff imaging, results showed that the VTI imaging algorithm generates images with higher resolution than the isotropic ones specifically in deeper parts, and the subsurface features are imaged correctly regarding their pattern and location. Moreover, study on error in anisotropy estimation demonstrated that up to 30 percent error can be tolerated in imaging but beyond that leads to considerable mispositioning which may affect the interpretation of the subsurface.

Competing Interests

The authors declare that they have no competing interests.


The authors gratefully acknowledge members of Centre of Seismic Imaging (CSI) at UTP for their helpful discussions. This work is funded by PETRONAS. Seismic package of Madagascar was used in developing the algorithm.