\documentclass{article}
\usepackage{fullpage}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage{xcolor}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage[natbibapa]{apacite}
\usepackage{eso-pic}
\AddToShipoutPictureBG{\AtPageLowerLeft{\includegraphics[scale=0.7]{powered-by-Authorea-watermark.png}}}

\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}

\usepackage{authblk}

\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}

\usepackage{listings}
\usepackage{natbib}
\usepackage{hyperref}

\lstset{ 
    showstringspaces=true, 
    backgroundcolor=\color{white}, % choose the background color 
    % basicstyle=\footnotesize, % size of fonts used for the code 
    breaklines=true, % automatic line breaking only at whitespace 
    keywordstyle=\color{black}, % keyword style 
    stringstyle=\color{black}, % string literal style 
    language=python, % Set your language (you can change the language for each code-block optionally) 
    frame=lrtb, % 
    xleftmargin=\fboxsep, % 
    xrightmargin=-\fboxsep, % 
}


\begin{document}

\title{Step-by-step NMO correction}


\author[ ]{Leonardo Uieda}

\affil[ ]{}
\vspace{-1em}


\date{}

\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup

\selectlanguage{english}
\begin{abstract}
Corresponding author: leouieda@gmail.com

This is a part of The Leading Edge ``Geophysical Tutorials'' series.
You can read more about it in \citet{Hall_2016a}.

Open any textbook about seismic data processing and you will inevitably find a
section about the normal moveout (NMO) correction.
There you'll see that we can correct the measured travel-time of a reflected
wave $t$ at a given offset $x$ to obtain the travel-time at normal
incidence $t_0$ by applying the following equation

\begin{equation}
t_0^2=t^2-\dfrac{x^2}{v_\mathrm{NMO}^2}
\label{eq:traveltime}
\end{equation}

in which $v_\mathrm{NMO}$ is the NMO velocity.
There are variants of this equation with different degrees of accuracy,
but we'll use this one for simplicity.

When applied to a common midpoint (CMP) section, the equation above is
supposed to turn the hyperbola associated with a reflection into a straight
horizontal line.
What most textbooks won't tell you is \textit{how, exactly, do you apply this
equation to the data}?

Read on and I'll explain step-by-step how the algorithm for NMO correction from
\citet{Yilmaz_2001} works and how to implement it in Python.
The accompanying Jupyter notebook \citep{Perez_2007} contains the full source
code, with documentation and tests for each function.
You can download the notebook at
\href{https://github.com/seg}{github.com/seg} or
\href{https://github.com/pinga-lab/nmo-tutorial}{github.com/pinga-lab/nmo-tutorial}.
%
\end{abstract}%



\section*{What the equation doesn't tell us}

Equation \ref{eq:traveltime} relates travel-times: the one we can measure ($t$)
and the one we want to know ($t_0$).
But the data in our CMP gather are actually a matrix of \textit{amplitudes}
measured as a function of time ($t$) and offset.
Our NMO corrected gather will also be a matrix of amplitudes as a function of
time ($t_0$) and offset.
So what we really have to do transform one matrix of amplitudes into the other.
\textit{But the equation has no amplitudes}!

This is a major divide between the formula we've all seen before and what
actually goes on in the software that implements it.
You have probably never thought about it -- I certainly hadn't --
so let's bridge this divide.
Next, I'll explain an algorithm that maps the amplitudes in the CMP to
amplitudes in an NMO corrected gather.


\section*{Doing it backwards}

It's surprisingly difficult to find a description of a method for calculating
the amplitudes in the NMO correction.
The only one I could find is a single paragraph in the book by
\citet{Yilmaz_2001} (available on the SEG Wiki at
\href{http://wiki.seg.org/wiki/NMO_for_a_flat_reflector}{wiki.seg.org/wiki/NMO\_for\_a\_flat\_reflector}):

\begin{quotation}
``The idea is to find the amplitude value at A' on the NMO-corrected gather
from the amplitude value at A on the original CMP gather. Given quantities
$t_0$, $x$, and $v_\mathrm{NMO}$, compute $t$ from equation (1). [...] The
amplitude value at this time can be computed using the amplitudes at the
neighboring integer sample values [...] This is done by an interpolation scheme
that involves the four samples.''
\end{quotation}

This paragraph is telling us to do the calculation backwards.
Instead of mapping where each point in the CMP goes in the NMO corrected
gather, we should map where each point in the NMO gather comes from in the CMP.
Figure 1 shows a sketch of the procedure to calculate the amplitude of a point
($t_0$, $x$) in the NMO gather.

Here is the full algorithm:

\begin{enumerate}
    \item Start with an NMO gather filled with zeros.
    \item For each point ($t_0, x$) in the NMO gather, do:
    \begin{enumerate}
        \item Calculate the reflection travel-time ($t$) given $v_\mathrm{NMO}$ using
              the equation in Figure 1.
        \item Go to the trace at offset $x$ in the CMP and find the two samples
              before and the two samples after time $t$.
        \item If $t$ is greater than the recording time or if it doesn't have
              two samples after it, skip the next two steps.
        \item Use the amplitude in these four samples to interpolate the
              amplitude at time $t$.
        \item Copy the interpolated amplitude to the NMO gather at ($t_0, x$).
    \end{enumerate}
\end{enumerate}

At the end of this algorithm, we will have a fully populated NMO gather with
the amplitudes sampled from the CMP.
Notice that we didn't actually use the equation for $t_0$.
Instead we calculate the reflection travel-time ($t$).
Good luck guessing that from the equation alone.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.7\columnwidth]{figures/sketch-interpolation/sketch-interpolation}
\caption{{Figure 1: Sketch of the algorithm for a single trace and $t_0$. To the left
is a trace from the CMP and to the right the corresponding trace from the NMO
corrected gather. The green square in the NMO is the amplitude  at $t_0$ that
we want to calculate. We apply the equation to find time $t$ in the CMP, then
interpolate the amplitude using the four samples around $t$ (orange
rectangles). This amplitude is then copied over to the NMO.%
}}
\end{center}
\end{figure}

\section*{The code for NMO}

Now I'll show how to implement the above algorithm in Python using the NumPy and SciPy libraries \citep{van_der_Walt_2011}.
We'll split the algorithm into three functions.
This is very important when programming any moderately complex code because it allows us to test each part of our code independently.
It also reduces the amount of code we have to search through to find that bug that is messing up our results.
Modular code is easier to understand and to reuse.

The first function I'll define performs the NMO correction on a given CMP gather.
We'll assume that the CMP gather is a 2D array of amplitudes and that the
velocities are a 1D array with $v_\mathrm{NMO}$ for each time sample.

\begin{lstlisting}[language=python]
    import numpy as np

    def nmo_correction(cmp, dt, offsets, velocities):
        nmo = np.zeros_like(cmp)
        nsamples = cmp.shape[0]
        times = np.arange(0, nsamples*dt, dt)
        for i, t0 in enumerate(times):
            for j, x in enumerate(offsets):
                t = reflection_time(t0, x, velocities[i])
                amplitude = sample_trace(cmp[:, j], t, dt)
                if amplitude is not None:
                    nmo[i, j] = amplitude
        return nmo
\end{lstlisting}

This function is essentially the algorithm above translated to Python with some
of the details pushed into the \texttt{reflection\_time} and
\texttt{sample\_trace} functions, which we will define below.

First, the function that calculates the reflection travel-time:

\begin{lstlisting}[language=python]
    def reflection_time(t0, x, vnmo):
        t = np.sqrt(t0**2 + x**2/vnmo**2)
        return t
\end{lstlisting}

For the \texttt{sample\_trace} function, we'll use cubic splines from the
\texttt{scipy.interpolate} package.
For more information on interpolation with scipy, see the tutorial by
\citet{Hall_2016}.

\begin{lstlisting}[language=python]
    from scipy.interpolate import CubicSpline

    def sample_trace(trace, time, dt):
        before = int(np.floor(time/dt))
        N = trace.size
        samples = np.arange(before - 1, before + 3)
        if any(samples < 0) or any(samples >= N):
            amplitude = None
        else:
            times = dt*samples
            amps = trace[samples]
            interpolator = CubicSpline(times, amps)
            amplitude = interpolator(time)
        return amplitude
\end{lstlisting}

The Jupyter notebook contains the full code for these functions, including
documentation through Python documentation strings or ``docstrings'' and code
that tests that the functions work as expected.
Also included is an application of our \texttt{nmo\_correction} function to a
synthetic CMP (Figure 2).\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.7\columnwidth]{figures/nmo-application/nmo-application}
\caption{{Figure 2: (Left) The $v_\mathrm{NMO}$ profile passed to
\texttt{nmo\_correction}.
The profile was interpolated on a line using the two picked velocities (black
squares). (Middle) A synthetic CMP gather. (Right) The output
from our \texttt{nmo\_correction} function.%
}}
\end{center}
\end{figure}

\section*{Acknowledgments}

My sincerest thanks to Evan Bianco, Gregorio Kawakami, Jesper Dramsch, and Matt Hall
for comments and suggestions
and to \selectlanguage{ngerman}Öz Yilmaz for generously making the full text of the
``Seismic Data Analysis'' book available for free and in the open.

\selectlanguage{english}
\FloatBarrier
\bibliographystyle{seg}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{CJK}\end{document}

