\chapter{Theory}

In this chapter, I present formal derivations of different implementations the JIRB algorithm. I first show the linearization from the most general case and the reason for its failure. Next, I present the original derivation from the FWI objective function. Next, I show the final nonlinear JIRB objective function. Finally, I present some implementation ideas in the extended domain. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Derivation of JIRB: Most general cases}
Let us recall equation \eqref{eq:lwi} for conventional LWI in model space:

\begin{equation}
\Phi({\bf r})=\frac{1}{2} \left\Vert{\bf H}({\bf b}){\bf r}-{\bf I}({\bf b}) \right\Vert_2^2. \label{eq:lwi_conv}
\end{equation}

The most general case of JIRB simply consist in the incorporation of ${\bf b}$ into the inversion:

\begin{equation}
\Phi({\bf r},{\bf b})=\frac{1}{2} \left\Vert{\bf H}({\bf b}){\bf r}-{\bf I}({\bf b}) \right\Vert_2^2. \label{eq:jirb_gen}
\end{equation}

Note that both the Hessian and the migration image depend on the background model ${\bf b}$, which hereafter will correspond to the {\it accurate} background model that we do not know a priori.

In equation \eqref{eq:jirb_gen}, the model space becomes expanded, and so it is the model null space \citep{Aster2013}. Therefore, it becomes necessary to incorporate another restriction to drive ${\bf b}$ away from spurious solutions. For zero subsurface offset imaging, we can maximize the energy of the migration image to constrain ${\bf b}$, obtaining

\begin{equation}
\Phi({\bf r},{\bf b})=\frac{1}{2} \left\Vert{\bf H}({\bf b}){\bf r}-{\bf I}({\bf b}) \right\Vert_2^2-\frac{\lambda}{2} \left\Vert {\bf I}({\bf b})\right\Vert_2^2, \label{eq:jinv-b}
\end{equation}

where $\lambda$ is a trade-off parameter. The minus sign allows maximizing the second functional within a minimization scheme. The interpretation of equation \eqref{eq:jinv-b} is that the reflectivity fits an updated migration image that gets improved as the background model is corrected. Now I substitute equation \eqref{eq:not2} into \eqref{eq:jinv-b} to obtain

\begin{equation}
\Phi({\bf r},{\bf b}_0+{\bf \Delta b})=\frac{1}{2} \left\Vert{\bf H}({\bf b}_0+{\bf \Delta b}){\bf r}-{\bf I}({\bf b}_0+{\bf \Delta b}) \right\Vert_2^2-\frac{\lambda}{2} \left\Vert {\bf I}({\bf b}_0+{\bf \Delta b})\right\Vert_2^2, \label{eq:jinv-db}
\end{equation}

where hereafter ${\bf b}_0$ will represent the {\it inaccurate} background model, which we assume it is available. Now it is the perturbation ${\bf \Delta b}$ what becomes a model parameter in the inversion. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Linearization attempt}

Equation \eqref{eq:jinv-db} is nonlinear with respect to ${\bf \Delta b}$ in both the Gauss-Newton Hessian and the migration image. We can make use of the following expansions:

\begin{subequations}
    \begin{eqnarray}
    {\bf H}({\bf b}_0+{\bf \Delta b}){\bf r}&=&{\bf H}({\bf b}_0){\bf r}+O({\bf \Delta b},{\bf r}), \label{eq:hess-approx} \\
    {\bf I}({\bf b}_0+{\bf \Delta b})&=&{\bf I}({\bf b}_0)+{\bf W}({\bf b}_0){\bf \Delta b}+O(\|{\bf \Delta b}\|^2). \label{eq:image-approx}
    \end{eqnarray}
\end{subequations}

For a sufficiently small ${\bf \Delta b}$ we can drop the higher order terms, ($O({\bf \Delta b},{\bf r})$ and $O(\|{\bf \Delta b}\|)$). After doing that, equation \eqref{eq:hess-approx} means that the action of the Hessian upon the reflectivity when it is evaluated at the correct background, is not significantly different from the same action of the Hessian, but evaluated at the incorrect background. In this case the dropped higher order terms are second order in ${\bf r}$ and ${\bf \Delta b}$ simultaneously, which are both small with respect to ${\bf b}_0$. To support this claim I performed numerical experiments with the sedimentary section of the Sigsbee 2A model. Figure \ref{fig:Ch2/Sigs-svmodel,Ch2/Sigs-gauss-vel,Ch2/Sigs-refl} show the true velocity, a perturbation in velocity, and the true reflectity. Figures \ref{fig:Ch2/hess-b-200,Ch2/hess-b0-200,Ch2/hess-diff-200}, \ref{fig:Ch2/hess-b-800,Ch2/hess-b0-800,Ch2/hess-diff-800}, and \ref{fig:Ch2/hess-b-1600,Ch2/hess-b0-1600,Ch2/hess-diff-1600} compare the terms of expression \eqref{eq:hess-approx} for Gaussian anomalies similar to that of Figure \ref{fig:Ch2/Sigs-gauss-vel}, but with peak values of -200, -800, and -1600 ft/s for the Gaussian anomaly. Notice that the difference becomes evident only at 1600 ft/s. This value represents an error of approximately 20\% with respect to the background average, which is larger than the residual anomalies intended for the either the WEMVA or the JIRB methods.

\multiplot{3}{Ch2/Sigs-svmodel,Ch2/Sigs-gauss-vel,Ch2/Sigs-refl}{width=0.47\columnwidth}{Sedimentary section of the Sigsbee model, showing: a) True velocity model. b) Gaussian anomaly. c) True reflectivity. \CR}

\multiplot{3}{Ch2/hess-b-200,Ch2/hess-b0-200,Ch2/hess-diff-200}{width=0.47\columnwidth}{Test on the validity of expression \eqref{eq:hess-approx} for the Sigsbee model sedimentary section, using a Gaussian velocity perturbation of 200 ft/s: a) ${\bf H}({\bf b}){\bf r}$; b) ${\bf H}({\bf b}_0){\bf r}$. c) Difference between images a) and b). \CR}

\multiplot{3}{Ch2/hess-b-800,Ch2/hess-b0-800,Ch2/hess-diff-800}{width=0.47\columnwidth}{Test on the validity of expression \eqref{eq:hess-approx} for the Sigsbee model sedimentary section, using a Gaussian velocity perturbation of 800 ft/s: a) ${\bf H}({\bf b}){\bf r}$; b) ${\bf H}({\bf b}_0){\bf r}$. c) Difference between images a) and b). \CR}

\multiplot{3}{Ch2/hess-b-1600,Ch2/hess-b0-1600,Ch2/hess-diff-1600}{width=0.47\columnwidth}{Test on the validity of expression \eqref{eq:hess-approx} for the Sigsbee model sedimentary section, using a Gaussian velocity perturbation of 1600 ft/s: a) ${\bf H}({\bf b}){\bf r}$; b) ${\bf H}({\bf b}_0){\bf r}$. c) Difference between images a) and b). \ER}

Expression \eqref{eq:image-approx} constitutes the Taylor's series expansion of the migration image around the wrong background model, which we can truncate after the first-order term in ${\bf \Delta b}$. ${\bf W}$ represents the forward WEMVA operator, which is the derivative of the migration image with respect to the background model. This expression is valid for perturbation values such that $\|{\bf \Delta b}\| \ll \|{\bf b}_0\|$.

Substituting equations \eqref{eq:hess-approx} and \eqref{eq:image-approx} into equation \eqref{eq:jinv-db} yields

\begin{equation}
\Phi({\bf r},{\bf \Delta b})=\frac{1}{2} \left\Vert{\bf H}({\bf b}_0){\bf r}-{\bf I}({\bf b}_0)-{\bf W}({\bf b}_0){\bf \Delta b} \right\Vert_2^2-\frac{\lambda}{2} \left\Vert {\bf I}({\bf b}_0)+{\bf W}({\bf b}_0){\bf \Delta b}\right\Vert_2^2. \label{eq:lwivu-zeroh}
\end{equation}

The first functional in equation \eqref{eq:lwivu-zeroh} has the expected behavior of minimizing the difference between the reflectivity under the action of the Hessian, and the migration image evaluated at the wrong background, corrected with the WEMVA term applied to the perturbation in the background. This term becomes a perturbation in the image that contains the corrections to the background model.

However, the second functional in equation \eqref{eq:lwivu-zeroh} does not have the expected behavior. The reason is because, whereas the linearization in expression \eqref{eq:image-approx} is valid for small ${\bf \Delta b}$ values, substituting such a linearized term into the norm gives

\begin{equation}
\left\Vert{\bf I}({\bf b}_0+{\bf \Delta b})\right\Vert_2^2 \approx \left\Vert{\bf I}({\bf b}_0)+{\bf W}({\bf b}_0){\bf \Delta b}\right\Vert_2^2, \label{eq:norm-image-approx}
\end{equation}

which is inaccurate. In fact, expression \eqref{eq:norm-image-approx} preserves a second order term in ${\bf \Delta b}$, but neglects another (see Appendix \ref{appendix:lin} for details). The consequence is that the linearized norm in the right-hand side of expression \eqref{eq:norm-image-approx} behaves as a quadratic function of the type $Ax^2+Bx+C$, with curvature term $A\geq0$ because it is represented by a positive semidefinite matrix \citep{Nocedal2006,Aster2013}. In other words, the linearized norm is {\it convex}, thus there is no maximum at the true background, as it is expected from the original functional (left-hand side of expression \eqref{eq:norm-image-approx}), which is {\it concave}.

To further illustrate the problem, I performed numerical tests on a flat-layer model and the sedimentary sector of the Sigsbee model (Figure \ref{fig:Ch2/Flat-svmodel,Ch2/Sigs-svmodel}). The flat-layer model includes a positive Gaussian anomaly in velocity. This model represents the true background model, ${\bf b}^{\text{true}}={\bf b}_0+{\bf \Delta b}^{\text{true}}$. The incorrect background model ${\bf b}_0$ does not contain the anomaly. Hence, this anomaly represents the true perturbation in the background, ${\bf \Delta b}^{\text{true}}$. On the other hand, the true Sigsbee background model does not contain a Gaussian anomaly. However, to construct an incorrect background model, I added a negative Gaussian anomaly in velocity to the true background model. The Gaussian anomalies for both models are shown in Figure \ref{fig:Ch2/Flat-gauss-vel,Ch2/Sigs-gauss-vel-cp}.

I synthesized true data $\mathbf{d}^{\text{true}}$ in both cases by applying Born modeling to the corresponding reflectivity models (Figure \ref{fig:Ch2/Flat-refl,Ch2/Sigs-refl}) using the true background models. For the flat-layer model, the acquisition geometry consists of $101$ shots regularly spaced every $140$ m, and receivers every $20$ m. For the Sigsbee model, the acquisition geometry consists of 54 shots spaced every $500$ ft, and receivers every $75$ ft.

The first test demonstrates that the image power maximization term in the linearized objective function (right-hand side of expression \eqref{eq:image-approx}) does not have the desired property of concavity. To do this, I created trial perturbations in the background model, given as ${\bf \Delta b}^{\text{true}}+\alpha{\bf \Delta b}_{\text{rand}}$, where ${\bf \Delta b}^{\text{true}}$ is the true anomaly, ${\bf \Delta b}_{\text{rand}}$ represents randomly distributed perturbations whose amplitudes are confined within the interval $\big[-\text{max}({\bf \Delta b}^{\text{true}}),\text{max}({\bf \Delta b}^{\text{true}})\big]$, and ${\alpha}$ is a scalar parameter with values $-1,-0.9,...,0.9,1$. Then I evaluated the functional, $\left \Vert{\bf I}({\bf b}_0)+ {\bf W}({\bf b}_0)\big[{\bf \Delta b}^{\text{true}} + \alpha{\bf \Delta b}_{\text{rand}} \big]\right \Vert_2^2$, using different random realizations of ${\bf \Delta b}_{\text{rand}}$, and exploring along each direction with the $\alpha$ values, so $\alpha=0$ corresponds to the evaluation of the objective function using the true background anomaly.

\multiplot{2}{Ch2/Flat-svmodel,Ch2/Sigs-svmodel}{width=0.47\columnwidth,trim=0in 0in 0in 0in, clip=true} {True background models in velocity. a) Flat model. b) Sigsbee model. \ER}

\multiplot{2}{Ch2/Flat-svmodel-pert,Ch2/Sigs-svmodel-pert}{width=0.47\columnwidth,trim=0in 0in 0in 0in, clip=true} {Wrong background models in velocity. a) Flat model. b) Sigsbee model. \ER}

\multiplot{2}{Ch2/Flat-gauss-vel,Ch2/Sigs-gauss-vel-cp}{width=0.47\columnwidth,trim=0in 0in 0in 0in, clip=true} {Perturbations in the background models expresed in velocity. a) for the flat model. b) For the Sigsbee model. \ER}

\multiplot{2}{Ch2/Flat-refl,Ch2/Sigs-refl}{width=0.47\columnwidth,trim=0in 0in 0in 0in, clip=true} {Reflectivity models in slowness squared. a) Flat model. b) Sigsbee model. \ER}

%\subsection{Testing the migration image power term in the linearized objective function}

Figure \ref{fig:Ch2/Flat-test-wemva-alpha,Ch2/SigsA-test-wemva-alpha} shows the result of 50 tests of the evaluation of the objective function using random directions, ${\bf \Delta b}_{\text{rand}}$, each one for the 21 values of $\alpha$. Notice that there is no upper bound, and although the minima of the curves cluster around the vicinity of the true anomaly, none of them corresponds to it. These results prove that this objective function is convex, as I discussed before.

\multiplot{2}{Ch2/Flat-test-wemva-alpha,Ch2/SigsA-test-wemva-alpha}{width=0.47\textwidth}{Experiment results of 50 random directions of the $\|{\bf I}({\bf b}_0)+{\bf W}({\bf b}_0){\bf \Delta b}^{\text{true}}+\alpha {\bf \Delta b}_{\text{rand}}\|_2^2$ functional in the a) flat model, and the b) Sigsbee model. Each random direction was explored with the $\alpha$ scalar. \NR}

\multiplot{2}{Ch2/Flat-test-rtm-alpha,Ch2/SigsA-test-rtm-alpha}{width=0.47\textwidth}{Experiment results of 50 random directions of the $\|{\bf I}({\bf b}^{\text{true}}+\alpha {\bf \Delta b}_{\text{rand}})\|_2^2$ functional in the a) flat model, and the b) Sigsbee model. Each random direction was explored with the $\alpha$ scalar. \NR}
%\subsection{Test on the fully nonlinear migration image power term}

I performed similar tests for the fully nonlinear objective function of image power maximization (left-hand side of expression \eqref{eq:image-approx}). Similar to the previous test, I evaluated this function at random directions. The results for both models are shown in Figure \ref{fig:Ch2/Flat-test-rtm-alpha,Ch2/SigsA-test-rtm-alpha}. These results confirm that this objective function is concave along the directions explored. While these experiments do not tell us for sure if it is, in fact, concave, we can still conclude that this may be an approximate heuristic. However, even though the maximum values of the curves cluster around $\alpha=0$, not all of them are exactly there.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Derivation from full-waveform inversion}
The first attempt to envision a JIRB method was through the simplification of the Newton's equation associated to the FWI iterative solution \citep[e.g.][]{Cabrales-Vargas.sep.165.alejandro1,Cabrales-Vargas.sep.172.alejandro1}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Full-waveform inversion \citep{Tarantola1984,Pratt1998,Virieux2009} is a nonlinear optimization technique that, in principle, can estimate all the wavenumber components of the subsurface model parameters, ${\bf m}$. However, in practice, it suffers from cycle skipping if the difference between the initial model and the actual model is more than half a wavelength in the data space. Multiscaling approaches are often employed to circumvent the local minima problem \citep[e.g.][]{Virieux2009}.

In a nutshell, FWI aims to minimize in the least-squares sense the misfit between the observed seismic data and synthetic data corresponding to model ${\bf m}$, which we can express as

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

where

\begin{equation}
    {\bf d}({\bf m})={\bf S}\big[\mathcal{L}({\bf m}) \big]
\end{equation}

is the modeled data calculated by using the full-wave propagation operator, $\mathcal{L}$, collected at the receivers' position by the sampling operator, ${\bf S}$.

The solution of the FWI optimization problem in equation \eqref{eq:fwi1} using gradient descent optimization schemes or Hessian based methods \citep{Nocedal2006} requires the computation of the gradient of the misfit function, which is given by

\begin{equation}
\nabla_{\bf m} \Phi({\bf m})=\big [\nabla_{\bf m}\mathcal{L}({\bf m}) \big]^{\text{T}}{\bf S}^{\text{T}} \big[ {\bf S}\mathcal{L}({\bf m})-{\bf d}_{\text{obs}} \big]={\bf L}({\bf m})^{\text{T}}{\bf S}^{\text{T}}\big[ {\bf S}\mathcal{L}({\bf m})-{\bf d}_{\text{obs}} \big], \label{eq:fwi2}
\end{equation}

where $\nabla_{\bf m}$ constitutes the gradient with respect to the subsurface model parameters, ${\bf S}^{\text{T}}$ is the adjoint operator of sampling (hence it injects data at the receivers' positions), and ${\bf S}\mathcal{L}({\bf m})-{\bf d}_{\text{obs}}$ represents the data residuals that are projected into the model space by means of the RTM operator, ${\bf L}({\bf m})^{\text{T}}=\big[\nabla_{\bf m} \mathcal{L}({\bf m})\big]^{\text{T}}$.

The FWI solution can be obtained by local linearization where the current subsurface model, ${\bf m}_{\text{i}}$, is updated by a direction vector ${\bf \delta m}$ such that $\|{\bf \delta m} \| << \|{\bf m}_{\text{i}} \|$, thereby exhibiting a local quadratic behavior. Now we can expand the gradient in Taylor's series around ${\bf m}_{\text{i}}$ and drop second- and higher-order terms, obtaining

\begin{equation}
\nabla_{\bf m} \Phi({\bf m}_{\text{i}}+{\bf \delta m}) \approx \nabla_{\bf m} \Phi({\bf m}_{\text{i}})+\nabla_{\bf m} \nabla_{\bf m} \Phi({\bf m}_{\text{i}}) {\bf \delta m} = \nabla_{\bf m} \Phi({\bf m}_{\text{i}})+{\bf H}_{\text{fwi}}({\bf m}_{\text{i}}) {\bf \delta m}, \label{eq:fwi3}
\end{equation}

where ${\bf H}_{\text{fwi}}({\bf m}_{\text{i}})$ is the FWI Hessian, not be confused with the FWI {\it Gauss-Newton} Hessian, ${\bf H}={\bf L}^{\text{T}}{\bf L}$, introduced in a previous section. As mentioned before, solving equation~\eqref{eq:fwi3} represents the optimization of a quadratic problem in ${\bf \delta m}$. Therefore, the gradient evaluated at the updated model, ${\bf m}_{{\text{i}}+1}={\bf m}_{\text{i}}+{\bf \delta m}$, becomes zero, and equation \eqref{eq:fwi3} reduces to

\begin{equation}
{\bf H}_{\text{fwi}}({\bf m}_{\text{i}}) {\bf \delta m}=-\nabla_{\bf m} \Phi({\bf m}_{\text{i}}). \label{eq:fwi4}
\end{equation}

Equation \eqref{eq:fwi4} represents the well-known Newton's equation \citep{Nocedal2006}.

Now I will incorporate equation \eqref{eq:not3} into Newton's equation by defining the current model and the model step as follows:

\begin{equation}
    \begin{aligned}
    {\bf m}_{\text{i}}&={\bf b}_0, \\
    {\bf \delta m}&={\bf \Delta b}+{\bf r}. \label{eq:newton1}
    \end{aligned}
\end{equation}

On the other hand, according to \citet{Biondi.sep.160.biondo1} we can split the FWI Hessian into the FWI Gauss-Newton Hessian and the adjoint of the WEMVA operator:

\begin{equation}
    {\bf H}_{\text{fwi}}={\bf H}+{\bf W}^{\text{T}}. \label{eq:newton2}
\end{equation}

Substituting equations \eqref{eq:newton1} and \eqref{eq:newton2} into the Newton's equation yields

\begin{equation}
{\bf H}_{\text{fwi}}({\bf m}_{\text{i}}) {\bf \delta m}=\big[{\bf H}({\bf b}_0)+{\bf W}({\bf b}_0)\big] ({\bf \Delta b}+{\bf r}) . \label{eq:newton3}
\end{equation}

As far as I am assuming the remaining inaccuracies in the background model to be small, it is valid to employ maximization of the image energy. Therefore, I dropped the transpose symbol in ${\bf W}$ because the WEMVA operator is self-adjoint in the zero-subsurface offset domain.

Re-arranging terms in the right-hand side of equation \eqref{eq:newton3} we obtain

\begin{equation}
\big[{\bf H}({\bf b}_0){\bf r}+{\bf W}({\bf b}_0){\bf \Delta b}\big]+
\big[{\bf H}({\bf b}_0){\bf \Delta b}+{\bf W}({\bf b}_0){\bf r})\big]. \label{eq:newton4}
\end{equation}

The first part inside brackets in expression \eqref{eq:newton4} contains the image term ${\bf H}({\bf b}_0){\bf r}$ employed in model-space LWI, and the term ${\bf W}({\bf b}_0) {\bf \Delta b}$ that represents the perturbation in the image regarding the correction of the background model. I keep both terms because they have a predictable impact on primary events related to the subsurface reflectivity. The second part in brackets very likely contributes to second-order scattering events $\big[{\bf W}({\bf b}_0) {\bf r}\big]$ and refinement of the low-wavenumber component $\big[{\bf H \Delta b}\big]$. I drop the first term because it is not related to the subsurface reflectivity's primary events, and I drop the second one because the method does not aim to retrieve a high-wavenumber background model.

Regarding the FWI gradient in the Newton's equation (right-hand side of equation \eqref{eq:fwi4}), it becomes $\nabla \Phi({\bf b}_0)$. Substituting this gradient in equation \eqref{eq:fwi2} yields

\begin{equation}
\nabla \Phi({\bf b}_0)=\big [\nabla\mathcal{L}({\bf b}_0) \big]^{\text{T}} 
{\bf S}^{\text{T}} \big[ {\bf S} \mathcal{L}({\bf b}_0)-{\bf d}_{\text{obs}} \big]. \label{eq:grad-b0}
\end{equation}

Nonlinear wave propagation at the background model (free of reflections) will model only direct waves and diving waves. Therefore, the term ${\bf S}^{\text{T}} \big[ {\bf S} \mathcal{L}({\bf b}_0)-{\bf d}_{\text{obs}} \big]$ simply represents the negative of the perturbation in the data, ${\bf \Delta d}_{\text{obs}}$, \citep{Barnier.sep.155.gbarnier1}, i.e., the data after removing the direct and diving waves. On the other hand, the adjoint of the gradient term represents Born adjoint modeling or RTM evaluated at ${\bf b}_0$. Thus, equation \eqref{eq:grad-b0} becomes the negative of the RTM image,

\begin{equation}
\nabla \Phi({\bf b}_0)=-{\bf L}({\bf b}_0)^{\text{T}} {\bf \Delta d}_{\text{obs}}=-{\bf I}({\bf b}_0). \label{eq:grad-rtm}
\end{equation}

Substituting the first term in braces in expression \eqref{eq:newton4} and equation \eqref{eq:grad-rtm} in the Newton's equation we obtain

\begin{equation}
{\bf H}({\bf b}_0){\bf r}+{\bf W}({\bf b}_0){\bf \Delta b}={\bf I}({\bf b}_0), \label{eq:orig_lwivu1}
\end{equation}

which I named Linearized Waveform Inversion with Velocity Updating (LWIVU), and is equivalent to minimizing the objective function

\begin{equation}
\Phi({\bf r},{\bf \Delta b})=\frac{1}{2} \left\Vert{\bf H}({\bf b}_0){\bf r}-{\bf I}({\bf b}_0)+{\bf W}({\bf b}_0){\bf \Delta b} \right\Vert_2^2. \label{eq:orig_lwivu2}
\end{equation}

Note that this objective function represents the first functional in equation \eqref{eq:lwivu-zeroh}. The only difference is in the positive sign for ${\bf W}({\bf b}_0){\bf \Delta b}$, which is not a problem because the perturbation in the background can take positive or negative values.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Nonlinear scheme}
Given the shortcoming of the objective function \eqref{eq:lwivu-zeroh}, \citet{Cabrales.sep.176.alejandro1} proposed to keep the nonlinear objective function for the migration image, still maintaining the assumption $||\mathbf{\Delta b}||_2 << ||\mathbf{b}_0||_2$, so that the approximation $\mathbf{H}(\mathbf{b}_0 + \mathbf{\Delta b}) \approx \mathbf{H}(\mathbf{b}_0)$ is accurate. Therefore, the JIRB objective function became

\begin{equation}
    \Phi(\mathbf{r}, \mathbf{\Delta b}) = \left \Vert \mathbf{H}(\mathbf{b}_0)\mathbf{r}
    - \mathbf{W}(\mathbf{b}_0)\mathbf{\Delta b} - \mathbf{I}(\mathbf{b}_0) \right \Vert _2^2
    - \lambda \left \Vert \mathbf{I}(\mathbf{b}_0 + \mathbf{\Delta b}) \right \Vert_2^2.
\label{eq:new_nlwivu}
\end{equation}

This change makes the first term a linear least-squares problem over both $\mathbf{\Delta b}$ and $\mathbf{r}$, while the second term is nonlinear and non-quadratic.

Although the objective function \eqref{eq:new_nlwivu} behaves as desired, it is more straightforward and sensible to implement the non-linearity in both objective functions:

\begin{equation}
    \Phi(\mathbf{r}, \mathbf{b}) = \left \Vert \mathbf{H}(\mathbf{b}_0)\mathbf{r}
    - \mathbf{I}(\mathbf{b}) \right \Vert _2^2
    - \lambda \left \Vert \mathbf{I}(\mathbf{b}) \right \Vert_2^2.
\label{eq:jirb_final}
\end{equation}

Objective function \eqref{eq:jirb_final} can be solved (for local minima) using gradient-based methods such as steepest descent, nonlinear conjugate gradients, or quasi-Newton \citep{Nocedal2006}. Notice that the Gauss-Newton Hessian is precomputed with the wrong background model, and does not require to be updated as it is the case for the migration image. Solving objective function \eqref{eq:jirb_final} is the core of the JIRB method. Chapter 4 shows a variation where I merely add preconditioning and weighting factors.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Extended domain}
The joint inversion of reflectivity and perturbation in the background can also be cast in the extended domain as

\begin{equation}
\Phi({\bf \tilde{r}},{\bf \Delta b})=\frac{1}{2} \left\Vert{\bf \tilde{H}}({\bf b}_0+{\bf \Delta b}){\bf \tilde{r}}-{\bf \tilde{I}}({\bf b}_0+{\bf \Delta b}) \right\Vert_2^2+\frac{\lambda}{2} \left\Vert {\bf D}{\bf \tilde{I}}({\bf b}_0+{\bf \Delta b})\right\Vert_2^2, \label{eq:jinv-nl-ext}
\end{equation}

where the tilde represents extension in subsurface offset and ${\bf D}$ represents the differential semblance optimization operator (DSO) \citep{Symes1991}. In this case, the second functional ought to be minimized to penalize non-zero subsurface offset energy originated from incorrect background model.

At first glance it appears that in the expanded case the second functional can be linearized by substituting ${\bf D}{\bf \tilde{I}}({\bf b}_0+{\bf \Delta b}) \approx {\bf D}\big[{\bf \tilde{I}}({\bf b}_0)+{\bf \tilde{W}}({\bf b}_0){\bf \Delta b}\big]$ in the norm, because both the nonlinear and the linearized functionals are convex, in contrast to the non-extended case where we maximize the image power, analyzed before. However, this approximation still lacks one second-order term that can lead to incorrect results (see Appendix \ref{appendix:lin}).

We can use a similar objective function as that of the zero subsurface offset domain:

\begin{equation}
\Phi({\bf \tilde{r}},{\bf b})=\frac{1}{2} \left\Vert{\bf \tilde{H}}({\bf b}_0){\bf \tilde{r}}-{\bf \tilde{I}}({\bf b}) \right\Vert_2^2+\frac{\lambda}{2} \left\Vert {\bf D}{\bf \tilde{I}}({\bf b}_0+{\bf \Delta b})\right\Vert_2^2. \label{eq:jinv-nl-ext-dso}
\end{equation}

The obvious caveat of equation \eqref{eq:jinv-nl-ext-dso} is that inversion in the extended domain is more costly than the inversion in the zero offset domain because the scattering operator and the imaging condition operator dominate the computations. This problem is particularly severe in 3D imaging. Some people have proposed strategies to alleviate the problem, such as the so-called ``Nyquist approach"\citep{Sun2013}, where we reduce to the minimum the number of times that the imaging condition is applied during RTM, depending on the frequency range of the data. However, in order to achieve enough extension, the computational burden is still considerable. The DSO operator often requires special preconditioning to avoid high-wavenumber artifacts in the gradient \citep[e.g.][]{Shen2008,Tang.sepphd.144,Weibull2013}. For such reasons, in this dissertation I implement the JIRB methodology exclusively in the zero subsurface offset domain.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section{Alternative implementation}

%There is still another alternative to invert for the reflectivity and the background components, which is worth to analyze. We can recast the Born modeling operator in data-space LWI as

%\begin{equation}
%\Phi({\bf r},{\bf \Delta b})=\frac{1}{2}\left\Vert {\bf L}({\bf b}_0+{\bf \Delta b}){\bf r}-{\bf d}_{\text{obs}} \right\Vert_2^2. \label{eq:born-jirb}
%\end{equation}

%We can linearize the Born operator applied to the reflectivity as follows

%\begin{equation*}
%{\bf L}({\bf b}_0+{\bf \Delta b}){\bf r}={\bf L}({\bf b}_0){\bf r}+\bigg[\frac{\partial {\bf L}({\bf b}_0)} {\partial \bf {b}} {\bf r}\bigg] {\bf \Delta b}+O(\|{\bf \Delta b}\|^2).
%\end{equation*}

%Dropping the higher-order terms in ${\bf \Delta b}$, and defining

%\begin{equation*}
%{\bf T}({\bf b}_0,{\bf r})=\frac{\partial {\bf L}({\bf b}_0)} {\partial \bf {b}} {\bf r},
%\end{equation*}

%we obtain

%\begin{equation}
%{\bf L}({\bf b}_0+{\bf \Delta b}){\bf r} \approx {\bf L}({\bf b}_0){\bf r}+{\bf T}({\bf b}_0,{\bf r}) {\bf \Delta b}, \label{eq:born-taylor}
%\end{equation}

%were ${\bf T}$ represents the tomographic operator, which relates perturbations in the background to perturbations in the data \citep{Barnier.sep.155.gbarnier1}. The counterpart is the WEMVA operator, which relates perturbations in the background to perturbations in the image.

%Substituting equation \eqref{eq:born-taylor} into equation \eqref{eq:born-jirb} yields

%\begin{equation}
%\Phi({\bf r},{\bf \Delta b})=\frac{1}{2}\left\Vert {\bf L}({\bf b}_0){\bf r}+{\bf T}({\bf b}_0,{\bf r}) {\bf \Delta b}-{\bf d}_{\text{obs}} \right\Vert_2^2. \label{eq:born-linear}
%\end{equation}

%Equation \eqref{eq:born-linear} has the inconvenience that the tomographic term contains both the reflectivity and the perturbation in the background, i.e., the inversion parameters. Moreover, the reflectivity is part of the tomographic operator, which makes it difficult to decouple it from the derivative of the Born modeling operator. Such a derivative is not a matrix, but a higher-order tensor (it becomes a matrix after pre-multiplying the reflectivity vector). For these reasons, in this dissertation I do not explore this possible variant of JIRB.