Pesquisa cientficiaa

Prévia do material em texto

FEDERAL UNIVERSITY OF PARÁ
GEOSCIENCE INSTITUTE
GRADUATE PROGRAM IN GEOPHYSICS
MASTER’S DISSERTATION
Wavefield Decomposition Based RTM using Low Rank
Approximation in TTI media
ANDREY MARCOS SOUZA DA SILVA DE LIMA
Belém-PA
2018
ANDREY MARCOS SOUZA DA SILVA DE LIMA
Wavefield Decomposition Based RTM using Low Rank
Approximation in TTI media
Dissertation presented to the Graduate Program in Geo-
physics from Federal University of Pará.
Concentration Area: Modeling and Inversion of Geo-
physical data
Research Line: Imaging of Geophysical Data
Advisor: Prof. Dr. Daniel Leal Macedo
Belém-PA
2018
Dados Internacionais de Catalogação na Publicação (CIP) de acordo com ISBD
Sistema de Bibliotecas da Universidade Federal do Pará
Gerada automaticamente pelo módulo Ficat, mediante os dados fornecidos pelo(a) autor(a)
L732w Lima, Andrey Marcos Souza da Silva de.
Wavefield decomposition based RTM using low rank approximation in TTI media / Andrey Marcos
Souza da Silva de Lima. — 2018.
47 f. : il. color.
Orientador(a): Prof. Dr. Daniel Leal Macedo
Dissertação (Mestrado) - Programa de Pós-Graduação Geofísica, Instituto de Geociências, Universidade
Federal do Pará, Belém, 2018.
1. Modelagem Sísmica. 2. Anisotropia. 3. Decomposição de Campos. 4. Baixo posto. 5. Imageamento
Sísmico. I. Título.
CDD 622.1592
I dedicate this work to my family who gave
me the necessary support to arrive where I
am now. I also dedicate this work to my wife
who is an incredible friend.
ACKNOWLEDGMENTS
I thank all my professors at UFPA specially Dr. Jessé and Dr. Daniel who helped me
a lot during both my undergraduate and graduate courses. I also thank the scholarships
and sponsors that I had throughout my academic life so far (PET-MEC, Capes-Science
without Borders at University of Oklahoma, SBGF-Scientific Initiation, CNPQ-Master’s
degree and Petrobras for funding the UFPA group).
RESUMO
Apresentamos uma metodologia combinando a aproximação de baixo posto para simular
a propagação de ondas em meios TTI com RTM usando decomposição de campos ascen-
dente e descendente para obter imagens com menos ruídos de baixa frequência em áreas
complexas. Comparamos três resultados principais: 1-aproximação de baixo posto sem
decomposição do campo de onda, 2-com decomposição do campo de onda e 3-uma imple-
mentação de diferenças finitas. Uma vez que o RTM usando decomposição precisa de ape-
nas duas transformadas de Hilbert extras, ele pode ser implementado facilmente usando
o algoritmo de baixo posto. Os resultados mostraram que o RTM usando decomposição
do campo de onda é capaz de remover artefatos e imagear ambientes anisotrópicos.
Palavras-chaves: Métodos Sísmicos. Modelagem Sísmica. Imageamento Sísmico.
Anisotropia. Decomposição de Campos. Baixo posto.
ABSTRACT
We have presented a methodology combining the low rank approximation to simulate the
wave propagation in TTI media with the wavefield decomposition based RTM to obtain
images with less low frequency noise in complex areas. We have compared three main
results: 1-low rank without wavefield decomposition, 2-with wavefield decomposition and
3-a finite differences implementation. Since the RTM using decomposition needs only two
extra Hilbert transforms, it can be easily implemented using the low rank algorithm. The
results showed that RTM using wavefield decomposition is capable of removing artifacts
and imaging anisotropic environments.
Keywords: Seismic Methods. Seismic Modeling. Seismic Imaging. Anisotropy.
Wavefield Decomposition. Lowrank.
LIST OF FIGURES
2.1 Illustration showing some anisotropic models and the parameters they rely
on. The more complex is the media, the more parameters are necessary to
extrapolate the wavefield and consequently to image the subsurface. Figure
from Martinez (2018). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1 Illustration of the extrapolation algorithm of the wavefield in anisotropic
media. Firstly, the wavefield at a time t (inner blue curve) is decomposed
in plane waves (red tangent lines). Secondly, each plane wave is propagated
during the interval ∆t with its respective phase velocity V (~x,~k) for each
point in the wavefront. Finally, the wavefield is synthesized (reconstructed)
at time t+∆t (outer blue curve). Figure from Baltazar and Carvalho (2018). 14
4.1 The low rank algorithm takes the extrapolation operator and decompose
into three smaller matrices. . . . . . . . . . . . . . . . . . . . . . . . . . . 17
5.1 Representation and example of the false image problem. The red arrows
indicate how these artifacts can contaminate the image if not removed.
Figures modified from Fei et al. (2015) . . . . . . . . . . . . . . . . . . . . 22
5.2 Representation of the low-frequency artifacts. The top of the salt body
presents smiles formed by low-frequencies waves cross-correlations. Figure
from Kaelin and Carvajal (2011). . . . . . . . . . . . . . . . . . . . . . . . 23
6.1 In (a) the back-propagated wavefield without decomposition, in (b) the
back-propagated wavefield decomposed, in (c) the forward wavefield with-
out decomposition and in (d) the forward wavefield decomposed. Both x
and z axis are in meters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6.2 All three results and the velocity model in (a), (b) is the RTM without the
decomposition, (c) is the RTM using finite differences and (d) the RTM
using the decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.3 In (a) is showed the RTM using the finite differences implementation and
in (b) the RTM using low rank without decomposition. . . . . . . . . . . . 27
6.4 In (a) is showed the RTM using low rank without decomposition and in
(b) the RTM using low rank with decomposition. . . . . . . . . . . . . . . 27
6.5 In (a) is showed the RTM using the decomposed wavefield and in (b) the
RTM using the finite differences implementation. . . . . . . . . . . . . . . 28
6.6 In (a) the back-propagated wavefield without decomposition, in (b) the
back-propagated wavefield decomposed, in (c) the forward wavefield with-
out decomposition and in (d) the forward wavefield decomposed . . . . . . 29
6.7 In (a) the velocity model, in (b) the epsilon model, in (c) the arbitrary
angles of symmetry model and in (d) the delta model. As we can see
from item (c) we did not select an area with strong variations on angle of
symmetry. We focused our analysis on the top of the salt to see if could
remove the low frequency artifacts. . . . . . . . . . . . . . . . . . . . . . . 30
6.8 In (a) the velocity model, in (b) the RTM not decomposed, in (c) the finite
differences implementation and in (d) the RTM decomposed. . . . . . . . . 31
6.9 In (a) the RTM using finite differences and in (b) the RTM without de-
composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.10 In (a) the RTM without decomposition and in (b) the RTM using decom-
position. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.11 In (a) RTM using decomposition and in (b) RTM using finite differences. . 33
CONTENTS
1 INTRODUCTION 1
2 BASIC THEORY IN ANISOTROPY 3
2.1 EQUATION OF MOTION, WAVE EQUATION AND HOOKE’S LAW . 3
2.2 PLANE WAVE SOLUTION . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 ANISOTROPIC SYSTEMS . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3 WAVEFIELD EXTRAPOLATION 9
3.1 PSEUDO ANALYTIC METHOD . . . . . . . . . . . . . . . . . . . . . . 9
3.2 WAVEFIELD EXTRAPOLATION OPERATOR . . . . . . . . . . . . . . 10
3.3 PHASE VELOCITY IN ANISOTROPIC MEDIA . . . . . . . . . . . . . 14
4 LOW RANK APPROXIMATION 16
4.1 LOW RANK MATRIX, W 1 AND W 2 . . . . . . . . . . . . . . . . . . . . 18
4.2 KACZMARZ ALGORITHM . . . . . . . . . . . . . . . . . . . . . . . . . 19
5 WAVEFIELD DECOMPOSITION BASED REVERSE TIMEMIGRA-
TION 21
5.1 RTM PROBLEMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2 WAVEFIELD DECOMPOSITION USING COMPLEX TRACES . . . . . 23
6 RESULTS25
6.1 MARMOUSI MODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6.1.1 Finite Differences x Low Rank Not Decomposed . . . . . . 27
6.1.2 Low Rank Not Decomposed x Decomposed . . . . . . . . . 27
6.1.3 Low Rank Decomposed x Finite Differences . . . . . . . . . 28
6.2 BP TTI MODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
6.2.1 Finite Differences x Low Rank Not Decomposed . . . . . . 32
6.2.2 Low Rank Not Decomposed x Decomposed . . . . . . . . . 32
6.2.3 Low Rank Decomposed x Finite Differences . . . . . . . . . 33
7 CONCLUSIONS 34
Bibliography 35
1 INTRODUCTION
Seismic imaging is a crucial part of any exploration project in the industry. The more
complex the geology, the bigger are the challenges to image the subsurface. Thanks to the
increasing computational power, Reverse Time Migration or RTM (Baysal et al., 1983;
Whitmore, 1983; McMechan, 1983) has become the standard tool to image these complex
environments.
The conventional imaging condition in RTM is the zero-lag crosscorrelation of the
forward-propagated source wavefield and the backward-propagated receiver wavefield.
This imaging condition provides the correct kinematics and is extremely simple to imple-
ment (Claerbout, 1985). However, undesired crosscorrelations along the wavepath produce
high amplitude, low-frequency noise and false images (Liu et al., 2011). Theses artifacts
are due to the strong velocity contrasts. When performing the Hilbert transform to the
source wavefield and decomposing the upgoing and downgoing wavefields, one can apply
the imaging condition to the appropriate combinations of the wavefield components (Liu
et al., 2007).
Wavefield extrapolation is the core of seismic modeling and consequently of seismic
imaging (RTM). There are many approaches to seismic modeling; we classify them into
three main categories: direct methods (finite differences-FD, pseudospectral-PS, and fi-
nite elements-FE), integral-equation methods (based on the Green’s function), and ray-
tracing methods (based on the eikonal and transport equations). Another way of see-
ing the wavefield extrapolation in time is to numerically approximate the mixed-domain
space-wavenumber operator (Wards et al., 2008). This work will discuss the low rank
approximation from Fomel et al. (2013), which basically selects a small set of representa-
tive spatial locations and a small set of representative wavenumbers of this operator, also
called pseudo-laplacian operator. This methodology allows the modeling of the quasi-P
(q-P) waves in anisotropic media with arbitrary symmetry axis and planes of symmetry
(Tsvankin, 2012). The resulting algorithm can be used for RTM in Tilted Transverse
Isotropic (TTI) media.
The application of seismic modeling in anisotropic media to generate seismic images
has some drawbacks. The correct treatment of the wavefields in anisotropic media needs
the solution of the elastodynamics equations. However, we conventionally record only
q-P waves and it is very computationally expensive to estimate the elastic model. These
two problems motivated the development of scalar wave equations that describe approx-
imately only the elastic q-P mode, called pseudo-acoustic equations (Alkhalifah, 1998).
Unfortunately, the solutions of these equations also contain spurious modes (pseudo-S)
due to the mathematical approximations. A lot of effort has been made to mitigate these
spurious modes. An overview of these techniques are found in Yan and Liu (2016) and
1
2
Schleicher and Costa (2016).
Etgen and Brandsberg-Dahl (2009) proposed the pseudo-analytical method to obtain
the approximate solutions to these pseudo-acoustic equations without the presence of these
spurious events in Tilted Transverse Isotropic media (TTI). The lowrank approximation
extends the pseudo-analytic method with optimally selected velocitites and weights to
any class of anisotropy and heterogeneity. The implementation of the lowrank algorithm
through least squares can present high storage demand and it can be very computation-
ally inefficient to applications in 3D seismic modeling. For these reasons, we will follow
the work of Baltazar and Carvalho (2018). They have used the Kaczmarz algorithm
(Kaczmarz, 1993) to solve the linear system of equations that appears after the pivoted
QR decompositions (Golub and Van Loan, 2012), used to find the best representatives
space and wavenumber locations. This methodology, combined with the separation of the
forward- and backward-propagated wavefields, was efficiently applied for imaging in the
pseudo-acoustic approximation both the Marmousi model (Versteeg, 2005) and the BP
model (Billette and Brandsberg-Dahl, 2005).
This work presents a methodology to separate the upgoing and downgoing wavefields
using the Hilbert transform (Zheng et al., 2018), combined with the lowrank approxima-
tion (Fomel et al., 2013), to produce better images when applying RTM in anisotropic
media. We have compared three main results: 1-low rank without wavefield decompo-
sition, 2-with wavefield decomposition and 3-a finite differences implementation. Since
the RTM using decomposition needs only two extra Hilbert transforms, it can be easily
implemented using the low rank algorithm. The results showed that RTM using wavefield
decomposition is capable of removing artifacts and imaging anisotropic settings.
2 BASIC THEORY IN ANISOTROPY
2.1 EQUATION OF MOTION, WAVE EQUATION AND HOOKE’S LAW
The Newton’s law or the equation of movement including an external force is
ρ
∂2ui
∂t2
− ∂σij
∂xj
= fi , (2.1)
where ρ is the density, ui = (u1, u2, u3) is the displacement vector, fi = (f1, f2, f3) is the
external body force per unit volume, t is time, xj are the Cartesian coordinates, and σij
is the stress tensor. Here we use the Einstein notation.
For a medium with certain density and certain spatial distribution of body forces,
fi, the equation 2.1 has two unknowns: the displacement field ui and the stress
tensor σij. Hence, the wave equation cannot be solved for the displacement unless
we have a constitutive relation between stress and deformation (or between stress and
displacement). In the specific case of small deformations (in the same scale of seismic wave
perturbations), the relationship between stress and deformation is linear and described
by the generalized Hooke’s law:
σij = cijkl ekl (2.2)
where cijkl is a fourth order tensor called elastic stiffness tensor (or elastic parameters
tensor) which reveals how each point of the medium reacts when deformed, and ekl is the
deformation tensor (also called strain tensor) defined as (with only the linear terms):
ekl =
1
2
(
∂uk
∂xl
+
∂ul
∂xk
)
(2.3)
Considering the symmetry of the stiffness tensor, cijkl = cijlk we can find the following
expression
cijklekl = cijklekl =
1
2
(
∂uk
∂xl
+
∂ul
∂xk
)
=
1
2
cijkl
∂uk
∂xl
+
1
2
cijkl
∂ul
∂xk
=
1
2
cijkl
∂uk
∂xl
+
1
2
cijlk
∂ul
∂xk
= cijkl
∂uk
∂xl
(2.4)
If we substitute 2.2 and 2.4 in 2.1 we will have
ρ
∂2ui
∂t2
− ∂(cijkl ekl)
∂xj
= fi (2.5)
ρ
∂2ui
∂t2
−
∂(cijkl
∂uk
∂xl
)
∂xj
= fi (2.6)
3
4
and considering that cijkl does not change with the position or it is very little dependent,
it can leave the spatial derivative
ρ
∂2ui
∂t2
− cijkl
∂(∂uk
∂xl
)
∂xj
= fi (2.7)
ρ
∂2ui
∂t2
− cijkl
∂2uk
∂xj∂xl
= fi (2.8)
By changing cijkl we find the solution for equation 2.8 specific for an elastic media. Cijkl
is a fourth order tensor, which in theory would generate a tensor with 81 independent
parameters. Due to the symmetry of the tensors σij = σji, ekl = elk and because the wave
propagation is considered an adiabatic process (which is translated in cijkl = cklij), the
stiffness tensor has 21 independent parameters which characterizes the most complex
elastic media (triclinic). Simplifying hypothesis in cijkl by c(trc), c(mnc), c(ort), c(vti), c(hti)
and c(iso) we will have the wave equation for the triclinic, monoclinic, orthorhombic, VTI,
HTI and isotropic media respectively. Theseanisotropic systems will be examined in the
next section.
2.2 PLANE WAVE SOLUTION
In order to give an analytic description of the plane waves in an anisotropic medium,
we evaluate equation 2.1 in the homogenous case by setting fi = 0
ρ
∂2ui
∂t2
− ∂σij
∂xj
= 0 (2.9)
Physically the homogeneous wave equations describe a medium without sources of elastic
energy. As a trivial solution for equation 2.9 we use the harmonic plane wave represented
by
uk = Uk e
iω(
njxj
V
−t) = Uk e
i(~k·~x−ωt) (2.10)
Where Uk are the vector components of the polarization vector ~U , ω is the angular
frequency, V is the propagation velocity of each plane wave (called phase velocity), and
nj is the unit vector orthogonal to the wavefront plane (the wavefront satisfy nj xj−V t =
constant). Another useful quantity in anisotropy is the slowness vector ~p = ~n/V where,
~p.ω = ~k. Substituting 2.10 and the Hooke’s law 2.2 in wave equation 2.9, we arrive at the
so called Christoffel equation to the phase velocity V and polarization vector ~U :Γ11 − ρV
2 Γ12 Γ13
Γ21 Γ22 − ρV 2 Γ23
Γ31 Γ32 Γ33 − ρV 2

U1U2
U3
 = 0 (2.11)
5
Where Γik = cijkl njnl is the Christoffel matrix that depends on the properties of the
medium and the direction that the wave propagates. Since the stiffness tensor is symmetric
(cijkl = cklij), the Christoffel matrix is also symmetric Γik = Γki. We can remove the
density ρ from Christoffel equation if we define [aijkl = cijkl/ρ] which is the normalized
stiffness tensor. If we use the Kronecker delta δik (where δik = 1 if i ≡ k, and δik = 0 if
i 6≡ k) the equation 2.11 can be rewritten in a more compact form
[Γik − ρV 2δik]Uk = 0 (2.12)
The Christoffel equation 2.12 describes an eigenvalue (ρV 2 ) - eigenvector (U) problem.
Γ is a symmetric and positive defined matrix. In this way the eigenvalues can be obtained
by calculating the determinant of the expression in the bracket since Uk cannot be zero
det[Γik − ρV 2δik] = 0 (2.13)
When we calculate the determinant of this matrix, we arrive at a third-degree polynomial
called the characteristic polynomial. For any direction ~n in an anisotropic medium the
Christoffel equation provides 3 possible velocity values of the phase velocity V correspond-
ing to the velocities of the wave P, SV and SH. The isotropic case would be a degenerated
case (with too much symmetry) of the anisotropic case where the two S-wave velocities
are the same (Tsvankin, 2012).
If we plot the phase velocity for a given wave mode as a vector-radius in all directions of
propagation ~n, we define the surface of the phase velocity. If we plot the inverse value 1/V
results in the slowness surface, also called Hamiltonian, where its shape is directly related
to the wave fronts of a point source in the presence of singularities in the S waves. The
direction of the radius (defined by the group velocity vector) generally deviates from the
direction of the phase velocity ~n, but the direction of the radius is always perpendicular to
the surface of slowness. If the medium is homogeneous and isotropic, the phase velocity
and slowness surfaces are spherical.
After the eigenvalues have been determined, their eigenvectors U can be found using
the equations 2.12. Although the magnitude of the eigenvectors is undefined, their ori-
entation determines the polarization of the plane waves 2.10 propagating in the direction
~n. The polarization vectors of the plane waves in isotropic media are either parallel (P
waves) or are perpendicular (S waves) to the slowness vector. However, when we have
anisotropy the polarization is governed by both the orientation of the vector ~n and by the
stiffness tensor that thus changes the shape of the Christoffel equation.
Since Γik is a real and symmetric matrix, the vectors of polarization of the 3 wave
modes eigenvalues) are orthogonal to each other, but not necessarily parallel or perpen-
dicular to ~n. Thus, apart from specific directions of propagation, there are no purely
longitudinal and shear waves in an anisotropic medium (Tsvankin, 2012). Thus, in the
6
anisotropic wave theory we refer to the faster mode as quasi-P or q-P wave, and for the
slower modes quasi-SV and q-SH or quasi-S1 and q-S2.
2.3 ANISOTROPIC SYSTEMS
Depending on the medium, the wave equation 2.6 will be solved using a different elastic
stiffness tensor cijkl. Below are the main anisotropic media which are classified based on
the symmetries of the stiffness tensors taken from (Tsvankin, 2012). Here we use the Voigt
notation where ij → α, 11→ 1; 22→ 2; 33→ 3; 23, 32→ 4; 13, 31→ 5 and 12, 11→ 6.
triclinic: 21 independent coefficients
c(trc) =

c11 c12 c13 c14 c15 c16
c12 c22 c23 c24 c25 c26
c13 c23 c33 c34 c35 c36
c14 c24 c34 c44 c45 c46
c15 c25 c35 c45 c55 c56
c16 c26 c36 c46 c56 c66

(2.14)
monoclinic: 13 independent coefficients
c(mnc) =

c11 c12 c13 0 0 c16
c12 c22 c23 0 0 c26
c13 c23 c33 0 0 c36
0 0 0 c44 c45 0
0 0 0 c45 c55 0
c16 c26 c36 0 0 c66

(2.15)
orthorhombic: 9 independent coefficients
c(ort) =

c11 c12 c13 0 0 0
c12 c22 c23 0 0 0
c13 c23 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66

(2.16)
7
Tranverse Isotropic - VTI: 5 independent coefficients
c(vti) =

c11 c11 − 2c66 c13 0 0 0
c11 − 2c66 c11 c13 0 0 0
c13 c13 c33 0 0 0
0 0 0 c55 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66

(2.17)
Transverse Isotropic - HTI: 5 independent coefficients
c(hti) =

c11 c13 c13 0 0 0
c13 c33 c33 − 2c44 0 0 0
c13 c33 − 2c44 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c55

(2.18)
isotropic: 2 independent coefficients
c(iso) =

λ− 2µ λ λ 0 0 0
λ λ− 2µ λ 0 0 0
λ λ λ− 2µ 0 0 0
0 0 0 µ 0 0
0 0 0 0 µ 0
0 0 0 0 0 µ

(2.19)
Where λ and µ are the Lamé’s parameters.
8
The anisotropic media used in this work was the TTI media, which is basically a VTI
media dependent on the angle propagation direction, not only vertical direction. An elastic
VTI medium is dependent on 5 parameters plus the angle of symmetry which is constant,
and in 3D it is dependent on azimuth (φ). After the pseudo acoustic approximation
(Alkhalifah, 1998), it becomes dependent only on 4 parameters: Vp, δ, ε, θ and φ. In
figure 2.1 we can see how a TTI media looks like and the parameters that this media relies
on. The more complex is the media, the more parameters are necessary to extrapolate
the wavefield and consequently to image the subsurface.
Figure 2.1: Illustration showing some anisotropic models and the parameters they rely
on. The more complex is the media, the more parameters are necessary to extrapolate
the wavefield and consequently to image the subsurface. Figure from Martinez (2018).
3 WAVEFIELD EXTRAPOLATION
The first step to RTM imaging using low rank modeling is to extrapolate the wavefront,
which basically assumes that we know the wavefield (acoustic, elastic, pseudo-acoustic,
etc) at a specific time. Then we decompose it in planes waves (going from space domain to
wavenumber domain). After the decomposition we obtain the spectrum of the wavefield
which represents the amplitude of each plane wave. Each plane wave is multiplied by an
extrapolation operator in order to propagate them. Finally we return to space domain and
reconstruct the wavefront in a new space location and with the interval evolution time.
We are going to present below the pseudo analytic methodology introduced by Etgen and
Brandsberg-Dahl (2009). The following section is devoted to calculate the extrapolation
operator, and in the end we will observe that it depends exclusively on the phase velocity.
The lowrank approach on the next section can be seen as an extension of this pseudo
analytic method.
3.1 PSEUDO ANALYTIC METHOD
The constant density acoustic wave equation without source is
∂2p(~x, t)
∂t2
= V (~x)2∇2p(~x, t) , (3.1)
where p(~x, t) is the pressure field, dependent on space ~x = (x1, x2, x3) and time t, V (~x) is
the wavepropagation velocity in the medium and ∇2 is the laplacian operator. If we use
the finite-difference approximation in the left part of equation 3.1 we have
p(~x, t+ ∆t)− 2p(~x, t) + p(~x, t−∆t)
∆t2
= V (~x)2∇2p(~x, t) (3.2)
Applying the space Fourier transform in all terms and assuming that the velocity is
constant and that the laplacian operator ∇2 is substituted by a pseudo laplacian operator
F (~k) we have
P (~k, t)e−iω∆t − 2P (~k, t) + P (~k, t)e+iω∆t = ∆t2V (~x)2F (~k)P (~k, t) (3.3)
and summing the similar terms we have
P (~k, t)
(
e−iω∆t + e+iω∆t − 2
)
= ∆t2V (~x)2F (~k)P (~k, t) (3.4)
which can be rewritten using Euler’s identity, cos(θ) = e−iθ+e+iθ
2
,
9
10
P (~k, t) (2cos(ω∆t)− 2) = ∆t2V (~x)2F (~k)P (~k, t) (3.5)
If we isolate the F (~k) term we have
F (~k) =
[2cos(ω∆t)− 2]
∆t2V (~x)2
(3.6)
where F (~k) is the pseudo laplacian operator. If we substitute the phase velocity,
Vphase(~x,~k) =
ω
~k
(3.7)
we have the expression for the pseudo laplacian operator.
F (~k) =
[
2cos(Vphase(~x,~k) · |~k|∆t)− 2
]
∆t2V (~x)2
(3.8)
The low rank approximation can be viewed as an extension of this pseudo analytic method
with optimally selected velocities and weights. However, it selects the best wavenumbers
~k and propagation directions ~x through the pivoted QR decomposition.
3.2 WAVEFIELD EXTRAPOLATION OPERATOR
Assuming that the pseudo-acoustic wavefield is known at time t0 in all propagation
domain, p(~x, t0), Fomel et al. (2013) proposed the following representation to the wavefield
extrapolation operator
p(~x, t0 + t) =
∫
P (~k, t0) e
iφ(~x,~k,t) d3~k (3.9)
where φ(~x,~k, t) is the phase operator defined in the mixed domain (~x−~k, space-wavenumber),
and it determines the evolution of the wave in the interval t while
P (~k, t0) =
1
(2π)3
∫
p(~x, t0) e
−i~k·~x d3~x (3.10)
indicates the direct Fourier transform in space of p(~x, t0) and
p(~x, t0) =
∫
P (~k, t0) e
i~k·~x d3~x (3.11)
defines the inverse Fourier transform in space of p(~x, t0).
If we substitute representation 3.9 in the acoustic wave equation,
∂2p
∂t2
− V 2∇2p = 0, (3.12)
11
we will need to evaluate the two expressions
∂2p
∂t2
= −
∫
P (~k, t0)
[(
∂φ
∂t
)2
+ i
∂2φ
∂t2
]
eiφ(~x,
~k,t)d3~k (3.13)
and
∇2p = −
∫
P (~k, t0)
[
∇φ · ∇φ+ i∇2φ
]
eiφ(~x,
~k,t) d3~k . (3.14)
The following step is to use the high frequency approximation from Cerveny (2005) to
the phase function that assumes ∣∣∣∣∂φ∂t
∣∣∣∣2 >> ∣∣∣∣∂2φ∂t2
∣∣∣∣ , (3.15)
|∇φ · ∇φ| >>
∣∣∇2φ∣∣ . (3.16)
Using these approximations we can write
∂2p
∂t2
≈ −
∫
P (~k, t0)
[(
∂φ
∂t
)2]
eiφ(~x,
~k,t)d3~k (3.17)
and
∇2p ≈ −
∫
P (~k, t0) (∇φ · ∇φ) eiφ(~x,
~k,t)d3~k. (3.18)
Consequently,
∂2p
∂t2
− V 2∇2p ≈ −
∫
P (~k, t0)
[(
∂φ
∂t
)2
− V 2∇φ · ∇φ
]
d3~k = 0 , (3.19)
we then obtain a differential equation for the phase function φ(~x,~k, t):(
∂φ
∂t
)2
− V 2(~x,~k)∇φ · ∇φ = 0. (3.20)
We can see from this equation that the phase velocity depends on the space and
wavenumber, V (~x,~k), which is essencial to the extrapolation operator p(~x, t0 + t) =∫
P (~k, t0) e
iφ(~x,~k,t) d3~k simulate the wave propagation in anisotropic models (Etgen and
Brandsberg-Dahl, 2009).
Assuming that the extrapolation interval t between two frames is small, φ(~x,~k, t) can
be approximated by its Taylor’s series,
φ(~x,~k, t) = φ(~x,~k, 0) +
∂
∂t
φ(~x,~k, 0)t+
1
2
∂2
∂t2
φ(~x,~k, 0)t2 +O(t3) . (3.21)
12
Looking at equation above and p(~x, t0) =
∫
P (~k, t0) e
i~x·~k d3~x, we can observe that when
the interval t = 0 we obtain
φ(~x,~k, 0) = ~k · ~x . (3.22)
If we take the gradient (∇) and the first order time derivative of 3.21 we obtain
∇φ(~x,~k, t) = ~k +∇∂φ
∂t
(~x,~k, 0) t+O(t2) (3.23)
and
∂φ
∂t
(~x,~k, t) =
∂
∂t
φ(~x,~k, 0) +
∂2
∂t2
φ(~x,~k, 0) t+O(t2) . (3.24)
The two expressions above will be used to approximate equation 3.20. Consequently,(
∂φ
∂t
(~x,~k, t)
)2
=
(
∂
∂t
φ(~x,~k, 0) +
∂2
∂t2
φ(~x,~k, 0) t
)2
=
(
∂
∂t
φ(~x,~k, 0)
)2
+ 2
(
∂
∂t
φ(~x,~k, 0)
)(
∂2
∂t2
φ(~x,~k, 0) t
)
+O(t2)
(3.25)
and
∇φ(~x,~k, t) · ∇φ(~x,~k, t) =
(
~k +∇∂φ(~x,
~k, 0)
∂t
t
)
·
(
~k +∇∂φ(~x,
~k, 0)
∂t
t
)
= ~k · ~k + ~k · ∇∂φ(~x,
~k, 0)
∂t
t+ ~k · ∇∂φ(~x,
~k, 0)
∂t
t+O(t2)
= ~k · ~k + 2~k · ∇∂φ(~x,
~k, 0)
∂t
t
(3.26)
If we substitute 3.25 and 3.26 in equation 3.20, we will have the following expression,(
∂
∂t
φ(~x,~k, 0)
)2
+ 2
(
∂
∂t
φ(~x,~k, 0)
)(
∂2
∂t2
φ(~x,~k, 0) t
)
+
(
∂2
∂t2
φ(~x,~k, 0) t
)2
− (V (~x,~k))2
[
~k · ~k + 2~k · ∇∂φ(~x,
~k, 0)
∂t
t
]
= 0
(3.27)
Matching the zero order terms (independent of t) we have(
∂
∂t
φ(~x,~k, 0)
)2
− (V (~x,~k))2 ~k · ~k = 0 (3.28)
∂
∂t
φ(~x,~k, 0) = ±V (~x,~k)~k (3.29)
13
Matching the first order terms (dependent of t) we have
2
(
∂
∂t
φ(~x,~k, 0)
)(
∂2
∂t2
φ(~x,~k, 0) t
)
− (V (~x,~k))2 2~k · ∇∂φ(~x,
~k, 0)
∂t
t = 0 (3.30)
±V (~x,~k)|~k| ∂
2
∂t2
φ(~x,~k, 0) = (V (~x,~k))2 ~k · ∇ ± V (~x,~k)~k (3.31)
∂2
∂t2
φ(~x,~k, 0) = ±V (~x,~k)∇V (~x,~k) · ~k (3.32)
Now, with both equation 3.29 and equation 3.32 we can construct the phase operator
φ(~x,~k, t) to obtain 3.33
φ(~x,~k, t) = ~k · ~x+ V (~x,~k)~k t+ 1
2
V (~x,~k)∇V (~x,~k) · ~k t2 +O(t3) (3.33)
If we substitute this phase operator φ(~x,~k, t) in the extrapolation operator, p(~x, t0 + t) =∫
P (~k, t0) e
iφ(~x,~k,t) d3~k, we obtain
p(~x, t0 + t) =
∫
P (~k, t0)e
i(~k·~x+V ~k t+ 12V∇V ·~k t2) d3~k (3.34)
which is the extrapolation operator using the second order approximation of phase’s ex-
pansion. If we use the first order approximation, we don’t need to calculate the ∇V and
obtain
p(~x, t0 + t) =
∫
P (~k, t0)e
+i(V ~k t)ei
~k·~x d3~k. (3.35)
This equation is called one step extrapolator. Another way of representing this equation
is to change the signal of t
p(~x, t0 − t) =
∫
P (~k, t0)e
−i(V ~k t)ei
~k·~x d3~k. (3.36)
and sum both equations to create the two step extrapolator
p(~x, t0 + t) + p(~x, t0 − t) = 2
∫
P (~k, t0)e
i~k·~xcos
([
V (~x,~k)~k
]
t
)
d~k (3.37)
Equation 3.37 is the one implemented on this work. As we can see, the extrapolation
operator depends on the phase velocity only. If we find an expression for this V (~x,~k) we
then have the velocity for each point of the model, ~x, to all propagation directions, ~k.
Figure 3.1 represents all the steps of the wavefield extrapolation algorithm in the mixed
domain ~x − ~k. Initially, the wavefield at a time t is decomposed in plane waves through
Fourier transform (first part of equation 3.35, P (~k, t0)). Later, the extrapolation operator
is applied to each plane wave, specifically the operator that evolute the phase is in the
middle of equation 3.35, ei(V (~x,~k)|~k| t). Finally, the wavefield at a time t+∆t is synthetized
14
by the inverse Fourier transform (third part of equation 3.35)
t 0
V(x,k) ∆ t
x
k
t 0 + t∆
Figure 3.1: Illustration of the extrapolation algorithm of the wavefield in anisotropic
media. Firstly, the wavefield at a time t (inner blue curve) is decomposed in plane waves
(red tangent lines). Secondly, each plane wave is propagated during the interval ∆t with
its respective phase velocity V (~x,~k) for each point in the wavefront. Finally, the wavefield
is synthesized (reconstructed) at time t + ∆t (outer blue curve). Figure from Baltazar
and Carvalho (2018).
3.3 PHASE VELOCITY IN ANISOTROPIC MEDIA
The phase velocity is obtained from the dispersion relation (relation between frequency
and wavenumber)
ω = ω(~k) . (3.38)
If we substitute the plane wave solution in any wave equation, we will have two equations:
eikonal and transport equations. The first one tell us where the wavefront is and the second
one tell us the intensity of the event. The dispersion relation can be found by substituting
the slowness vector ~p = ~k
ω
in the eikonal equation (Musgrave (1970); Cerveny (2005)).
Once we have the dispersion relation we can find the phase velocity for an specific medium
V (~k, ω) =
ω(~k)
~k
. (3.39)
In this work, the algorithm was applied to simulate the wave propagationin Tilted Trans-
verse Isotropic (TTI). The phase velocity for q-P waves in TTI media can be precisely
calculated through expression:
V 4−(a11(~n ·~n−(~n ·~ν)2)+a33(~n ·~ν)2)V 2 +(a11a33−(a13 +2a55)2)(~n ·~n−(~n ·~ν)2)(~n ·~ν)2 = 0 ,
(3.40)
where aIJ = CIJ/ρ indicates the elastic parameters normalized by the density (Musgrave
(1970); Cerveny (2005)), ~ν is the unitary vector pointing in the direction of symmetry and
15
~n indicates the direction of the wavenumber vector, in other words ~k //~n. Considering
the pseudo-acoustic approximation (Alkhalifah, 1998), where a55 → 0 keeping the sum
a13 + 2a55 fixed, and using the Thomsen’s parameters (Thomsen, 1986), ε ≡ a11−a332a33 ,
δ ≡ (a13+a55)
2−(a33−a55)2
2a33(a33−a55) , Co ≡ a33, Cε ≡ Co(1 + 2ε), Cδ ≡ Co(1 + 2δ), equation 3.40 (also
called phase velocity surface) can be rewritten in the following way
V 2− (Cε(~n · ~n− (~n · ~ν)2) +Co(~n · ~ν)2)V +Co(Cε−Cδ)(~n · ~n− (~n · ~ν)2)(~n · ~ν)2 = 0 (3.41)
If we solve this second order polynomial equation using Bhaskara’s formula, we will
have an expression for the phase velocity at each point ~x and direction ~k
V (~x,~k) =
1
2
[
Cε (~n · ~n)− (~n · ~ν)2 + Co (~n · ~ν)2
]
±
√
1
4
[
Cε (~n · ~n)− (~n · ~ν)2 + Co (~n · ~ν)2
]2 − [Co (Cε − Cδ) (~n · ~n)− (~n · ~ν)2]
(3.42)
Equation 3.42 combined with
p(~x, t0 + t) + p(~x, t0 − t) = 2
∫
P (~k, t0)e
i~k~xcos
([
V (~x,~k) . |~k|
]
. t
)
d~k (3.43)
are the two main equations of the implemented algorithm. The next sections are devoted
to explain how the low rank approximation uses both equations.
4 LOW RANK APPROXIMATION
The low rank approximation (Fomel et al., 2013) starts by assuming we can express
the extrapolation operator
p(~x, t+ ∆t) =
∫
P (~k, t)eiφ(~x,
~k,∆t)d3~k (4.1)
in the following way
pm(~x, t+ ∆t) ≈
∑
n
Wm,n Pn( ~kn, t) e
i~kn·~xm∆k1∆k2 , (4.2)
where ∆k1 and ∆k2 are the uniform sampling interval in the first and second dimension
respectively in the wavenumber domain, pm(~x, t + ∆t) is the discretized wavefield in the
space domain, Pn( ~kn, t) is the discretized wavefield in the wavenumber domain and Wm,n
is the discretization of the evolution operator, where the blue operator will be called
extrapolation matrix defined as
W (~x,~k) = eiφ(~x,
~k,∆t) = Wm,n = e
iV (m,n)| ~kn|∆t . (4.3)
The storage of this matrix is prohibited even in 2D models. For example if we use a model
with dimensions nx = 451 and nz = 901, and 50 points on the absorption boundary we
will have the extended dimensions nxx = 551 and nzz = 1001. Since it needs to be
transformed by FFT , the Fourier dimensions would be nxft = 1024 and nzft = 1024.
Then, the matrix Wm,n would have dimensions m = nxx ∗ nzz and n = nxft ∗ nzft, so
m = 551551 and n = 1048576.
The low rank approximation tries to represent the blue operator by three new matrices
and to factorate the extrapolation operator into
W (~x,~k) ≈
M∑
m=1
N∑
n=1
W (~x, ~km)amnW ( ~xn, ~k) . (4.4)
Wm,n ≈
M∑
m=1
N∑
n=1
W 1m,MaM,NW
2
N,n . (4.5)
where M << m and N << n. M and N are the ranks chosen from the pivoted QR
decompositions after we have chosen how many columns would be selected. For example,
we may enter with a first guess rank equal to 10 and after the pivoted QR decomposition,
the algorithm outputs a rank equal to 7. Since it has an harmonic structure and redundant
information, the rank depends on the ∆t, so the smaller the ∆t, the less ranks is needed
16
17
from this matrix. This is the trade-off between accuracy and efficiency.
Assuming that this decomposition can be done, the evolution of the wavefield in the
discretized grid can be written in the form
pm(t+ ∆t) =
∑
M
W 1m,M
∑
N
aM,N
[∑
n
W 2N,n Pn(t) e
i~kn·~xm∆k1∆k2
]
. (4.6)
The sum inside the brackets in equation above represents the inverse Fourier transform
and it can be calculated by the Fast Fourier transform algorithm (Press et al., 2007).
So, to evaluate the wavefield pm(t + ∆t) is necessary the calculus of N inverse Fourier
transforms. Observing that Pn(t) is the space Fourier transform of pm(t),
Pn(t) =
1
(2π)2
∆x1∆x2
∑
q
pq(t)e
−i~kn·~xq , (4.7)
so, is necessary to calculate (N+1)FFTs. The algorithm to seismic modeling in anisotropic
media consists in applying equation 4.6 recursively during the total time of simulation.
The figure 4.1 summarize the idea of the algorithm.
Figure 4.1: The low rank algorithm takes the extrapolation operator and decompose into
three smaller matrices.
The way to construct both W (~x, ~kM) = W 1 and W ( ~xN , ~k) = W 2 is described in the
the original algorithm found in Fomel et al., 2013. Let ε be a prescribed accuracy of this
separated representation and rε be the numerical rank ofW . The algorithm for computing
W ≈ W 1 aW 2 takes the following steps:
1. Pick a uniformly random set S of W . This set needs to have β · rε columns of W ,
so an integer multiple of the columns. In practice β is chosen to be 3 or 4. Perform
the pivoted QR factorization of (W (:, S))T = (W (:, S))∗, (Golub and Van Loan,
2012). The first rε pivoted columns correspond to rε rows of the matrix W (:, S).
18
Define W 1 to be the submatrix W that consists in these rows and set x1, ..., xn ,
with N = rε to be the corresponding x values of these rows (see section 4.1).
2. Pick a uniformly random set T of W . This set needs to have α · rε rows of W , so an
integer multiple of rows. Perform the pivoted QR factorization of W (T, :). Define
W 2 to be the submatrix ofW that consists in these columns and set k1, ..., km , with
M = rε, to be the corresponding k values of these columns (see section 4.1).
3. Set the middle matrix A = W †( ~xN , ~kM), with 1 ≤ N ≤ n , 1 ≤M ≤ n (see 4.2)
4. Combine the result of the previous three steps to obtain the required separated
representation W ≈ W 1 aW 2.
4.1 LOW RANK MATRIX, W 1 AND W 2
Given a complex matrix W = [Wm,n], with m ∈ 1, ...,ml and n ∈ 1, ..., nc, its low rank
approximation is defined by
WLow = W
1aW2
H (4.8)
where H indicates the complex conjugate transpose, the columns of W 1m,M and W 2n,N are
orthornormals, aM,N has a complete rank equal to minMl, Nc, while M ∈ {1, ...,Ml} and
N ∈ {1, ..., Nc}; such that Ml, Nc << ml, nc and L2-norm ‖W −WLow‖ < ε.
The calculation of matrixW can be implemented using the QR decomposition, (Golub
and Van Loan, 2012), applied sequentially to the column space and row space of W
in other to estimate the minimum values of Ml and Nc that satisfy the criteria of the
approximation and thus determine the matrices W1 and W2.
The low rank algorithm described by Fomel et al. (2013), was used to determine the
Ml and Nc. The algorithm starts with randomly choosing the N columns ofW. In general
N ≈ 3Nc. The Nc is the expected rank to the column space of W with the same rank as
W1.
The minimum rank, Nc is estimated observing the absolute value of pivots during the
QR decomposition. If the ratio between the biggest pivot and the lowest one is less than
10−5 the QR decomposition is stopied and the Nc value is the number of pivots computed
immediately before the criteria is satisfied. After the iterations we obtain the following
approximation to the matrix W
Wc = Q
1Rc = W
1
cRc (4.9)
where W1c is equal to the first Ml columns of matrix Q1. The same similar process is
applied to WH which gives us the following approximation
WHl = Q
2Rl = W
2H
l R
H
l (4.10)
19
where W2Hl is equal to the first Nl columns of matrix Q2. Once we have the low rank
generators of the column space (inside the columns of W2) and row space (inside the rows
of W1) we need to determine the low rank matrix a using equation 4.8.
4.2 KACZMARZ ALGORITHM
In order to reduce the storage cost to compute matrix a we use the the Kaczmarz
algorithm (Kaczmarz, 1993). This algorithm has been used to solutions of robust sparse
linear systems. In each iteration the matrix a is updated from a unique value sampled
randomly from Wm,n. In each iteration the following problem of projection is applied:
argmin{ak}
1
2
(
ak − ak−1
)H (
ak − ak−1
)
(4.11)constrained by Wi,j =
∑
m
∑
n
W 1i,ma
k
m,nW
2∗
n,j (4.12)
where ∗ indicates the complex conjugate. The Lagrangian associated to the problem is
L(akp,q, λ) =
Objective Function
1
2
∑
m
∑
n
(akm,n − ak−1m,n)∗(akm,n − ak−1m,n)+
AdjointOperator
λ∗i,j
Constrain(
Wi,j −
∑
m
∑
n
W 1i,ma
k
m,nW
2∗
n,j
)
.
(4.13)
The extreme points are given by the condition
∂L
∂akp,q
=
(
akp,q − ak−1p,q
)∗ − λ∗i,jW 1i,pW 2∗q,j = 0 . (4.14)
This results in the following expression to the Lagrange multipliers
λi,j =
Wi,j −
∑
m
∑
nW
1
i,ma
k−1
m,nW
2∗
n,j(∑
mW
1
i,mW
1∗
i,m
) (∑
nW
2
n,jW
2∗
n,j
) (4.15)
resulting in the following expression to update the matrix
akp,q = a
k−1
p,q + λi,jW
1∗
i,pW
2
q,j . (4.16)
In the implementation, we randomly select Ml logml rows and Nc log nc columns of
matrix a. This prescription greatly reduce the complexity of the algorithm. We use
Ml logml x Nc log nc random samples to estimate a. Baltazar and Carvalho (2018), have
written the pseudo-code described below in algorithm 1:
20
Algorithm 1: Kaczmarz algorithm to compute low-rank matrix a.
input : Irank, Jrank, Ml, Nc, W 1, W 2, Operator(i, j)
output: a
/* Initialization : */
Niter ← 10MlNc;
error ← 0;
a← 0;
TOL← 10−6;
/* compute auxiliar arrays : */
W 1norm(1 : Ml)← 0;
W 2norm(1 : Nc)← 0;
for m← 1 to Ml do
for i← 1 to Irank do
W 1norm(m)← W 1norm(m) +W 1i,mW 1
∗
i,m
for n← 1 to Nc do
for j ← 1 to Jrank do
W 2norm(m)← W 2norm(m) +W 2i,mW 2
∗
i,m
/* Kaczmarz iterations : */
for iter ← 1 to Niter do
errorp ← error;
error ← 0;
for i← 1 to Irank do
for j ← 1 to Jrank do
Wi,j ← Operator(i, j);
res = Wi,j;
for m← 1 to Ml do
for n← 1 to Nc do
res← res−W 1i,mam,nW 2
∗
n,j;
if |res| > 0 then
error ← max(error, |res|);
λ← res
W 1norm(m)W
2
norm(n)
for m← 1 to Ml do
for n← 1 to Nc do
am,n ← am,n + λW 1∗i,mW 2n,j;
if |error − errorp| < TOL× errorp then
EXIT
5 WAVEFIELD DECOMPOSITION BASED REVERSE TIME
MIGRATION
Reverse Time Migration (Baysal et al., 1983; Whitmore, 1983; McMechan, 1983) is
now a standard imaging technique in complex environments thanks to the computational
improvement capabilities. The basic idea of this technique is the zero-lag crosscorrelation
of the forward-propagated source wavefield and the backward-propagated receiver wave-
field. A practical overview of how it works is described in Robein (2016), and it can be
found some of the advantages of RTM over ray based methods and one-way wave-equation
methods such as: stability in complex models, no dip limitations, addresses multipathing,
handles prism waves and turning waves. This last one can become a problem because
all waves are cross-correlated during the conventional imaging condition creating a low
frequency noise and false images. To overcome these problems Fletcher et al. (2006)
proposed a way to reduce the low-frequency noise. Yoon and Marfurt (2006) proposed
the Poynting vector to improve image quality. Fei et al. (2010) proposed to separate the
up- and downgoing waves in the propagated source and receiver wavefields. The main
idea is to apply the imaging condition only to the desired wave components, for example,
between the downgoing source and upgoing receiver wavefields. This method can be com-
putationally inefficient if not done properly. Liu et al. (2011) proposed a computational
efficient method to separate up- from down-going waves based on the Hilbert transform
which will be discussed below.
5.1 RTM PROBLEMS
The image I(~x) obtained from conventional imaging condition is defined by
I(~x) =
∫ tmax
0
s(~x, t)r(~x, t)dt . (5.1)
Where ~x = (x, y, z), tmax is the total time of the seismic experiment, s(~x, t) is the modeled
wavefield called source wavefield and r(~x, t) is the recorded data propagated backward in
time from receiver locations generating the receiver wavefield. Since the RTM uses two-
way wave equation, s(~x, t) and r(~x, t) contain both the downgoing sd(~x, t) and rd(~x, t)
and upgoing su(~x, t) and ru(~x, t) wavefields. Thus, the source and receiver wavefields can
be represented as
s(~x, t) = sd(~x, t) + su(~x, t) (5.2)
and
r(~x, t) = rd(~x, t) + ru(~x, t) (5.3)
21
22
When we perform the conventional imaging condition 5.1, all these wavefields are cross-
correlated as indicated below
I(~x) =
∫ tmax
0
sd(~x, t)ru(~x, t)dt+
∫ tmax
0
su(~x, t)rd(~x, t)dt
+
∫ tmax
0
sd(~x, t)rd(~x, t)dt+
∫ tmax
0
su(~x, t)ru(~x, t)dt ;
(5.4)
I(~x) = I1(~x)
Desired Image
+ I2(~x)
False Image
+ I3(~x)
Low-freq-artifact
+ I4(~x)
Low-freq-artifact
. (5.5)
Fei et al. (2015) explains the false imaging problem indicated by figure 5.1. This can be
more harmful than low frequency artifacts because they can generate false reflectors that
can be misinterpreted as geology.
(a) Representation of the false image prob-
lem. When we cross-correlate the complete
wavefield (a), we observe false images rep-
resented by the red arrows in (b)
(b) The red arrows indicate the false im-
ages. The wavefield decomposition heals
the false image splitting at yellow arrow lo-
cation.
Figure 5.1: Representation and example of the false image problem. The red arrows
indicate how these artifacts can contaminate the image if not removed. Figures modified
from Fei et al. (2015)
Kaelin and Carvajal (2011) exemplify the low-frequency artifacts that appear when
cross-correlating all the wavefields. They look like smiles on the top of the salt and can
pollute the image above it.
23
Figure 5.2: Representation of the low-frequency artifacts. The top of the salt body
presents smiles formed by low-frequencies waves cross-correlations. Figure from Kaelin
and Carvajal (2011).
In the next section we will see how to overcome these two problems by separating the
wavefields using the Hilbert transform and complex traces.
5.2 WAVEFIELD DECOMPOSITION USING COMPLEX TRACES
A complex analytic trace or extended complex trace εz [s(z, t)] is composed with the
original signal s(z, t) as the real part and the Hilbert transform H of the s(z, t) as the
imaginary part
εz [s(z, t)] = s(z, t) + iHz [s(z, t)] . (5.6)
Where s(z, t) stands for the source wavefield and r(z, t) stands for the receiver wavefield
which can also be extended using operator ε to build the receiver complex wavefield
εz [r(z, t)] = r(z, t) + iHz [r(z, t)] , (5.7)
using the relationship between Fourier transform F and Hilbert transformH, e.g., F [H[u(z or t)]] =
− i sgn(kz or ω)F [u(z or t)], where sgn(kz or ω) can be defined according to the domain
that the transformation is applied
24
sgn(kz or ω) =

1, kz or ω > 0
0, kz or ω = 0
−1, kz or ω < 0
(5.8)
When using a complex pulse to simulate the wave propagation, the source wavefield is
complex. The main advantage of using complex traces is that the spectrum has only
positive wavenumbers. Mathematically this can be expressed by
Fz {εz [s(z, t)]} =
{
2S(kz, t), kz ≥ 0
0, kz < 0
(5.9)
and the same principle works for the receiver wavefield
Fz {εz [r(z, t)]} =
{
2R(kz, t), kz ≥ 0
0, kz < 0
(5.10)
The decomposition of the wavefields in up- and down-components is very easily im-
plemented in space Fourier domain. The only requirement is to select either the positive
wavenumber (going downwards) or negative wavenumbers (going upwards). Since
the low rank approximation already executes Fourier transform from space domain (~x) to
wavenumber domain (~k) the expression for both up- and down-components of the source
wavefield are
Sdown(kz, t) =
{
S(kz, t), kz ≥ 0
0, kz < 0
(5.11)
and
Sup(kz, t) =
{
S(kz, t), kz < 0
0, kz ≥ 0
(5.12)
Applying the same idea, the receiver wavefield can be decomposed by selecting the signal
of the wavenumbers
Rdown(kz, t) =
{
R(kz, t), kz ≥ 0
0, kz < 0
(5.13)
Rup(kz, t) =
{
R(kz, t), kz < 0
0, kz ≥ 0
(5.14)
As we have seen above, the best image is obtained when cross-correlated only the down-
going source wavefield with the upgoing receiver wavefield
I(~x) =
∫ tmax
0
sdown(~x, t)rup(~x, t)dt . (5.15)
6 RESULTS
We have tested our methodology intwo models: the Marmousi model and the BP
TTI model. On both models we have compared the results of three different algorithms:
the image without wavefield decomposition, the image with wavefield decomposition and
a finite differences implementation. The difference between the first two images is an
additional Hilbert transform in time to make each trace of the back-propagated wave-
field become an analytic trace and an Hilbert transform in z direction to decompose the
forward-propagated wavefield.
6.1 MARMOUSI MODEL
Below are four different screen shots of the wavefield. In (a) and (b) are the back-
propagated wavefields injected at the receiver locations, without and with decomposition,
respectively. In (c) and (d) are the forward-propagated wavefields injected at the source
location, without and with decomposition.
(a) Backpropagated without decomposition. (b) Backpropagated with decomposition.
(c) Forward propagated without decomposition. (d) Forward propagated with decomposition.
Figure 6.1: In (a) the back-propagated wavefield without decomposition, in (b) the back-
propagated wavefield decomposed, in (c) the forward wavefield without decomposition
and in (d) the forward wavefield decomposed. Both x and z axis are in meters.
25
26
The Marmousi model represents a complex geological setting, containing horizontally
layered horizons and a series of normal faults. Besides that its tilted blocks complicates
the model towards its center.
The figures below shows all three results obtained with RTM imaging plus the velocity
model. Overall the three methodologies recovered the main reflectors efficiently, despite
some punctual differences that will be discussed in the next subsections.
0
1000
2000
D
e
p
th
 (
m
)
0 2000 4000 6000 8000
Inline (m)
Marmousi Model
2000
3000
4000
5000
6000
m
/s
(a)
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Not Decomposed
0
2000
4000
6000
(b)
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Finite Differences
-5
0
5
x107
(c)
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Decomposed
-1
0
1
x104
(d)
Figure 6.2: All three results and the velocity model in (a), (b) is the RTM without
the decomposition, (c) is the RTM using finite differences and (d) the RTM using the
decomposition.
27
6.1.1 Finite Differences x Low Rank Not Decomposed
A point to highlight is the magnitude and the range of the amplitude values, which
make the image appear brighter. Even without the decomposition of the wavefield, the
low rank algorithm could remove the low frequency effects close to the source and receivers
positions. The yellow arrow shows a discontinuity in the reflector which matches the true
model, while the finite differences approach registers it wrongly as a continuous reflector.
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Finite Differences
-5
0
5
x107
(a)
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Not Decomposed
0
2000
4000
6000
(b)
Figure 6.3: In (a) is showed the RTM using the finite differences implementation and in
(b) the RTM using low rank without decomposition.
6.1.2 Low Rank Not Decomposed x Decomposed
As we can see the magnitude of the amplitudes for item (a) are not well distributed
as in (b). The blue circle shows a cleaner image for the result using the decomposed
wavefield. The biggest problem in (b) is the top of the model around x=6 km. Below the
yellow arrow in (a) we can see a better discontinuity than the reflectors in (b).
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Not Decomposed
0
2000
4000
6000
(a)
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Decomposed
-1
0
1
x104
(b)
Figure 6.4: In (a) is showed the RTM using low rank without decomposition and in (b)
the RTM using low rank with decomposition.
28
6.1.3 Low Rank Decomposed x Finite Differences
The decomposed low rank image on item (a) shows a better resolution of the reflectors
specially if we compare with the purple box on item (b). On top of item (a) the blue
circle shows an image without strong artifacts compared with item (b). The image using
decomposition could capture the top and bottom of the tilted blocks crossed by the faults.
Also, look down in the deep parts of the model where we can see that when we apply the
laplacian we loose the low frequency content. Only in few parts on top of (b) shows a
better resolution than (a), but it is contaminated with artifacts specially between inlines
4000 m and 6000 m.
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Decomposed
-1
0
1
x104
(a)
0
1000
2000
D
e
p
th
 (
m
)
2000 4000 6000 8000
Inline (m)
Marmousi Finite Differences
-5
0
5
x107
(b)
Figure 6.5: In (a) is showed the RTM using the decomposed wavefield and in (b) the
RTM using the finite differences implementation.
29
6.2 BP TTI MODEL
The 2007 BP TTI Velocity-Analysis Benchmark dataset was created by Hemang Shah
and provided as courtesy of BP Exploration Operation Company Limited ("BP"). The
region we have selected comprehends the first 23 km of the original velocity model. Below
are four different screen shots of the wavefield. In (a) and (b) are the back-propagated
wavefields injected at the receiver locations, without and with decomposition, respectively.
In (c) and (d) are the forward-propagated wavefields injected at the source location,
without and with decomposition.
(a) Backpropagated without decomposition. (b) Backpropagated with decomposition.
(c) Forward propagated without decomposition. (d) Forward propagated with decomposition.
Figure 6.6: In (a) the back-propagated wavefield without decomposition, in (b) the back-
propagated wavefield decomposed, in (c) the forward wavefield without decomposition
and in (d) the forward wavefield decomposed
The main difference between tests using the Marmousi and the BP is that in the first
we did not considered anisotropy, we considered only the velocity model. After the pseudo-
acoustic approximation described previously we ended up with four models parameters to
describe the wave propagation in the anisotropic TTI medium which are the q-P velocity,
epsilon, delta and the angle of symmetry.
30
Below is showed the velocity model in (a), the Thomsen’s parameters epsilon in (b)
and delta in (d) and the tilt angle in (c). We have selected the first 23 km of the original
model because the focus was to compare the results for the top of the salt body.
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
0 5 10 15 20
Inline (km)
BP 2007 VELOCITY MODEL 0km-23km
2
3
4
k
m
/s
(a)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
0 5 10 15 20
Inline (km)
BP 2007 EPSILON MODEL 0km-23km
0
0.05
0.10
0.15
0.20
(b)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
0 5 10 15 20
Inline (km)
BP 2007 THETA MODEL 0km-23km
-20
-10
0
10
(c)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
0 5 10 15 20
Inline (km)
BP 2007 DELTA MODEL 0km-23km
0
0.02
0.04
0.06
0.08
0.10
(d)
Figure 6.7: In (a) the velocity model, in (b) the epsilon model, in (c) the arbitrary angles
of symmetry model and in (d) the delta model. As we can see from item (c) we did not
select an area with strong variations on angle of symmetry. We focused our analysis on
the top of the salt to see if could remove the low frequency artifacts.
31
Below we see all three implementations compared with the velocity model. Since we
focused our analysis on the top of the salt, we ran shots from 205 to 268 on (b) and
(c), while in (d) we ran shots from 205 to 400. It was enough to show the improvement
not only in illumination, but also the decomposed RTM significantly removed the low
frequency noise.
Besides that, we can see that even though we ran less shots in (d), it illuminated
better the region below salt indicated by the yellow arrow. In the following subsections
we comparethemselves and show the main differences.
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
0 5 10 15 20
Inline (km)
BP 2007 VELOCITY MODEL 0km-23km
2
3
4
k
m
/s
(a)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
5 10 15 20
Inline (km)
BP Not Decomposed
0
100
200
(b)
0
0.2
0.4
0.6
0.8
1.0
x104
D
e
p
th
 (
k
m
)
0.5 1.0 1.5 2.0
x104Inline (km)
BP Finite Differences
-5
0
5
x107
Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org) (c)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
5 10 15 20
Inline (km)
BP Decomposed
-100
0
100
(d)
Figure 6.8: In (a) the velocity model, in (b) the RTM not decomposed, in (c) the finite
differences implementation and in (d) the RTM decomposed.
32
6.2.1 Finite Differences x Low Rank Not Decomposed
The main difference between these two implementations is the lack of distribution for
amplitudes for (b) which makes the image to look brighter. One thing in common in these
two implementations is the low frequency noise indicated by the green and red boxes. The
top of the salt is less contaminated in item (a).
0
0.2
0.4
0.6
0.8
1.0
x104
D
e
p
th
 (
k
m
)
0.5 1.0 1.5 2.0
x104Inline (km)
BP Finite Differences
-5
0
5
x107
Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org) (a)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
5 10 15 20
Inline (km)
BP Not Decomposed
0
100
200
(b)
Figure 6.9: In (a) the RTM using finite differences and in (b) the RTM without decom-
position
6.2.2 Low Rank Not Decomposed x Decomposed
Besides the difference in the range of amplitudes, we can compare the top of the salt
on both images. We can see that the low frequency noise on top of salt in item (a) con-
taminated the image with high amplitude values while most of the data is comprehended
close to zero values indicated by the yellow arrow in (a).
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
5 10 15 20
Inline (km)
BP Not Decomposed
0
100
200
(a)
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
5 10 15 20
Inline (km)
BP Decomposed
-100
0
100
(b)
Figure 6.10: In (a) the RTM without decomposition and in (b) the RTM using decompo-
sition.
33
6.2.3 Low Rank Decomposed x Finite Differences
On this example we can see the power of RTM using decomposition of the wavefields.
It significantly removed the low frequency noise if we compare the green box in both
images. The structures above the salt in item (b) is contaminated with migration smiles.
Even with fewer shots we can see an improvement in illumination specially indicated by
the yellow arrow in (a). Also, the right corner of the salt is shorter in item (b) compared
to (a).
0
2
4
6
8
10
D
e
p
th
 (
k
m
)
5 10 15 20
Inline (km)
BP Decomposed
-100
0
100
(a)
0
0.2
0.4
0.6
0.8
1.0
x104
D
e
p
th
 (
k
m
)
0.5 1.0 1.5 2.0
x104Inline (km)
BP Finite Differences
-5
0
5
x107
Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org) (b)
Figure 6.11: In (a) RTM using decomposition and in (b) RTM using finite differences.
7 CONCLUSIONS
We have presented a methodology combining the low rank approximation to simulate
the wave propagation in TTI media with the power of RTM to obtain a better image with
less low frequency noise in complex areas.
We have compared three main results: low rank without decomposition, with decom-
position and a finite differences implementation. The results showed that RTM using
wavefield decomposition is capable of improving image quality and imaging anisotropic
settings. The main reason was the use of the Hilbert transform to create an complex
wavefield, so we could decompose it before applying the cross-correlation.
Also, since it needs only two extra Hilbert transforms (one outside the loop of the time)
it can be implemented easily in the low rank algorithm. The low rank formulation can be
extended to other classes of anisotropy and a next step would be the implementation 3D
multi-GPU or the orthorhombic case.
34
BIBLIOGRAPHY
Alkhalifah, T., 1998, Acoustic approximations for processing in transversely isotropic
media: Geophysics, 63, 623–631.
Baltazar, L. P. R., and J. C. Carvalho, 2018, Modelagem sísmica em meios anisotrópicos
usando aproximações de baixo posto.
Baysal, E., D. D. Kosloff, and J. W. Sherwood, 1983, Reverse time migration: Geophysics,
48, 1514–1524.
Billette, F., and S. Brandsberg-Dahl, 2005, The 2004 bp velocity benchmark, in 67th
EAGE Conference & Exhibition: British Petroleum, 10–11.
Cerveny, V., 2005, Seismic ray theory: Cambridge university press.
Claerbout, J. F., 1985, Imaging the earth’s interior: Blackwell scientific publications
Oxford, 1.
Etgen, J. T., and S. Brandsberg-Dahl, 2009, The pseudo-analytical method: Application
of pseudo-laplacians to acoustic and acoustic anisotropic wave propagation, in SEG
Technical Program Expanded Abstracts 2009: Society of Exploration Geophysicists,
2552–2556.
Fei, T. W., Y. Luo, S. Aramco, and G. T. Schuster, 2010, De-blending reverse-time
migration, in SEG Technical Program Expanded Abstracts 2010: Society of Exploration
Geophysicists, 3130–3134.
Fei, T. W., Y. Luo, J. Yang, H. Liu, and F. Qin, 2015, Removing false images in re-
verse time migration: The concept of de-primaryde-primary reverse time migration:
Geophysics, 80, S237–S244.
Fletcher, R. P., P. J. Fowler, P. Kitchenside, and U. Albertin, 2006, Suppressing unwanted
internal reflections in prestack reverse-time migration: Geophysics, 71, E79–E82.
Fomel, S., L. Ying, and X. Song, 2013, Seismic wave extrapolation using lowrank symbol
approximation: Geophysical Prospecting, 61, 526–536.
Golub, G. H., and C. F. Van Loan, 2012, Matrix computations: JHU Press, 3.
Kaczmarz, S., 1993, Approximate solution of systems of linear equations: International
Journal of Control, 57, 1269–1271.
Kaelin, B., and C. Carvajal, 2011, Eliminating imaging artifacts in rtm using pre-stack
gathers, in SEG Technical Program Expanded Abstracts 2011: Society of Exploration
Geophysicists, 3125–3129.
Liu, F., G. Zhang, S. A. Morton, and J. P. Leveille, 2007, Reverse-time migration using
one-way wavefield imaging condition, in SEG Technical Program Expanded Abstracts
2007: Society of Exploration Geophysicists, 2170–2174.
——–, 2011, An effective imaging condition for reverse-time migration using wavefield
decomposition: Geophysics, 76, S29–S39.
35
36
Martinez, R., 2018, How recent advances in seismic depth imaging can enhance prospect
identification and appraisal, in SEG Honorary Lecture in Latin America: Society of
Exploration Geophysicists, 0–10.
McMechan, G. A., 1983, Migration by extrapolation of time-dependent boundary values:
Geophysical Prospecting, 31, 413–420.
Musgrave, M. J. P., 1970, Crystal acoustics: Introduction to the study of elastic waves
and vibrations in crystal: Holden-Day.
Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, 2007, Numerical
recipes 3rd edition: The art of scientific computing: Cambridge university press.
Robein, E., 2016, Reverse time migration:how does it work,when to use it, in Youtube:
EAGE: E-lectures, 15/11/2016.
Schleicher, J., and J. C. Costa, 2016, A separable strong-anisotropy approximation for
pure qp-wave propagation in transversely isotropic mediapure qp-wave propagation:
Geophysics, 81, C337–C354.
Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954–1966.
Tsvankin, I., 2012, Seismic signatures and analysis of reflection data in anisotropic media:
Society of Exploration Geophysicists.
Versteeg, R., 2005, Marmousi, model and data, in EAGE Conference & Exhibition:
EAGE, 10.3997/2214–4609.201411190.
Wards, B. D., G. F. Margrave, and M. P. Lamoureux, 2008, Phase-shift time-stepping for
reverse-time migration, in SEG Technical Program Expanded Abstracts 2008: Society
of Exploration Geophysicists, 2262–2266.
Whitmore, N. D., 1983, Iterative depth migration by backward time propagation, in SEG
Technical Program Expanded Abstracts 1983:Society of Exploration Geophysicists,
382–385.
Yan, J., and H. Liu, 2016, Modeling of pure acoustic wave in tilted transversely isotropic
media using optimized pseudo-differential operators: Geophysics, 81, T91–T106.
Yoon, K., and K. J. Marfurt, 2006, Reverse-time migration using the poynting vector:
Exploration Geophysics, 37, 102–107.
Zheng, Y., Y. Wang, and X. Chang, 2018, 3d forward modeling of upgoing and downgoing
wavefields using hilbert transform: Geophysics, 83, F1–F8.