EAGE2016.tex

% 78th EAGE Conference & Exhibition 2016
% Vienna, Austria
% 30 May - 2 June 2016
%-----------------------------------------------------------

% Authors:
% ----------------------------------------------------------
% Wanderson C. Ferreira
% Henrique Bueno dos Santos

% December 2015 - January 2016

% ----------------------------------------------------------
\documentclass{vie16}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage{wrapfig}
\usepackage{float}
\usepackage{amsmath}
\usepackage[table,xcdraw]{xcolor}
\usepackage{booktabs}
\usepackage{multirow}
\usepackage{adjustbox}
\usepackage{subfigure}

% Path to figures
\graphicspath{{Fig/}}

% COMENTAR ANTES DE SUBMETER!
% Comando temporarios
%\usepackage{color}
%\usepackage[normalem]{ulem}
%\newcommand{\new}[1]{\textcolor{blue}{#1}}
%\newcommand{\old}[1]{\textcolor{MyDarkViolet}{\sout{#1}}}
%\newcommand{\att}[1]{\textcolor{red}{#1}}
%\newcommand{\js}[1]{\textcolor{darkgreen}{#1}}
%\newcommand{\obs}[1]{\textcolor{orange}{#1}}
%\definecolor{MyDarkViolet}{rgb}{0.602,0,0.602}
%\definecolor{orange}{RGB}{255,127,0}
%\definecolor{darkgreen}{RGB}{0,127,0}

% For pdf hyperlinks
%\usepackage{hyperref} \hypersetup{colorlinks=true}
%\hypersetup{citecolor=blue} \hypersetup{pdftitle={Global
%optimization for AVO inversion: a genetic algorithm using a table-based
%ray-theory algorithm}}
%\hypersetup{pdfauthor={W.~C.~Ferreira, F.~Hilterman,
%L.~A.~Diogo, H.~B.~Santos, J.~Schleicher and A.~Novais}}
% For making a print without color links
%\hypersetup{linkcolor=black,citecolor=black,filecolor=black,urlcolor=black}
%\hypersetup{linkcolor=blue,citecolor=blue,filecolor=blue,urlcolor=blue}
%\usepackage{url}


% Criar fluxogramas
\usepackage[latin1]{inputenc}
\usepackage{tikz}
\usetikzlibrary{shapes,arrows} 

% definindo blocos para o fluxograma
\tikzstyle{decision} = [diamond, draw, fill=blue!20, text width=6.5em, text 
badly centered, node distance=3cm, inner  sep=0pt]
\tikzstyle{decisionmi} = [diamond, draw, fill=blue!20, 
text width=1.5em, text badly centered, node distance=3cm, inner sep=0pt]
\tikzstyle{block} = [rectangle, draw, fill=blue!20, 
text width=16em, text centered, rounded corners, minimum height=4em]
\tikzstyle{blockm} = [rectangle, draw, fill=blue!20, 
text width=6em, text centered, rounded corners, minimum height=4em]
\tikzstyle{blockmg} = [rectangle, draw, fill=blue!20, 
text width=10em, text centered, rounded corners, minimum height=2em]
\tikzstyle{blockmi} = [rectangle, draw, fill=blue!20, 
text width=3em, text centered, rounded corners, minimum height=2em]
\tikzstyle{blockf} = [rectangle, draw, fill=blue!20, 
text width=16em, text centered, rounded corners, minimum height=4em]
\tikzstyle{line} = [draw, -latex']
\tikzstyle{cloud} = [draw, ellipse,fill=blue!20, node distance=3cm,
minimum height=3em]
\tikzstyle{cloudr} = [draw, ellipse,fill=blue!20, node distance=3cm,
minimum height=3em]

\begin{document}
\title{Global optimization for AVO inversion: a genetic algorithm using
a table-based ray-theory algorithm}
\author{W.~C.~Ferreira, F.~Hilterman, L.~A.~Diogo, H.~B.~Santos,
J.~Schleicher and A.~Novais}
\maketitle

\begin{abstract}
Amplitude Variation with Offset (AVO) inversion provides
estimates of the P-wave velocity, S-wave velocity and density of
a stratified medium. Global optimization is desirable
for the inversion to account for the multi-parametric behaviour
of the AVO inversion which is strongly affected by the initial
estimates of the model rock properties.  We carried out an analysis to verify
the dependency between P-wave, S-wave velocity and density in
the recovered parameters using empirical
relations as constraints. In inversion schemes, the forward
modelling is often the most time consuming process. To reduce
computation time, we have implemented a genetic algorithm using a
\textit{table-based} ray-theory algorithm to allow for 
a large amount of models in the global search. Our
results show that the genetic algorithm was capable of recovering
the physical parameters with good agreement for examples
using the empirical constraints. However, it sometimes
converged to solutions which were far from the correct answer,
but were good models to explain the observed dataset. The
forward modelling algorithm has shown excellent performance to
be used in global optimization schemes, because it allows the use of a
large number of members in the population of the genetic
algorithm.
\end{abstract}

\newpage
\section{Introduction}
The methodology of AVO has been widely used in the industry
as a direct hydrocarbon indicator. A more quantitative
interpretation can be achieved by means of AVO inversion for
the rock properties. For the forward modeling part, one
classically calculates the reflection coefficients for plane
waves as a function of incident angle (offset) using the
formulas of \cite{Knott1899} and \cite{Zoeppritz1919}. For
the reverse model, \cite{Rosa1976} derived and verified the
ill-posed nature of the inversion of Zoeppritz equations for
rock properties, meaning that we have multiple combinations
of input to produce the same output. Besides, due to its
multi-parametric formulation, the solution space is very
complex with many local minima, which makes it difficult to
find correct solutions. For those reasons,
\cite{Stoffa1991}, \cite{Mallick1995} and others suggest
global optimization schemes to treat the AVO
inversion. However, the process of global optimizations
requires many forward models which considerably increase the
computation time. Many of the works published in literature
use forward modelling algorithms based on the wave equation,
which demand great computational power. This study is based
on the premise that unconsolidated sediments with small
rock-property changes between layers can be modeled with ray
theory without a loss of resolution. For computational ease
we have developed a forward modeling algorithm called
\textit{table-based ray theory} in order to speed up the
computation of synthetic seismograms for global
optimization.  The genetic algorithm is based on findings of
\cite{Stoffa1991} and \cite{Sen1992} and uses the
\cite{Gardner1974} and \cite{Castagna1985} empirical
relationships as constraints to verify the dependency among
the parameters.

\section{Synthetic common-midpoint seismogram: forward
modelling} 

The forward modelling algorithm implemented is the \textit{table-based
ray tracing}. It assumes an Earth model composed of $n$ horizontal
isotropic layers. The parametric equations to compute the two-way
traveltimes and offset as a function of a constant ray parameter are
given by \citep{Slotnick1959}
%!
\begin{eqnarray} x & = & 2\sum_{i=1}^{n} \frac {h_{i}pv_{i}} {(1
		- p^{2}v_{i}^{2})^{1/2}} \ ,\label{eq.x} \\ t & = &
		2\sum_{i=1}^{n} \frac {h_{i}} {v_{i}(1 -
		p^{2}v_{i}^{2})^{1/2}} \ , \label{eq.t} \end{eqnarray}
%
where $h_{i}$, $v_{i}$ and $p$ are, respectively, the
layer thickness and velocity of the $i$th layer, and the ray parameter 
defined by
%!
\begin{equation}
p = \frac{\sin(\theta_{i})}{v_{i}}\ .\label{eq.p}
\end{equation}
%!
Here $\theta_i$ is the angle between the seismic ray and the
vertical in the $i$th layer. Using a sonic log or other
vertical velocity information, we build three different
two-dimensional tables (\textit{offset, traveltime and
reflection coefficient}) starting with equations
\ref{eq.x} and \ref{eq.t}.

To start the process, the P-wave velocity as a function of depth
is transformed to time. The same applies to the S-wave velocity
and density. For the 2D traveltime table, there are 90 columns
associated with the ray parameter when the incident angle in the
first medium varies from $0^{\circ}$ to $89^{\circ}$. The rows are
associated with the $t_{o}$ times in increments of the time
sample interval. Now, the 2D tables for offset and traveltime
can easily be determined using equations \ref{eq.x} and
\ref{eq.t} respectively. Since the velocity function is sampled
in equal increments of the sample rate the same way the 2D
tables are, the computation of the incident angle for each cell
in the 2D table is possible with equation \ref{eq.p}. With the
incidence angle known, the next step is to compute the 2D table
of reflection coefficients at each interface. This was
accomplished with Shuey's approximation \citep{Shuey1985} given by
%!
\begin{equation}
R_{P} (\phi) = A + B\sin^{2}(\theta) + C\sin^{2}(\theta)\tan^{2}(\theta)
\ , \label{eq.shuey}
\end{equation}
%!
where $\phi$ is the incident angle, $\theta$ is the average of
the incident and transmitted angles, and the terms $A$, $B$ and
$C$ are defined as
%?!
\begin{equation}
\begin{split}
A & = \frac{1}{2} \left(\frac{\triangle V_{P}}{V_{P}}     +
\frac{\triangle \rho}{\rho}     \right)\ , \\
B & = \frac{1}{2} \frac{\triangle V_{P}}{V_{P}}  - 2\left(\frac{
V_{S}}{V_{P}}\right)^{2}  \left(\frac{2\triangle V_{S}}{V_{S}} +
\frac{\triangle \rho}{\rho}\right)   \ , \\
C & = \frac{1}{2} \left(\frac{\triangle V_{P}}{V_{P}} \right) \ .
\end{split}
\end{equation}
%!
In these expressions, $\triangle V_{P} = V_{P2} - V_{P1}$,
$\triangle V_{S} = V_{S2} - V_{S1}$, $\triangle \rho = \rho_{2}
- \rho_{1}$ are the differences between the P and S-wave velocity and
density values, respectively, across an interface and $V_{P}$, $V_{S}$
and $\rho$ are their arithmetic averages.

In order to build a synthetic common-midpoint (CMP) section,
we define which are the offsets to model and start a
search row-by-row inside the \textit{offset table} looking for
the requested offsets.  Once the cell is found for the specified
offset, we go into the same cell position inside the
\textit{traveltime} and \textit{reflectivity tables} to place
the reflection coefficient from that position at the correct
time arrival inside the trace.

Because the algorithm is very simple, it provides the potential of being
adapted in the future to study more complex effects such as the stretch
in normal moveout (NMO) correction, NMO without stretch
(ray-trace NMO correction), array effects which affect the
amplitude of the wave received by the streamer, 
extraction of AVO attributes, polar anisotropic media, etc.
 
\section{Global optimization approach: Genetic Algorithm}
\begin{wrapfigure}{l}{0.49\textwidth}
 \centering
 \begin{tikzpicture}[thick,scale=0.53, every 
 node/.style={scale=0.53},node distance = 2.4cm, auto]
    % Place nodes
    \node [block] (init) {Build initial population};
    \node [cloud, above of=init, node distance=2cm] (expert) {Genetic 
    Algorithm};
    \node [block, below of=init] (identify) {Apply \textbf{Forward modelling}};
    \node [blockf, below of=identify] (pop) {Apply \textbf{objective 
    function}}; 
    \node [decision, below of=pop] (evaluate) {Converged?};
    \node [blockm, below of=evaluate, node distance=3cm] (select) 
    {Selection};
    \node [blockm, below of=select] (reprod) {Reproduction};
    \node [cloudr, right of=reprod, node distance=6cm] (crossover) 
    {Crossover};
    \node [cloudr, above of=crossover, node distance=2cm] (mutation) 
    {Mutation};
    \node [cloudr, above of=mutation, node distance=2cm] (update) {Update};
    \node [blockm, above of=update, node distance=2cm] (newpop) {New 
    population};
    \node [blockmi, below of=reprod, node distance=2cm] (stop) {Stop};
    \node [blockmi, left of=select, node distance=4cm] (sim) {Yes};

    % Draw edges
    \path [line] (init) -- (identify);
    \path [line] (identify) -- (pop);
    \path [line] (pop) -- (evaluate);
    \path [line] (expert) -- (init);
    \path [line] (evaluate) --  node {no} (select);
    \path [line] (select) -- (reprod);
    \path [line] (reprod) -- (crossover);
    \path [line] (crossover) -- (mutation);
    \path [line] (mutation) -- (update);
    \path [line] (update) -- (newpop);  
    \path [line] (newpop) |- (identify);
    \path [line] (evaluate) -| (sim);
    \path [line] (sim) |- (stop);
\end{tikzpicture}
\caption{Flowchart of the genetic algorithm.}
% implemented in FORTRAN 77/90.}
\label{fig.fluxo}
\end{wrapfigure}

The genetic algorithm (Figure~\ref{fig.fluxo}) is a technique
used to perform global optimizations based on the natural
selection process \citep{Holland1975}. In order to simulate this
process, the algorithm constantly modifies the initial
pseudo-random population in order to reach local minima
positions.  Initially, each member of the pseudo-random
population is a potential solution to the problem and after some
iterations the genetic algorithm guides the population to
the best fit positions.

The algorithm needs to \textit{select} a percentage of the
population to start the reproduction scheme where the
\textit{crossover} and the \textit{mutation} are the two classic
techniques to exchange genetic information between the members
and also randomly change the genetic information to provide
exploration of the solution space. In order to select the
members of the population we need to measure the fitness between
each potential solution and the data we want to optimize. For this
purpose, we used the objective function of \cite{Porsani2000}
given by
%!
\begin{equation}
h = \frac{2y^{T} x} {y^{T}y + x^{T}x} \ ,
\label{eq.Porsani}
\end{equation}
%!
where, $x$ and $y$ are the observed and
modelled data in the time domain and  $x^{T}$, $y^{T}$ are their
transposes. In order to constrain the inverse problem due to the
dependency between S-wave velocity and density with P-wave
velocity, we used the Gardner \citep{Gardner1974} and
Castagna \citep{Castagna1985} relations to write S-wave velocity and
density as a function of P-wave velocity.

\section{Results}
We used a synthetic model in order to verify the global
optimization method for AVO inversion and its relationships with
the prior information, i.e., the \cite{Gardner1974} and
\cite{Castagna1985} constraints, and also to analyze the
performance of the investigated forward modelling algorithm. The
physical parameters for a simple synthetic model are described
in the Table~\ref{tab.model}. The thickness between the two
inner layers was chosen to be small enough to allow interaction
between the wavelets.

% Please add the following required packages to your document preamble:
% \usepackage{booktabs}
% \usepackage[table,xcdraw]{xcolor}
% If you use beamer only pass "xcolor=table" option, i.e. 
%\documentclass[xcolor=table]{beamer}
%\begin{table}[]
%\centering
%\caption{My caption}
%\label{my-label}
\begin{wraptable}{l}{0.3\textwidth}
\caption{Model parameters to test the inversion scheme.}
\begin{adjustbox}{max width=0.3\textwidth}
\label{tab.model}
\begin{tabular}{@{}
>{\columncolor[HTML]{FFFFFF}}c 
>{\columncolor[HTML]{FFFFFF}}c 
>{\columncolor[HTML]{FFFFFF}}c 
>{\columncolor[HTML]{FFFFFF}}c 
>{\columncolor[HTML]{FFFFFF}}c @{}}
\toprule
\textbf{Layers} & \textbf{\begin{tabular}[c]{@{}c@{}}Thickness\\ 
(m)\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Vp \\ 
(m/s)\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Vs\\ 
(m/s)\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Density\\ 
(g/$\mathbf{c^{3}}$)\end{tabular}} \\ \midrule
1 & 1000 & 2000 & 551.72 & 2.07 \\
2 & 50 & 2800 & 1241.38 & 2.25 \\
3 & 50 & 2300 & 810.34 & 2.14 \\
4 & 500 & 3000 & 1413.80 & 2.29 \\ \bottomrule
\end{tabular}
%\end{table}
\end{adjustbox}
\end{wraptable}

\vspace{-0.4cm}
A synthetic data set was modelled with the table-based ray theory
algorithm using the parameters from Table \ref{tab.model},
%using XXX modeling %AQUI FALTA DIZER COMO FOI MODELADO ???  Foi
%modelado com o algoritmo de tabelas descrito no começo
a zero-phase Ricker wavelet with central frequency of 40~Hz, and a time
sample rate of 2~ms. In order to perform the global optimization, each
CMP gather modelled during each genetic algorithm run is NMO corrected.
%QUAL E ESSE MODELO ALEATORIO? ELE NAO FOI MENCIONADO ANTES. ???  !
% Aquela parte estava sobrando no texto. Devo ter escrito e apagado
% algo e ela sobrou. O que eu queria dizer nessa frase era que cada
% modelo gerado dentro do algoritmo genetico eh corrigido de NMO.
The parameters used in the \textit{genetic algorithm} are given
by a \textit{selection rate} of $50\%$, $P_{mutation}$~=~0.1,
$P_{update}$~=~0.47, $P_{crossover}$~=~0.90,~1000 members in the
population and a convergence criterion of 200 iterations. To study the quality of the inversion, we concentrate on the
parameters of the fourth
layer of the model in Table~\ref{tab.model}.
Figure~\ref{fig2} shows the evolution of the genetic
algorithm solutions for the inverted P-wave velocity parameter in the
fourth layer when running with and without the
\cite{Gardner1974} and \cite{Castagna1985} constraints. As we
can see in Figure~\ref{fig2-with}, after a few iterations the
algorithm is already guiding all the members of the population
close to the correct position in the solution space.
%!
However, Figure~\ref{fig2-without} shows that without the
additional constraints to restrict the solutions, the
algorithm is leading a significant part of the population to particular solutions
that are not very close to the correct one but are good answers
to the objective function of the inverse problem as we can see in
Table~\ref{tab.recovered}.

\begin{figure}[H]
%\subfigcapskip-5pt
\subfigtopskip-5pt
\subfigbottomskip0pt
\subfigure[]{
    \includegraphics[width=0.49\textwidth]{evol-sol-constrains-edit}
    \label{fig2-with} 
}
\subfigure[]{
    \includegraphics[width=0.49\textwidth]{evol-sol-no-constrains-edit}
    \label{fig2-without} 
}
\caption{Evolution of the genetic algorithm solutions for the P-wave 
velocity parameter in the fourth layer when running with \subref{fig2-with} 
and without \subref{fig2-without} the \cite{Gardner1974} and 
\cite{Castagna1985} constraints.
}
%% Wanderson, assim que tivermos as figuras finais a gente tem que voltar 
%%aqui e melhorar a descricao
\label{fig2}
\end{figure}
\vspace{-0.8cm}
In Table~\ref{tab.recovered} we have compiled the recovered parameters for the
%best models
models with best value of the objective function
%QUAL E A AFIRMACAO CORRETA AQUI? Acredito que seja a segunda. ???
% A segunda afirmaçao eh mesmo a correta.
%recovered 
with and without the constraints on the
parameters. The algorithm generated 301.000 forward synthetics in
25 minutes. Even with more than $20\%$ of relative error in the
average, the best model recovered without the constraints still
produce a very similar CMP gather as we can verify by the
correlation coefficient value.

% Please add the following required packages to your document preamble:
% \usepackage{booktabs}
% \usepackage{multirow}

\section{Conclusions}
The forward modelling algorithm presented has shown excellent
potential to be used in global optimization schemes by providing
good performance with restricted computer power. The parameter constraints
were helpful to reach the correct global minimum. However, the
result without constraints indicate that the procedure can still converge
to incorrect results that satisfy the objective function. In our
future research, the genetic algorithm will be formulated inside
the Bayesian framework to provide solid statistical information
about the distribution of the solutions during the iterations.
It is expected to be able to visualize distinguishable
concentrations of particular solutions which could be useful to
determine possible scenarios for the recovered parameters.
%Finally, for future work in this project will also introduce
%errors in the models and an application to real data to verify
%the reliability of this technique.



\begin{table}[h!]
	\centering
%\caption{Recovered parameters of the model in the fourth layer when 
%running with and without the Gardner and Castagna constrains.}
\caption{Recovered parameters of the model with and without constraints.}
\begin{adjustbox}{scale=0.7}
\label{tab.recovered}
\begin{tabular}{@{}cccccc@{}}
\toprule
\multirow{2}{*}{\textbf{Parameters}} & \multirow{2}{*}{\textbf{Synthetic 
Model}} & \multicolumn{2}{c}{\textbf{\begin{tabular}[c]{@{}c@{}}Recovered 
Model\\ (with constraints)\end{tabular}}} & 
\multicolumn{2}{c}{\textbf{\begin{tabular}[c]{@{}c@{}}Recovered Model\\ 
(without constraints)\end{tabular}}} \\ \cmidrule(l){3-6} 
 &  & \textbf{Parameters} & \textbf{Relative Error (\%)} & \textbf{Parameters} 
 & \textbf{Relative Error (\%)} \\ \cmidrule(r){1-2}
\multirow{4}{*}{\begin{tabular}[c]{@{}c@{}}P-wave\\ (m/s)\end{tabular}} & 
2000 & 2000 & 0.00 & 2002 & -0.10 \\
 & 2800 & 2788 & 0.43 & 2699 & 3.61 \\
 & 2300 & 2280 & 0.87 & 2394 & -4.09 \\
 & 3000 & 2971 & 0.97 & 2658 & 11.40 \\
\cmidrule{2-6}
\multirow{4}{*}{\begin{tabular}[c]{@{}c@{}}S-wave\\ (m/s)\end{tabular}} & 
551.72 & 552.71 & -0.07 & 778 & -41.01 \\
 & 1241.38 & 1231 & 0.84 & 1327 & -6.90 \\
 & 810.34 & 793 & 2.14 & 1106 & -36.49 \\
 & 1413 & 1389 & 1.75 & 1358 & 3.95 \\
\cmidrule{2-6}
\multirow{4}{*}{\begin{tabular}[c]{@{}c@{}}Density\\
	(g/$c^{3}$)\end{tabular}} 
& 2.07 & 2.069 & 0.01 & 1.535 & 25.83 \\
 & 2.25 & 2.247 & 0.13 & 1.739 & 22.71 \\
 & 2.14 & 2.134 & 0.24 & 1.541 & 28.01 \\
 & 2.29 & 2.284 & 0.25 & 1.934 & 15.31 \\
\cmidrule{1-6}
\multicolumn{1}{l}{\textbf{Correlation Coef.}} & 1 & 
\multicolumn{2}{c}{0.999887645} & \multicolumn{2}{c}{0.999847412} \\ 
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}
\section{Acknowledgements}
This research was supported by the Brazilian research agencies CAPES, CNPq, 
FAPESP and FINEP. The first author W. C. Ferreira thanks to CAPES for the 
Science Without Borders fellowship, Geokinetics and the Allied Geophysical 
Laboratories. H. B. Santos is grateful to CGG-Brazil, Petrobras, ANP and 
PRH-PB230 for his fellowships. Additional support for the authors was 
provided by the sponsors of the \textit{Wave Inversion Technology (WIT) 
Consortium}.

\bibliography{EAGE2016}
\end{document}