\chapter{Introduction}

Until recently, worldwide crescent demand for hydrocarbons imposed oil \& gas geoscientists the task of continuous improvements in exploration techniques. Nowadays, strategies have diversified and pointed toward more efficient recovery and less environmental impact by better drilling planning. Of particular importance is to have robust algorithms for accurately imaging the subsurface. We typically evaluate the robustness of imaging algorithms according to their ability to obtain 1) reliable subsurface velocity/density models, and 2) accurate estimation of the subsurface reflectivity.

We achieve the first objective when the estimated velocity and density models accurately position subsurface reflectivity contrasts during the imaging process. Early estimations of velocity make use of ray tracing tomography (still the workhorse of the industry), which resorts to the wave equation's high-frequency approximation. As such, wavepaths are represented as rays (Figure \ref{fig:Ch1/wavepaths}). However, such a representation frequently fails in the presence of sharp velocity contrasts, often leaving vast ``shadows zones", i.e., areas of poor or null ray coverage. Better results can be obtained using wave-equation based velocity estimation methods, such as wave-equation migration velocity analysis (WEMVA) \citep{Biondi1999,Biondo2006} and full-waveform inversion (FWI) \citep{Tarantola1984,Virieux2009}, which perform the optimization in the image domain and the data domain, respectively. These methods exploit the band-limited character of the wavepaths that makes them sensitive to broader areas. They constitute the well-known ``banana-doughnut" sensitivity kernels \citep{Tromp2005} that reduce the occurrence of shadow zones, even in the presence of complex geology, as illustrated in Figure \ref{fig:Ch1/wavepaths}.

\multiplot{1}{Ch1/wavepaths}{width=0.95\textwidth}{Wavepaths vs. rays. Rays operate under the high-frequency approximation of the wave equation, whereas wavepaths are band-limited entities. Notice how the area of influence of the wavepaths (known as ``banana-doughnut" sensitivity kernels, because of their shape in 3D) is significantly larger than that of the rays, which often leave shadow zones in complex areas such as underneath salt bodies. Modified from \citet{Biondo2006}. \NR}

Regarding the second objective, and assuming that we fulfilled the first one, we say that the reflectivity amplitude is accurate when it truly represents lithology and fluid content heterogeneity in the subsurface \citep[e.g.][]{Ha2017}. If we achieve this objective, we can rely on amplitude variation to perform qualitative and quantitative estimation of reservoir properties.

After estimating the velocity and density subsurface models (most of the time, we regard the latter as constant, for convenience), we can confidently image the subsurface seismic events by employing conventional imaging methods such as Kirchhoff migration and wave-equation migration. Seismic migration is {\it kinematically} correct. However, it does not generally recover amplitude information that serves the purpose of interpreting the rocks and fluids contrasts. In the case of Kirchhoff migration, seismic amplitudes can be corrected using weighting coefficients derived from approximate solutions to the transport equation \citep{Schleicher2007}. For wave-equation migration, such as reverse-time migration (RTM), we can balance the amplitudes to some extent by utilizing the deconvolution imaging condition \citep{Valenciano2003}.

None of these conventional imaging methods addresses illumination problems. Likewise, they often have less than optimal seismic resolution because of the inherently blurring effect of seismic migration \citep{Yu2006}. These tasks can be achieved by means of least-squares migration (LSM) \citep{Nemeth1999}, also known as linearized waveform inversion (LWI) \citep{Tarantola1984b}. LWI can be implemented in {\it data space} \citep[e.g.][]{Nemeth1999,Duquet2000,Dai2013b,Mandy2015,Fletcher2016}, i.e., inverting for a reflectivity model, ${\bf r}$, from which data are synthesized by using the Born modeling operator, ${\bf L}$. The reflectivity model is iteratively updated until the misfit between the corresponding synthetic data, ${\bf d}({\bf b}) = {\bf L}({\bf b}){\bf r}$, and the recorded field data, ${\bf d}_{obs}$, becomes small in the least-squares sense.

We typically pose LWI in data domain as follows,

\begin{equation}
\Phi({\bf r})=\frac{1}{2} \|{\bf L}({\bf b}){\bf r}-{\bf d}_{\text{obs}} \|_2^2, \label{eq:lsqrtm}
\end{equation}

where ${\Phi}$ is the misfit function, and ${\bf b}$ is the background component of the subsurface model, either velocity or slowness or slowness squared.

LWI can also be implemented in {\it model space} \citep[e.g.][]{Valenciano2006,Valenciano.sepphd.133,Valenciano2009,Fletcher2016}. We invert for a reflectivity model that matches the migration image, ${\bf I}$, by means of the action of the Gauss-Newton Hessian of FWI, $\bf H$, as follows:

\begin{equation}
\Phi({\bf r})=\frac{1}{2} \|{\bf H}({\bf b}){\bf r}-{\bf I}({\bf b}) \|_2^2, \label{eq:lwi}
\end{equation}

where ${\bf H}={\bf L}^{\text{T}}{\bf L}$, and ${\bf L}^{\text{T}}$ represents the {\it adjoint} of the Born modeling operator, a.k.a. the RTM operator. The advantage of implementing LWI in model space in comparison with the implementation in data space is that, once we estimated the Hessian, the inversion consists of matrix-like multiplications. However, the main drawback is that the Gauss-Newton Hessian needs to be pre-computed and stored. We often resort to affordable approximations such as point-spread functions (PSF) \citep{Valenciano.sepphd.133,Fletcher2016} for the Hessian estimation.

\section{Setting the problem}

Reflectivity estimation using LWI relies upon the availability of an accurate subsurface background model. Otherwise, the image will not be well focused and the inversion can exhibit slow convergence. There have been efforts to nullify this limitation, such as extending the domain in subsurface offset \citep[e.g.][]{Elita2018}. However, remaining errors in the velocity model still make seismic reflections defocus. 
In this dissertation, I analyze the inaccuracy in the background model from a different point of view. In some cases, the background model can be accurate enough for seismic events to become reasonably well positioned after LWI. However, the accumulation of small inaccuracies in the background model can distort the reflectivity amplitude in some areas. This situation can be critical when we need such a reflectivity to estimate reservoir properties.

Figure \ref{fig:Ch1/Reservoir-a,Ch1/Reservoir-b,Ch1/Reservoir-c,Ch1/Reservoir-d} portrays a hypothetical (though entirely plausible) situation in reservoir characterization works. In some geologic environments, it makes sense to interpret amplitude anomalies to map and delimit reservoir facies, under the assumption that their constituent rocks exhibit amplitude contrast with respect to rocks of non-reservoir facies. If seismic amplitudes become altered because of various small inaccuracies in the background model, the seismic interpreter may not detect prospect opportunities. Moreover, in some cases, he or she can end up proposing drilling targets based on false amplitude anomalies, with a dry hole as the ultimate consequence.

\multiplot{4}{Ch1/Reservoir-a,Ch1/Reservoir-b,Ch1/Reservoir-c,Ch1/Reservoir-d}{width=0.45\columnwidth} {Impact of inaccuracies in the background subsurface model on the reflectivity estimation. a) Productive facies in the subsurface. b) After the seismic acquisition, processing, and reflectivity inversion, we can interpret the amplitude anomalies to map the reservoir rocks. c) Same as a), but including an unaccounted anomaly in the background. d) Reflectivity amplitude becomes altered (red traces), leading to incorrect mapping of the reservoir rocks. In more realistic scenarios, amplitudes would be affected by the combined effect of more than one anomaly. \NR}

\section{Proposed solution}

In this dissertation, I explore joint inversion of the subsurface reflectivity and background components, $\mathbf{r}$ and $\mathbf{b}$, respectively. Hereafter I will refer to such joint inversion as JIRB, standing for ``Joint Inversion for Reflectivity and Background". The idea about incorporating $\mathbf{b}$ is to account for the inaccuracies of the background model during the reflectivity inversion. Hereafter I denominate the inaccurate background model as $\mathbf{b}_0$ so that the accurate background model is given as ${\mathbf b}=\mathbf{b}_0+\mathbf{\Delta b}$, where $\mathbf{\Delta b}$ is the perturbation in the background. I incorporate the background model component into the inversion by using the WEMVA operator. The migration image, which depends on the background model, gets corrected as the inversion progresses. Simultaneously, the reflectivity is fitted to such an improved migration image using the action of the Gauss-Newton Hessian, similarly as LWI in model space operates. The reflectivity estimated with JIRB is expected to be more reliable than that obtained with conventional LWI using an inaccurate background model, ${\bf b}_0$.

\section{Opportunity area for the new method}

One drawback of the JIRB method proposed in this dissertation is the much longer computational time in comparison to conventional reflectivity inversion, mainly because of the addition of the WEMVA operator. In fact, why not performing velocity estimation separately and then going for conventional reflectivity inversion? My answer to this question is that jointly inverting for reflectivity and the background is not intended for large-scale exploration imaging, but detailed works in small areas for reservoir characterization purposes. In reconnaissance and exploration works, seismic interpreters ``comb" the subsurface in the quest of oil and gas plays and opportunities. During this stage, drilling prospects are not examined in great detail because they are not yet proved to be productive. On the contrary, once hydrocarbon accumulations have been discovered, careful and detailed interpretation is performed at the reservoir scale to propose a drilling plan for optimal exploitation. It is often required amplitude-preserving processing focused on the reservoir area. We may use LWI in model space to accomplish this task by pre-computing the Gauss-Newton Hessian and performing the inversion. However, if the background model is discovered not to be accurate enough, we may re-use the existing Gauss-Newton Hessian to run JIRB, rather than re-computing the Hessian with a previously re-estimated velocity.

There is another point of view regarding the implementation of the JIRB method. I hypothesize that the reflectivity estimation would improve if jointly performed with the background model. Conversely, the background model would benefit from the incorporation of the reflectivity in the inversion. For such a reason, in the numerical tests, I also compare the velocity estimation using WEMVA and JIRB.

\section{Notation convention}

Here I formally introduce the notation that I will use in this dissertation. Following \citet{Barnier.sep.155.gbarnier1}, let ${\bf m}$ represent the subsurface model parameters, hereafter consisting of slowness squared. Based upon the ``scale dichotomy" \citep{Shen2008}, let us assume that ${\bf m}$ can be split into the contribution of a background model (low-wavenumber component), ${\bf b}$, and a reflectivity model (high-wavenumber component), ${\bf r}$,
\begin{equation}
{\bf m}={\bf b}+{\bf r}. \label{eq:not1}
\end{equation}

As aforementioned, we can further split the background component into the superposition of an inaccurate background, ${\bf b}_0$ and perturbation effect, ${\bf \Delta b}$,

\begin{equation}
{\bf b}={\bf b}_0+{\bf \Delta b}. \label{eq:not2}
\end{equation}

Therefore, substituting equation \eqref{eq:not2} into equation \eqref{eq:not1} we obtain

\begin{equation}
{\bf m}={\bf b}_0+{\bf \Delta b}+{\bf r}. \label{eq:not3}
\end{equation}

The above assumption is fundamental in seismic imaging for techniques, including LWI and WEMVA. Hereafter I assume that the observed field data are well represented as Born linearized data, which I denote as $\bf{d}_{\text{obs}}$. Using the notation $\bf{L}(\bf{b})$ to denote the Born linearization operator evaluated at the background model, $\bf{b}$, and for the true background model and reflectivity model, ${\bf b}_{\text{true}}$ and ${\bf r}_{\text{true}}$, respectively, we obtain

\begin{equation}
    \bf d_{\text{obs}} = \bf{L}(\bf{b}_{\text{true}}) \bf r_{\text{true}}.
    \label{eq:dtrue}
\end{equation}

The migration image evaluated at the background model $\bf{b}$ will be denoted as $\bf{I}(\bf{b})$, and it is obtained as the adjoint of the Born linearization operator applied to the observed Born data $\bf{d}_{\text{obs}}$, constituting reverse-time migration,

\begin{equation}
    \bf{I}(\bf{b}) = \bf{L}(\bf{b})^{\text{T}} d_{\text{obs}}.
    \label{eq:image-def}
\end{equation}

The full-waveform inversion Gauss-Newton Hessian is similarly defined for a background model $\bf{b}$ as ${\bf{H}(\bf{b})}={\bf L(\bf{b})^{\text{T}} L(\bf{b})}$, while the wave equation migration velocity analysis operator, denoted as $\bf{W}(\bf{b})$, is defined as the derivative of the migration image with respect to the background model,

\begin{equation}
    \bf{W}(\bf{b}) = \frac{\partial \bf{I}(\bf{b})}{\partial \bf{b}}.
    \label{eq:wemva}
\end{equation}

%The tomographic operator, denoted as $\bf{T}(\bf{b})$, constitutes the derivative of the Born data with respect to the background model \citep{Almomin.sep.149.ali2, Barnier.sep.155.gbarnier1},

%\begin{equation}
%    \bf{T}(\bf{b}) = \frac{\partial \bf{d}}{\partial \bf{b}} (\bf{b}) = \frac{\partial \bf{L}}{\partial \bf{b}} (\bf{b},\bf{r}).
%    \label{eq:tomo}
%\end{equation}

%Note that the derivative of the Born modeling operator with respect to the background model is not a matrix operator by itself. It is its action upon the reflectivity that becomes a matrix operator.


\section{Thesis overview}

\subsection{Chapter 2: Theory}

In this chapter, I formally derive the JIRB method. I discuss the initial stage, where I had set the problem as a linear optimization process. Then I show that the linearization of the migration image with respect to the background model makes this approach fail. Finally, I show the current nonlinear implementation.

An additional term is needed to drive the background model to solutions that correct and focus the migration image. In the zero subsurface offset domain, we can enforce this requirement by maximizing the energy of the migration image \citep{Tang.sep.145.yaxun1}. This method is similar to semblance maximization \citep{Soubaras2007} and stack power maximization \citep{Zhang2013}. There are two choices in the extended domain: maximization of the stacking power (for angle-domain common image gathers) or minimization of the image after applying differential semblance optimization (for offset-domain common image gathers).

\subsection{Chapter 3: Random boundary condition}

The random boundary condition (RBC) was originally proposed by \citet{Clapp2009} for RTM, to avoid the storage of the whole source wavefield time history. It can also be used for WEMVA to avoid the storage of source and receiver wavefields time histories, enabling the application of the method under limited-memory computational settings, at the cost of additional wavefields propagations. In this chapter, I elaborate on the implementation of the RBC as part of the JIRB constituents, in particular for the WEMVA operator.

\subsection{Chapter 4: Application to synthetic 2D data}

In this chapter, I show the method's implementation in the sedimentary section of the Sigsbee model. I synthesized the data using the acoustic scalar Born modeling operator. I prepared five inversions. The first one is conventional LWI using as input the migration image with the wrong background. The second one is conventional LWI using the migration image with the correct background model. The third one is the JIRB test using the wrong background model. I evaluated the effect of the wrong background model in the estimated reflectivity, compared with the correct result, and to which extent JIRB is capable of estimating the correct background while correcting the reflectivity during the inversion. The fourth and fifth inversions consist of WEMVA and JIRB, respectively, whose background estimations are followed by conventional LWI for comparison purposes.

\subsection{Chapter 5: 3D Application to 3D ocean-bottom node data from the Gulf of Mexico}

In the last chapter, I apply the JIRB method to a 3D ocean-bottom node dataset provided by Shell. I perform the imaging experiments using the downgoing component as input field data, which has been previously obtained using the PZ-summation technique. I employ the well-known technique of mirror imaging. I initially planned on performing two inversions. The first one is conventional LWI using the migration image with the available velocity model. The second one is the JIRB method. The objective is to evaluate the ability to focus residual inaccuracies in the estimated reflectivity. However, the estimated reflectivity using JIRB significantly differs from that of LWI, rather than just constituting a better version of the latter. Therefore, I added numerical tests comparing the background models obtained with WEMVA and JIRB in RTM volumes produced with refined grids to better appreciate changes in the stratigraphy.