# Influence of Error in Estimating Anisotropy Parameters on VTI Depth Imaging

Authors: S. Y. Moussavi Alashloo; D. P. Ghosh; Y. Bashir; W. Y. Wan Ismail
Journal: International Journal of Geophysics
 (2016)

---

## Abstract

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.


---

## Body

## 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:(1)
∂
t
∂
x
2
+
∂
t
∂
y
2
+
∂
t
∂
z
2
=
1
v
2
θ
,where t is the travel time, x,  y,  z are the Cartesian coordinates, and v
θ 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 qP 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:(2)
v
2
θ
≈
1
2
e
θ
+
1
2
e
2
θ
+
4
q
-
1
a
c
sin
2
⁡
θ
cos
2
⁡
θ
,where a
=
c
11,  c
=
c
33, where c
i
j
X are the density-normalized components of the elastic tensor, θ is the phase angle, and q and e
(
θ
) are the anellipticity coefficient and the elliptical component of the velocity, respectively, defined by (3)
q
=
1
+
2
δ
1
+
2
ε
=
1
1
+
2
η
,
(4)
e
θ
=
a
sin
2
⁡
θ
+
c
cos
2
⁡
θ
,where ε and δ are Thomsen’s parameters, and η
=
ε
-
δ
/
1
+
2
δ. 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 (5)
1
V
2
Θ
≈
1
+
2
Q
2
1
+
Q
E
Θ
+
1
2
1
+
Q
E
2
Θ
+
4
Q
2
-
1
A
C
sin
2
⁡
Θ
cos
2
⁡
Θ
,where Q
=
1
/
q,  A
=
1
/
a,  C
=
1
/
c,  Θ is the group angle, and E
(
Θ
) is the elliptical part:(6)
E
Θ
=
A
sin
2
⁡
Θ
+
C
cos
2
⁡
Θ
.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; andFarAway, which remains untouched ahead of the wavefront (Figure 1). In other words, the source positions at the beginning of the evaluation are consideredAlive. Given that they are initial points, their travel time is zero. All points that are one grid point away are taken asNarrowBand, and their travel times are computed analytically. All other grid points are marked asFarAway 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 theNarrowBand points.
(2)
Tag the point asAlive and remove it fromNarrowBand.
(3)
Check the neighbours of the minimum point that are notAlive. If any of them are categorized asFarAway, update them asNarrowBand. 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, (7)
∑
j
t
i
-
t
j
Δ
x
i
j
2
=
s
i
2
,should be solved for t
i which is the updated travel time value. t
j is the travel time at the neighbouring points, s
i is the slowness (s
=
1
/
v) at the point i, and Δ
x
i
j is the grid size in i
j direction.Figure 1
Schematic illustration of fast marching method.


### 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 by(8)
I
ξ
=
∫
Ω
ξ
W
ξ
,
m
,
h
D
t
=
t
D
ξ
,
m
,
h
,
m
,
h
d
m
d
h
.The image I
ξ, given in a two-dimensional space ξ
=
(
z
ξ
,
x
ξ
), can be determined from the integral on the data values D
(
t
,
m
,
h
) at the time t
D
(
ξ
,
m
,
h
) and weighted by an appropriate element W
ξ
,
m
,
h. The time factor t
D
(
ξ
,
m
,
h
) is defined as the total computed time while the waves travel from the source s to the image point I
(
ξ
) and propagate back to the receiver point r 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.


## 2.1. VTI Fast Marching Eikonal Solver

The anisotropic wavefield propagation under a high frequency assumption is defined by the eikonal equation:(1)
∂
t
∂
x
2
+
∂
t
∂
y
2
+
∂
t
∂
z
2
=
1
v
2
θ
,where t is the travel time, x,  y,  z are the Cartesian coordinates, and v
θ 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 qP 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:(2)
v
2
θ
≈
1
2
e
θ
+
1
2
e
2
θ
+
4
q
-
1
a
c
sin
2
⁡
θ
cos
2
⁡
θ
,where a
=
c
11,  c
=
c
33, where c
i
j
X are the density-normalized components of the elastic tensor, θ is the phase angle, and q and e
(
θ
) are the anellipticity coefficient and the elliptical component of the velocity, respectively, defined by (3)
q
=
1
+
2
δ
1
+
2
ε
=
1
1
+
2
η
,
(4)
e
θ
=
a
sin
2
⁡
θ
+
c
cos
2
⁡
θ
,where ε and δ are Thomsen’s parameters, and η
=
ε
-
δ
/
1
+
2
δ. 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 (5)
1
V
2
Θ
≈
1
+
2
Q
2
1
+
Q
E
Θ
+
1
2
1
+
Q
E
2
Θ
+
4
Q
2
-
1
A
C
sin
2
⁡
Θ
cos
2
⁡
Θ
,where Q
=
1
/
q,  A
=
1
/
a,  C
=
1
/
c,  Θ is the group angle, and E
(
Θ
) is the elliptical part:(6)
E
Θ
=
A
sin
2
⁡
Θ
+
C
cos
2
⁡
Θ
.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; andFarAway, which remains untouched ahead of the wavefront (Figure 1). In other words, the source positions at the beginning of the evaluation are consideredAlive. Given that they are initial points, their travel time is zero. All points that are one grid point away are taken asNarrowBand, and their travel times are computed analytically. All other grid points are marked asFarAway 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 theNarrowBand points.
(2)
Tag the point asAlive and remove it fromNarrowBand.
(3)
Check the neighbours of the minimum point that are notAlive. If any of them are categorized asFarAway, update them asNarrowBand. 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, (7)
∑
j
t
i
-
t
j
Δ
x
i
j
2
=
s
i
2
,should be solved for t
i which is the updated travel time value. t
j is the travel time at the neighbouring points, s
i is the slowness (s
=
1
/
v) at the point i, and Δ
x
i
j is the grid size in i
j direction.Figure 1
Schematic illustration of fast marching method.


## 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 by(8)
I
ξ
=
∫
Ω
ξ
W
ξ
,
m
,
h
D
t
=
t
D
ξ
,
m
,
h
,
m
,
h
d
m
d
h
.The image I
ξ, given in a two-dimensional space ξ
=
(
z
ξ
,
x
ξ
), can be determined from the integral on the data values D
(
t
,
m
,
h
) at the time t
D
(
ξ
,
m
,
h
) and weighted by an appropriate element W
ξ
,
m
,
h. The time factor t
D
(
ξ
,
m
,
h
) is defined as the total computed time while the waves travel from the source s to the image point I
(
ξ
) and propagate back to the receiver point r 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 (Figure2). 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 x
=
5 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 2
(a) Marmousi velocity model and (b) Marmousi synthetic shot gathers.
(a)
(b)Figure 3
Marmousi anisotropy model: (a) trueη model, (b) model with 30% reduced η value, (c) model with 30% added η value, and (d) model with 40% reduced η value.
(a)
(b)
(c)
(d)Figure 4
Fast marching travel time results of Marmousi model for isotropy (dashed green line), anisotropy with trueη value (light blue line), anisotropy with 30% error less than true η value (yellow line), anisotropy with 30% error more than true η value (magenta line), and anisotropy with 40% error less than true η value (dark blue line) for a source at x
=
5 km.Figure5 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.Figure 5
Comparison between (a) isotropic and (b) VTI Marmousi images. Superposition of reflectivity on images generated by (c) isotropic model, (d) true eta model, (e) eta model with 30% decreased value, (f) eta model with 30% increased value, and (g) eta model with 40% decreased value.
(a)
(b)
(c)
(d)
(e)
(f)
(g)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 Figure5(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.



---

*Source: 2848750-2016-05-15.xml*