\chapter{Application to 3D marine data}

In this chapter, I illustrate the JIRB method by presenting an application to an ocean bottom node (OBN) 3D dataset from the Gulf of Mexico, provided by Shell. I followed similar steps as in the synthetic tests of the previous chapter. With the tests shown in the current chapter, I evaluate the plausibility and performance of the JIRB method in real data, free of the inverse crime.

For comparison purposes, I prepared conventional LWI reflectivity. Similar to the previous chapter, I implemented JIRB and analyzed the corresponding background and reflectivity components. I also prepared an LWI test using the JIRB background model. Finally, I further evaluate the performance of JIRB versus WEMVA by utilizing the corresponding background models as input for fine-scale RTM.

I judged the results according to two criteria: 1) image focusing, and 2) enhancement of geological features. The seismic interpretation resorts to the limited available knowledge of the geology of the area. More accurate interpretations require more regional assessment of the geology, and incorporation of borehole data.

\section{Preliminary steps}

\subsection{Data preparation}

Figure \ref{fig:allshots,allrecs} shows the areal distribution of the shots and the receivers (nodes) for the OBN acquisition campaign. The target consists of sedimentary reservoir rocks pinching out against the wall of a salt diapir, with apparent episodes of passive and active growths \citep{Hudec2017}. Stratigraphic traps are also present, in particular related to deepwater deposits. The squares in both maps encompass the imaging area that I chose for the numerical tests. It contains part of the sedimentary volume and the diapir. There are $226$ nodes contributing to the image within this area. The original recording time is $14$ s sampled at $2$ ms. For the numerical tests shown in this chapter, I only employed the first $8$ s.

\multiplot[t]{2}{allshots,allrecs}{trim=0.1cm 3cm 0.1cm 3cm, clip=true,width=0.47\textwidth}{Areal distribution of a) shots, b) receivers (nodes). The small squares indicate the image area for the numerical tests. \CR}

\multiplot[t]{2}{wavelet,wavelet-spec}{width=0.47\textwidth}{a) Source wavelet and b) the corresponding amplitude spectrum. \ER}

\citet{Dahlke.sepphd.175} performed the PZ-summation to separate the upgoing and the downgoing wavefield components. This procedure consists of the adaptive addition and subtraction of the geophone vertical component and the hydrophone component of the data \citep{Soubaras1996,Pacal2015}. \citet{Dahlke.sepphd.175} also reshaped the wavelet and removed the bubble effect. I utilized these pre-processed data for the numerical tests shown in this chapter. Although the upgoing component produces sharper reflectivity contrasts, I chose the downgoing component because it has better coverage of the subsurface, and therefore, less tendency to produce aliasing artifacts in the images \citep{Pacal2015}. I employed the standard mirror technique \citep{Grion2007,Pacal2015} for the imaging processes. I also resorted to the principle of reciprocity \citep{Knopoff1959}, injecting the source wavelet at the nodes' positions, whereas I injected the data traces at the source positions in the corresponding common-receiver gathers (CRG).

I binned the downgoing component data (which I henceforth refer to simply as ``the data") to obtain a regular grid of $25\times25$ m, which is suitable for imaging frequencies up to $30$ Hz. For simplicity, I constructed the binned CRGs spanning the whole study area. I employed a $10.5$ Hz dominant frequency Ricker wavelet as a source function, bandpassed to a maximum frequency of $30$ Hz. Figure \ref{fig:wavelet,wavelet-spec} shows the normalized wavelet and the corresponding amplitude spectrum.

\subsection{Reverse-time migration}

I prepared conventional RTM volumes to have the first view of the subsurface in the study area with my computational code [RTM images of the area are also available in \citet{Dahlke.sepphd.175}]. It allows the evaluation of my imaging toolbox on this OBN dataset. The RTM image also constitutes the input data for model-space LWI and JIRB.

I prepared the velocity model to image using the downgoing component, as shown in Figure \ref{fig:mirror1,mirror2,mirror3}. First, I numerically expanded the water layer to enact the mirror imaging strategy. Then, I defined a datum corresponding to a numerical ocean floor mirrored about the water surface, and where the mirrored nodes will lie. Finally, as aforementioned, using reciprocity, I inject the source wavelet at the position of such mirrored nodes and the data at the source positions.

\multiplot[t]{3}{mirror1,mirror2,mirror3}{width=0.75\textwidth}{a) Ocean-bottom node acquisition cartoon, illustrating the water-surface multiple trajectory for a source-node pair. b) Numerical extension of the water velocity and re-datuming of the nodes, showing the equivalent source-node trajectory. c) After the application of the reciprocity principle to a source-node pair. \NR}

The migration aperture is $1250$m ($50$ samples) in both $x$ and $y$ directions. For the random boundary condition \citep{Clapp2009} I surrounded the volume with a halo of $50$ samples. Like the synthetic data in Chapter 4, I included an extra layer of $50$ samples on the top boundary to avoid placing the ``sources" (which are the nodes acting as such, for reciprocity) upon the random layer's border, thereby preventing numerical artifacts. Figure \ref{fig:RBC-vel} shows one velocity volume prepared for the migration of a single CRG. The lines indicate planes whereby the inline, the crossline, and the depth slice sections project to the lateral sides and top of the volume. Note that the randomness is mild at the volumes' interior, and increases toward the computational boundaries. I pre-computed $226$ velocity volumes for different realizations of the engine code to generate different random halos. This strategy allows one to cancel out the incoherent random artifacts while reinforcing the subsurface signal in the image during stacking.

\multiplot[t]{0}{RBC-vel}{width=0.99\textwidth}{Random boundary frame for a single common-receiver gather. \CR}

Figure \ref{fig:rtm3d} shows inline, crossline, and depth slice sections of the RTM volume. I applied a Laplacian filter to attenuate the low-wavenumber artifacts. The reflectors are almost flat except in the vicinity of the salt intrusion. The sedimentary units are tilted by the diapir growth, which is still deforming the ocean floor. I interpreted the presence of a possible channel complex at the indicated position. The presence of this geological element can be expected in turbidite systems \citep[e.g.][]{Huang2018}. See Figure \ref{fig:turbidite}. The diapir's walls are virtually vertical and are barely imaged with the chosen aperture. In other tests (not shown), I recovered such a diapir wall by using a larger migration aperture, but at the expense of introducing migration ``smiles" that contaminated the reflectors. The purpose of the inversion tests is to properly image the sedimentary structures, rather than imaging the salt flanks. Hence, I keep the migration aperture small.

\multiplot[t]{0}{turbidite}{width=\textwidth}{ Illustration of turbidite system. Taken from \citet{Huang2018}.  \NR}

For the discussion on the image differences and improvements, I mainly focused on sedimentary features, such as those indicated by arrows in Figure \ref{fig:rtm3d-arrows}, which probably correspond to a confined channel complex \citep[e.g.][]{Huang2018}.

\multiplot[t]{2}{rtm3d,rtm3d-arrows}{width=0.75\textwidth}{a) RTM image of the OBN data. b) Same as a), with arrows indicating the base and the lateral bounds of a possible channel complex. \CR}

%\plot[t]{rtm3d-arrows}{width=0.97\textwidth}{Same as Figure \ref{fig:rtm3d-rot}. The arrows indicate the base and the lateral bounds of a possible channel complex. \NR}

\FloatBarrier

\subsection{Point-spread functions}

Like the synthetic test of the previous chapter, I employed PSFs to estimate the Gauss-Newton Hessian for the 3D data tests. I seeded spikes every $15$ gridpoint and cascaded Born modeling followed by RTM to compute such PSFs, shown in Figure \ref{fig:psf3d}. From this volume, the PSFs are read and interpolated ``on the fly" to estimate Hessian terms. Notice that at the salt diapir's position, the energy of the PSFs becomes more smeared and faint than in the sedimentary part. Although Figure \ref{fig:psf3d} shows the distribution of PSFs only in the model space, they span all across the computational space (i.e., including the boundaries halo).

\multiplot{0}{psf3d}{width=0.9\textwidth}{PSF distribution in the image space. \CR}

%\FloatBarrier

\section{Inversion results}

\subsection{Reflectivity comparison}

In this section, I compare the result of conventional LWI and JIRB for the reflectivity estimation. It is homologous to the numerical test in section \ref{ssec:jirb-vs-lwi}. 

Figure \ref{fig:hess3d} shows the result of conventional LWI in model space, i.e., the result of solving equation \eqref{eq:lwi}. The background model is the original provided by Shell and is not updated during the LWI process. The Gauss-Newton Hessian was approximated by an ``on the fly" interpolation of the PSFs shown in Figure \ref{fig:psf3d}. The inversion ran for 20 iterations, which was enough to flatten the objective function. The inversion recovered some high wavenumber components. It also introduced vertical artifacts in the estimated reflectivity image, even destroying the top of the salt diapir.

Figure \ref{fig:jirb3d-comp1} shows the result of JIRB reflectivity after $10$ iterations, i.e., the result of solving equation \eqref{eq:jirb_final_prec}. Notice that the events of the JIRB image are more focused. Such is a consequence of updating the background component at the same time as the reflectivity. However, this result is different from merely an improved version of the LWI result (Figure \ref{fig:hess3d}). Logically, there is a more complex interaction between the model parameters (the reflectivity and background) than what I expected when I initially formulated the method as a linear problem. 

%The result, however, is significantly different from what I expected when I had formulated the problem as a linear one (equation \ref{eq:orig_lwivu2}), i.e., a reflectivity estimation similar to Figure \ref{fig:hess3d}, but with more focused seismic events. The reason for the difference could be that the nonlinear interaction between the model parameters is more complex than the simplistic, linear formulation, had the latter worked.

\multiplot{2}{hess3d,jirb3d-comp1}{width=0.75\textwidth}{a) Reflectivity obtained with conventional LWI. b) Reflectivity component of JIRB. \CR}

%\plot[t]{hess3d-rot}{width=0.99\textwidth}{a) Reflectivity obtained with conventional LWI. \CR}

%\plot[t]{jirb3d-comp1-rot}{width=0.96\textwidth}{Reflectivity component obtained with JIRB. \CR}

\FloatBarrier

\subsection{Background component comparison}

\subsubsection{Background component inversions}

Figures \ref{fig:slow3d}, \ref{fig:wemva3d}, and \ref{fig:jirb3d-comp2} show the initial or original background model, the inverted background model using WEMVA, and the inverted background model using JIRB for $\lambda=8$, respectively. During the inversions, I applied a mask to the gradient to prevent updating the squared slowness of the salt body and the water layer (see equation \eqref{eq:grad-wemva-prec-w}). Notice that the initial model appears to be laterally smooth and relatively featureless in the sedimentary part, whereas the updated models exhibit vertical and lateral variations. This is better observed in the difference between the inverted and the initial background models in Figures \ref{fig:wemva-slow-dB3d} and \ref{fig:jirb-slow-dB3d}. The additional detail obtained with the inversions should contribute to focus the image better. We observe that improvement in the estimated JIRB reflectivity image (Figure \ref{fig:jirb3d-comp1}). On the other hand, geologically speaking, it makes perfect sense to expect such heterogeneity in a deepwater depositional environment. In both inversions, JIRB and WEMVA, I preconditioned the background model using B-splines (see Appendix \ref{appendix:prec}), which helped prevent the occurrence of spurious, high wavenumber, artifacts.

\multiplot[t]{0}{slow3d}{width=\textwidth}{Original background model. ER}

\multiplot[t]{0}{wemva3d}{width=\textwidth}{Background model obtained with WEMVA. ER}

\multiplot[t]{0}{jirb3d-comp2}{width=\textwidth}{Background model obtained with JIRB. \CR}

\multiplot[t]{0}{wemva-slow-dB3d}{width=\textwidth}{Difference between WEMVA and initial background models. \CR}

\multiplot[t]{0}{jirb-slow-dB3d}{width=\textwidth}{Difference between JIRB and initial background models. \CR}

\FloatBarrier

\subsubsection{Refined RTM images}

It remains in a quandary which one of the inversions achieved better result, WEMVA, or JIRB. It is very difficult to make a sensible judgement based solely on the results of Figures \ref{fig:slow3d} to \ref{fig:jirb-slow-dB3d}. One alternative is to proceed similarly as in Chapter 4: performing LWI using the corresponding inverted background models. However, for this 3D case, I performed refined RTM tests using the inverted models, reducing the grid size to $12.5\times12.5\times12.5$ m, i.e., half the original grid size and thereby duplicating the wavenumber of the migration image. This refinement allows for better imaging of the stratigraphic features that I want to enhance. The drawback of such refinement is the increase in computational time and storage by a factor of 8 or 16 in 3D.

Figure \ref{fig:rtm-orig-fine} shows the refined RTM image obtained with the original background model (Figure \ref{fig:slow3d-rot}). I zoomed into the sedimentary section to better appreciate the stratigraphy, excluding the water layer and the salt diapir. I show the same inline and crossline that I showed before in the previous results, but a shallower depth slice (approximately 1100 m). The latter exhibits channel signatures, which are very common in deep water environments. For comparison, Figure \ref{fig:rtm-wemva-fine} shows the RTM refined image using the WEMVA inverted background model (Figure \ref{fig:wemva3d-rot}). Notice how the stratigraphic features of the inline and crossline sections become more focused in the WEMVA volume, e.g., the channel features immediately below the reflector at 2000 m. However, the comparison is not valid for the depth slice because the velocity change introduced by the inversion significantly shifted the seismic events upward. Thus, while in the depth slice of Figure \ref{fig:rtm-orig-fine} we interpret channels, in the depth slice of \ref{fig:rtm-wemva-fine} we can interpret other geological elements, a fan constituted by the white feature at the top of the depth slice, fed by a channel, and either crevasse splays or minor fan lobes originating from the same channel.

\multiplot{2}{rtm-orig-fine,rtm-wemva-fine}{width=0.8\textwidth}{Refined RTM volume for a) Original model, b) WEMVA model. \CR}

In the following, I compare the WEMVA- and the JIRB-based RTM images. Their differences are much more subtle than those with the RTM for the original background model. I highlighted some of such differences; others are difficult to observe without toggling the images.

Figures \ref{fig:rtm-wemva-fine-1-a} and \ref{fig:rtm-jirb-fine-1-a} show the refined RTM image obtained with the WEMVA and the JIRB background models, respectively (Figures \ref{fig:wemva3d-rot} and \ref{fig:jirb3d-comp2-rot}). The depth slice shows improvement in a part of the fan and one of the channels, as indicated by the red ellipses and the red arrow. In the sections, the yellow ellipses indicate differences in the reflectors that suggest more focusing. There are also vertical adjustments of the seismic events due to the velocity discrepancy between the background models.

Figures \ref{fig:rtm-wemva-fine-2-a} and \ref{fig:rtm-jirb-fine-2-a} show the comparison for a deeper level (around 1300 m). The left half of the depth slice is featureless, possibly an abyssal plain. The right half shows a straight channel that bifurcates into two fans with a dendritic pattern. The rightmost feature appears to be another fan lobe, with internal channelization. The main improvements in the JIRB volume that I can identify are the better definition of the channels inside the first lobe (red oval) and the sharper border of the second one, both indicated by red ellipses. On the contrary, the straight channel, indicated by the red arrows, is better defined in the WEMVA volume. In the section, I highlighted some channel features that appear to be better focused when using JIRB.

Figure \ref{fig:rtm-wemva-fine-3a-a,rtm-jirb-fine-3a-a} shows a deeper level, around $1800$m, where some of the changes are indicated in red. It corresponds to the top of a stratigraphic unit between $1800$ and $2000$m that I interpret as a possible mass-transport complex that spans beyond the limits of the study area. In the vertical sections, we observe shingle geometries (short reflector units beginning in toplaps and ending in downlaps), and possible channels and fans deposited on top. Although there are features in the sections that noticeably become more focused in the JIRB volume (indicated by the yellow ellipses), the improvements are much more difficult to observe in the depth slice.

Figures \ref{fig:rtm-wemva-fine-4-a} and \ref{fig:rtm-jirb-fine-4-a} show a channelized level at 2000 m. The channels' signatures are distinguishable in the sections. In the depth slice, we observe that they constitute part of the channel complex in white upon the dark gray plain, although some other, grayish, channels, are visible at the lower part of the depth slice. The red ellipses indicate features that are better observed in the JIRB volume. Most of them are obscured in the WEMVA volume because the corresponding background model was either inaccurately focus them, and/or because they are slightly shifted with respect to the depth slice. The correction in the JIRB volume brought them back to the stratigraphic level. Additional features that are slightly corrected by JIRB are indicated in the sections by the yellow ellipses and arrow.

\multiplot{2}{rtm-wemva-fine-1-a,rtm-jirb-fine-1-a}{width=0.8\textwidth}{Refined RTM volume for a) WEMVA model, b) JIRB model. \CR}

\multiplot{2}{rtm-wemva-fine-2-a,rtm-jirb-fine-2-a}{width=0.8\textwidth}{Refined RTM volume for a) WEMVA model, b) JIRB model. \CR}

\multiplot{2}{rtm-wemva-fine-3a-a,rtm-jirb-fine-3a-a}{width=0.8\textwidth}{Refined RTM volume for a) WEMVA model, b) JIRB model. \CR}

\multiplot{2}{rtm-wemva-fine-4-a,rtm-jirb-fine-4-a}{width=0.8\textwidth}{Refined RTM volume for a) WEMVA model, b) JIRB model. \CR}

\section{Conclusions}

The JIRB method results on the OBN 3D dataset show some subtle improvements in comparison with the WEMVA results. It demonstrates that the joint inversion of the reflectivity and the background model by combining the least-squares RTM and the WEMVA operators can potentially correct residual inaccuracies in velocity, which in turn benefits the estimated reflectivity. The JIRB method yielded a better estimation of the reflectivity than LWI, but at a much higher cost. I had initially sought the improvement of the reflectivity by allowing updating the background model. However, as my project progressed, I also had to look at the method from a different perspective: to improve the background model by incorporating the reflectivity, instead of inverting for the background model alone. As so, I compared the JIRB vs. the WEMVA inversion, which are much more comparable in cost. There are small yet noticeable improvements with the JIRB method. Such improvements give a firm clue that incorporating the reflectivity into the background model inversion contributes to the correction of small inaccuracies in the latter. There is still plenty of room for future research on this topic. Of particular interest is obtaining an appropriate value of the $\lambda$ parameter, which in this dissertation, I had to estimate by trial and error.


