# Robust Full Waveform Inversion with deep Hessian deblurring

Mustafa Alfarhan
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Matteo Ravasi
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Fuqiang Chen
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Tariq Alkhalifah
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>

###### Abstract

Full Waveform Inversion (FWI) is a technique widely used in geophysics to
obtain high-resolution subsurface velocity models from waveform seismic data.
Due to its large computation cost, most flavors of FWI rely only on the
computation of the gradient of the loss function to estimate the update
direction, therefore ignoring the contribution of the Hessian. Depending on
the level of computational resources one can afford, an approximate of the
inverse of the Hessian can be calculated and used to speed up the convergence
of FWI towards the global (or a plausible local) minimum. In this work, we
propose to use an approximate Hessian computed from a linearization of the
wave-equation as commonly done in Least-Squares Migration (LSM). More
precisely, we rely on the link between a migrated image and a doubly migrated
image (i.e., an image obtained by demigration-migration of the migrated image)
to estimate the inverse of the Hessian. However, instead of using non-
stationary compact filters to link the two images and approximate the Hessian,
we propose to use a deep neural network to directly learn the mapping between
the FWI gradient (output) and its Hessian (blurred) counterpart (input). By
doing so, the network learns to act as an approximate inverse Hessian: as
such, when the trained network is applied to the FWI gradient, an enhanced
update direction is obtained, which is shown to be beneficial for the
convergence of FWI. The weights of the trained (deblurring) network are then
transferred to the next FWI iteration to expedite convergence. We demonstrate
the effectiveness of the proposed approach on two synthetic datasets and a
field dataset.

_Keywords_ Full waveform inversion $\cdot$ Imaging $\cdot$ Born modeling
$\cdot$ Deep learning $\cdot$ Neural networks

##  1 Introduction

Seismic data, commonly acquired at the surface of the Earth, contain useful
information about the Earth’s structure. Seismic waves travel in the
subsurface, interact with it, and some of them eventually return to the
surface to be recorded by an array of receivers. However, inferring the
subsurface structure and its medium properties from seismic data is a
challenging task. This difficulty arises mostly from the fact that the
subsurface is accessed from only one side, and therefore some areas are not
covered well by the recorded waves, or in other words, they are poorly
illuminated. Moreover, multiply scattered waves (i.e., free-surface and
internal multiples) interfere with primary arrivals, which present a direct
linear link with the subsurface properties of interest. As a consequence,
geophysical research over the years has been focused on developing methods
that can separate different parts of the seismic data and optimally utilize
them to image the subsurface. Full Waveform Inversion (FWI) was introduced by
Lailly and Bednar (1983); Tarantola (1984) as an alternative approach that
aims to utilize the entire seismic waveform data at once to construct a high-
resolution subsurface model. More precisely, FWI casts the problem of
transforming seismic records into a subsurface model as an optimization
problem, where a measure of the difference between modeled and observed data
is minimized. Given an initial subsurface model, obtained, for example,
through traveltime tomography Bording et al. (1987) or migration velocity
analysis Symes (2008), synthetic seismic data are calculated with a wave-
equation based forward modeling engine. The data misfit between modeled and
observed data is then utilized to update the model parameters and the process
is repeated until convergence.

As a consequence of its pursuit in utilizing the entire seismic data content
at once, FWI is known to be a highly non-linear inverse problem. Whether or
not FWI converges to the right model depends on the initial model: more
specifically, if the travel-time difference between the modeled and observed
data is more than half a cycle, we will converge to a local minimum model, a
problem known as cycle-skipping. To mitigate cycle-skipping, FWI is usually
performed in a multi-scale fashion Bunks et al. (1995); Miller and Willsky
(1996); Fichtner et al. (2013); Alkhalifah (2016); Zhang et al. (2019) where
the inversion is applied to the low-frequency components of the data first.
Only after the kinematic part of the subsurface model is correctly inverted
for to ensure that we are within the half-cycle difference for the higher
frequencies that the high-frequency components of the data are incorporated to
add high-resolution features to the inverted model. Although there exist both
local Lailly and Bednar (1983); Tarantola (1984) and global Jin and Madariaga
(1994) approaches to solving FWI, the former is usually preferred: each model
update involves calculating the first and (possibly) second derivatives of the
data-misfit with respect to the current model, which are called the gradient
and the Hessian, respectively. Whilst the gradient provides the direction of
the model update, the inverse of the Hessian alters the direction slightly in
compliance with the curvature information and introduces a parameter-dependent
scaling to the gradient. From a physical standpoint, the inverse of the
Hessian accounts for losses in amplitude due to geometrical spreading and
transmission losses occurring in the wave equation, as well as it provides a
trade-off between different model parameters in multi-parameter FWI Alkhalifah
(2016). However, calculating the inverse of the Hessian for large subsurface
models is computationally expensive.

Traditionally, gradient-based methods have been popularized in FWI
applications as they only rely on the gradient for the model update
calculation. Among gradient-based methods, the most popular ones are the
steepest-descent method, whose update is given by a scaled version of the
negative of the gradient, and the non-linear conjugate-gradient method, whose
update is a combination of the gradient at the current iteration and the
previous update. On the other hand, second-order methods such as the Newton
method rely on the Hessian in addition to the gradient to define the model
update. A special case, when second-order derivative terms in the Hessian
(which account for second-order scattering) are neglected, is represented by
the Gauss-Newton (GN) method Pratt et al. (1998). Different strategies have
emerged over the years to approximate the GN Hessian. For example, the
gradient can be preconditioned by the diagonal elements of the Hessian or the
pseudo-Hessian as proposed by Shin et al. (2001). On the other hand, the model
and gradient from previous iterations can be used to iteratively update the
inverse of the Hessian, as done in the family of quasi-Newton methods such as
the limited-memory BFGS algorithm Liu and Nocedal (1989). Another alternative
to estimating the Hessian is the family of truncated-Newton methods, which do
not need the Hessian explicitly but rather compute the Hessian-vector product
Métivier et al. (2014); AlTheyab et al. (2013). An equivalently accurate
approach to approximate the GN Hessian, and is the inspiration to the approach
developed here, is a least-squares optimization of the gradient Wu and
Alkhalifah (2015). In general, the inclusion of the approximate inverse of the
Hessian provides faster convergence than gradient-based methods.

Seismic imaging is another area where approximating the Hessian is of great
importance as the resolution of the resulting seismic images can be enhanced
and uneven illumination resulting from imperfect acquisition geometries can be
mitigated. More specifically, Least-Squares Migration (LSM) is a technique
that aims to mitigate the degrading effects of the Hessian. In the image
domain, the Hessian effect can be considered a form of blurring operator.
Therefore, LSM can be formulated as an image deblurring problem where the
inverse of the Hessian is the deblurring operator. Using a migrated and doubly
migrated image, the inverse of the Hessian can be approximated with non-
stationary matching filters as done in Guitton (2004). Similarly, Aoki and
Schuster (2009) use deblurring filters to accelerate LSM in a regularized or
preconditioning fashion. Alternatively, as initially proposed by Valenciano
(2008), one can rely on point spread functions (PSFs) to construct the
deblurring filters. In recent years, various deep learning methods have been
proposed for image deblurring (which are equivalently applicable to the
problem of approximating the action of the inverse of the Hessian matrix). Liu
et al. (2022) suggested to employ a convolutional neural network to learn the
deconvolution operator of the PSFs in least-squares reverse time migration
(LSRTM). Kumar et al. (2022) uses a neural network to learn the inverse of the
Hessian using the original migration image and a remigrated image. Finally,
Peng et al. (2023) suggest to use a loop-unrolled gradient descent method to
include the contribution of the known PSF in the learning process. In all
studies, it was concluded that using a network to approximate the action of
the inverse of the Hessian is faster than other state-of-the-art methods and
can more stably deconvolve seismic images.

In this work, we propose to embed a similar approach to that of Kumar et al.
(2022) within FWI. More specifically, Born modeling and its adjoint are used
to create our Hessian blurred gradient. In this case, the Hessian blurred
gradient is approximately related to the FWI gradient through the inverse of
the Hessian, whose mapping is learned using a deep neural network. After
training, the network is applied to the gradient to obtain a Hessian corrected
update that we demonstrate can speed up FWI convergence. We apply our approach
to two synthetic data to showcase its benefits, and to the Volve field data to
validate its applicability under field conditions.

##  2 Theoretical Framework

###  2.1 Full waveform inversion

The constant-density acoustic wave equation is given by

$\frac{1}{m^{2}(\mathbf{x})}\frac{\partial^{2}u(\mathbf{x},t)}{\partial
t^{2}}-\nabla^{2}u(\mathbf{x},t)=f(\mathbf{x},t)$ (1)

where $m(\mathbf{x})$ is the medium velocity, $u(\mathbf{x},t)$ is the
wavefield, $f(\mathbf{x},t)$ is the source function, all defined in space
$\mathbf{x}$ and time $t$, and $\nabla^{2}$ is the Laplacian operator.
Starting from an initial velocity model, the wave equation is solved using a
forward modeling scheme such as the finite difference method to obtain the
simulated data, $\mathbf{d}_{s}$. Then, the simulated data are compared to the
observed data, $\mathbf{d}_{o}$, using a misfit function such as the $L_{2}$
norm:

$J(\mathbf{m})=\frac{1}{2}\Delta\mathbf{d}^{T}\Delta\mathbf{d}$ (2)

where $\Delta\mathbf{d}=\mathbf{d}_{o}-\mathbf{d}_{s}$. Note that in the above
equation we have expressed both datasets and the model in vector form. From
now on, the use of regular or bold symbols will discriminate whether we are
considering a scalar or a vector for our equations.

Using the generalized Taylor series expansion around a background model
$\mathbf{m}_{o}$, the misfit function can be approximated by

$J(\mathbf{m})\approx J(\mathbf{m}_{o})+\left(\frac{\partial
J(\mathbf{m}_{o})}{\partial\mathbf{m}}\right)^{T}(\mathbf{m}-\mathbf{m}_{o})+\frac{1}{2}(\mathbf{m}-\mathbf{m}_{o})^{T}\frac{\partial^{2}J(\mathbf{m}_{o})}{\partial\mathbf{m}^{2}}(\mathbf{m}-\mathbf{m}_{o})$
(3)

Taking the derivative of equation (3) with respect to
$\boldsymbol{\delta}\mathbf{m}=\mathbf{m}-\mathbf{m}_{o}$ and setting it to
zero yields

$\boldsymbol{\delta}\mathbf{m}=-\left(\frac{\partial^{2}J(\mathbf{m}_{o})}{\partial\mathbf{m}^{2}}\right)^{-1}\frac{\partial
J(\mathbf{m}_{o})}{\partial\mathbf{m}}$ (4)

where $\boldsymbol{\delta}\mathbf{m}$ is the model update, $\partial
J(\mathbf{m}_{o})/\partial\mathbf{m}$ is the gradient, and
$\left(\partial^{2}J(\mathbf{m}_{o})/\partial\mathbf{m}^{2}\right)^{-1}$ is
the inverse of the Hessian. FWI is an optimization process where the current
model $\mathbf{m}$ is iteratively updated by generating the data
$\mathbf{d}_{s}$ and comparing them against the observed data to compute the
model update $\boldsymbol{\delta}\mathbf{m}$ accordingly until the data
residual is sufficiently small.

###  2.2 Born modeling

In Born modeling, the model $m(\mathbf{x})$ is commonly split into two parts:
background model and perturbation. The background model is often smooth, while
the perturbation may lead to scattering in the wavefield. Therefore, the
velocity $m(\mathbf{x})$ can be written as a combination of a background
velocity $m_{o}(\mathbf{x})$ and a perturbation $\delta m(\mathbf{x})$

$m(\mathbf{x})=m_{o}(\mathbf{x})+\delta m((\mathbf{x})$ (5)

In addition, we can approximate $1/{m^{2}(\mathbf{x})}$ using Taylor series as

$\frac{1}{m^{2}(\mathbf{x})}\approx\frac{1}{m_{o}^{2}(\mathbf{x})}-2\frac{\delta
m(\mathbf{x})}{m_{o}^{3}(\mathbf{x})}$ (6)

Similarly, the seismic wavefield $u(\mathbf{x},t)$ can be split into a
background wavefield $u_{o}(\mathbf{x},t)$, satisfying the wave equation for
$m_{o}(\mathbf{x})$ and a scattering wavefiled $\delta u(\mathbf{x},t)$ such
that

$u(\mathbf{x},t)=u_{o}(\mathbf{x},t)+\delta u(\mathbf{x},t)$ (7)

Substituting equations (6) and (7) into equation (1) yields

$\frac{1}{m_{o}^{2}(\mathbf{x})}\frac{\partial^{2}u_{o}(\mathbf{x},t)}{\partial
t^{2}}-\nabla
u_{o}(\mathbf{x},t)+\frac{1}{m_{o}^{2}(\mathbf{x})}\frac{\partial^{2}\delta
u(\mathbf{x},t)}{\partial t^{2}}\\\ -\nabla\delta
u(\mathbf{x},t)-2\frac{\delta
m(\mathbf{x})}{m_{o}^{3}(\mathbf{x})}\frac{\partial^{2}(u_{o}(\mathbf{x},t)+\delta
u(\mathbf{x},t))}{\partial t^{2}}=f(\mathbf{x},t)$ (8)

which can be split into the background and the scattering parts, respectively:

$\frac{1}{m_{o}^{2}(\mathbf{x})}\frac{\partial^{2}u_{o}(\mathbf{x},t)}{\partial
t^{2}}-\nabla u_{o}(\mathbf{x},t)=f(\mathbf{x},t)$ (9)
$\frac{1}{m_{o}^{2}(\mathbf{x})}\frac{\partial^{2}\delta
u(\mathbf{x},t)}{\partial t^{2}}-\nabla\delta u(\mathbf{x},t)\\\
=2\frac{\delta
m(\mathbf{x})}{m_{o}^{3}(\mathbf{x})}\frac{\partial^{2}(u_{o}(\mathbf{x},t)+\delta
u(\mathbf{x},t))}{\partial t^{2}}$ (10)

Equation (10) is implicit since both the left- and right-hand sides depends on
the scattering wavefield. To make it explicit, we use the Born approximation,
which assumes that the scattered wavefield is small compared to the background
wavefield,

$\frac{1}{m_{o}^{2}(\mathbf{x})}\frac{\partial^{2}\delta
u(\mathbf{x},t)}{\partial t^{2}}-\nabla\delta u(\mathbf{x},t)=2\frac{\delta
m(\mathbf{x})}{m_{o}^{3}(\mathbf{x})}\frac{\partial^{2}u_{o}(\mathbf{x},t)}{\partial
t^{2}}$ (11)

Consequently, Born modeling can be performed by first modeling the background
wavefield $u_{o}(\mathbf{x},t)$ in the background model $m_{o}(\mathbf{x})$
with source function $f(\mathbf{x},t)$ using equation (9). This wavefield is
then used to ignite a secondary source in equation (10); in other words, a
second step of modeling is therefore performed in the background model to
compute the scattered wavefield using this secondary source. Note that, in
practice, equations (9) and (10) are jointly solved so that there is no need
to store the entire background wavefield from the first modeling to be
injected in the second one.

In compact notation, we can define the linear Born modeling as

$\boldsymbol{\delta}\mathbf{u}=\mathbf{L}\boldsymbol{\delta}\mathbf{m}$ (12)

with
$\boldsymbol{\delta}\mathbf{m}=\boldsymbol{\delta}\mathbf{m}/\mathbf{m}_{o}$.
Finally, an additional sampling operator $\mathbf{P}$ is used to extract the
scattered wavefield at the location of the receivers and produce the surface
seismic data:

$\mathbf{d}=\mathbf{PL}\boldsymbol{\delta}\mathbf{m}=\tilde{\mathbf{L}}\boldsymbol{\delta}\mathbf{m}$
(13)

where, for simplicity, we define $\tilde{\mathbf{L}}$ from here onwards as the
Born modeling operator.

###  2.3 Approximating the inverse of the Hessian for FWI with deep learning

The full model update for FWI as shown in equation (4) requires the
calculation of the gradient and the inverse of the Hessian. However, the
calculation of the inverse of the Hessian is highly expensive: this has led to
the development of methods to estimate an approximate version of such an
inverse with affordable computational cost. In this work, we propose a method
to approximate the Hessian that relies on deep learning. The full Hessian is
composed of two terms: a single-scattering and a double-scattering terms. In
the Gauss-Newton method, the second term is neglected and only the first term
is used. Multiplying the gradient by the inverse of the Hessian can be viewed
as a deblurring operation that corrects the gradient for the effects of the
Earth such as geometric spreading. We present an approach for approximating
the inverse of the Hessian for FWI in the following sub-section.

####  2.3.1 The Born modeling and its adjoint approach

One way to approximate the Hessian is via a migration/demigration approach,
or, after linearizing FWI, Born modeling and its adjoint. The gradient of FWI
can be summarized as follows

$\mathbf{g}=\tilde{\mathbf{L}}^{T}\Delta\mathbf{d}$ (14)

where $\tilde{\mathbf{L}}^{T}$ is the adjoint of the Born modeling operator,
and $\Delta\mathbf{d}$ is the data residual. Moreover, considering the true
Earth perturbation is $\boldsymbol{\delta}\mathbf{m}$,
$\Delta\mathbf{d}=\tilde{\mathbf{L}}\boldsymbol{\delta}\mathbf{m}$, equation
(14) can be equivalently written as

$\mathbf{g}=\tilde{\mathbf{L}}^{T}\tilde{\mathbf{L}}\boldsymbol{\delta}\mathbf{m}$
(15)

which reveals that the gradient $\mathbf{g}$ and the true perturbation
$\boldsymbol{\delta}\mathbf{m}$ are linked via a pair of Born modeling and its
adjoint (demigration/migration) operations. A blurred gradient (referred in
this work as doubly migrated image) can be obtained from the gradient by
applying once again Born modeling and its adjoint

$\boldsymbol{\delta}\mathbf{m}_{1}=\tilde{\mathbf{L}}^{T}\tilde{\mathbf{L}}\mathbf{g}$
(16)

Given availability of both $\mathbf{g}$ and
$\boldsymbol{\delta}\mathbf{m}_{1}$, one can approximate
$(\tilde{\mathbf{L}}^{T}\tilde{\mathbf{L}})^{-1}$ (the inverse of the Hessian)
and apply it to $\mathbf{m}$ to obtain an estimate of
$\boldsymbol{\delta}\mathbf{m}$, which is the objective of LSM. In our
approach, $\mathbf{m}$ is the FWI gradient, $\mathbf{g}$. In other words, by
estimating and applying the inverse of the Hessian to the FWI gradient, we aim
to retrieve an approximation of the scattering potential as used in Born
modeling.

In this work, we suggest to train a neural network to learn the mapping from
$\boldsymbol{\delta}\mathbf{m}_{1}$ to $\mathbf{g}$, which ultimately
represents the action of the inverse of the Hessian on a vector. Also, since
the FWI update is related to the gradient through the Hessian (see equation
(4)), the trained network can be applied to the gradient to obtain an improved
FWI update; this can ultimately lead to faster convergence of the FWI process.
The same process is repeated at every iteration, except for the fact that the
neural network is not trained from scratch: instead, we leverage transfer
learning and start from the weights learned in the previous iteration. The
proposed method is depicted in Figure 1.

Figure 1: Proposed method to approximate the inverse of the Hessian using deep
learning.

###  2.4 Neural-network based inverse Hessian approximation

In this work, we suggest to train a neural network to approximate the inverse
of the Hessian in FWI. Given that learning the action of the inverse of the
Hessian can be framed as a domain translation process (from the doubly
migrated to the migrated image), our network presents a UNet-like architecture
Ronneberger et al. (2015). More specifically, the network architecture is
composed of four blocks in both the encoder and decoder parts. The input is
passed to an input block which consists of a double convolution layer which
itself contains two sequences of operations. Each sequence involves a 2D
convolution with a kernel size of $3\times 3$, stride of 1, and padding of 1,
followed by a Leaky ReLU activation function with a negative slope of 0.2. The
first convolutional layer transforms the input channel to 128 channels, and
the second maintains these 128 channels. There are four downsampling blocks in
the encoder, each consisting of an average pooling operation with a $3\times
3$ kernel and a stride of 2 followed by a double convolution block (similar to
the input block). These blocks progressively increase the number of channels
while reducing the spatial dimensions of the feature maps. The channel sizes
double at each block, going from 128 to 256, then to 512, 1024, and finally
2048. Corresponding to the downsampling blocks, there are four upsampling
blocks in the decoder. Each block starts with a transposed 2D convolution
which halves the number of channels, and is followed by a double convolution
block (similar structure as the encoder). These blocks incrementally reduce
the number of channels while increasing the spatial dimensions of the feature
maps, effectively recovering the spatial resolution that was lost during
downsampling. The final part of the network is an output convolution layer. It
consists of a convolution layer with a kernel size of $1\times 1$, used to
transform the 128 channels from the last upsampling block to the desired
number of output channels (in this case, 1 channel). The network architecture
is summarized in Table 1 in Appendix A. The Adam optimizer is used for
optimizing the network’s parameters with a cosine annealing scheduler for the
learning rate. The MSE loss function is used to measure the misfit between the
network’s output and the FWI gradient. The training process can be carried out
using patches or the full image at once. Whilst both approaches have similar
performance, we chose the latter because it makes training faster.

##  3 Numerical Examples

In this section, we test the approach on two synthetic examples and one field
data used to validate the proposed methodology.

###  3.1 Marmousi model

Due to its intricate geological composition, the synthetic Marmousi model
Brougois et al. (1990) serves as a benchmark for assessing algorithms in
seismic imaging and inversion. The model used in this example spans $601$ grid
points horizontally and $221$ grid points vertically, with a grid size of $15\
m$ in both directions. Figure 2 illustrates both the Marmousi model and its
smoothed counterpart, utilized as the starting model for FWI. The observed
data is modeled using a variable-velocity, constant-density acoustic finite-
difference modeling engine, using 30 sources and 300 receivers equally spaced
at zero depth. A Ricker wavelet of 5 Hz peak frequency is used for the source
signature.

Figure 2: a) the true Marmousi model and b) its smoothed version. Figure 3:
The components of our approach: a) the FWI gradient, b) $\delta m_{1}$, c) the
predicted gradient by the network, and d) $\delta m$.

To begin with, we consider the first iteration of the FWI process to detail
the proposed approach. Starting from the FWI gradient (Figure 3a), we apply
Born modeling and its adjoint to obtain the doubly-migrated image $\delta
m_{1}$ in Figure 3b. Comparing the gradient to the doubly migrated image, we
can see that the latter is of lower resolution and presents a stronger
amplitude decay over depth; this is expected as the effect of geometrical
spreading is accumulated multiple times as part of the Born modeling and its
adjoint. Figure 3c shows the network output after training. This result
highlights that the model has learned to map $\delta m_{1}$ to the gradient.
Also, we show the normalized training loss for the first three iterations of
FWI in Figure 4 which highlights that the network is learning to map $\delta
m_{1}$ to the gradient. Note that in the first FWI iteration the network is
trained for 1000 iterations, but we only show the first 300 iterations just to
make the plot more presentable. In addition, we should point out that the
fitting of the gradient is not perfect (e.g., there is an amplitude loss)
which could be a result of not having enough information to go from $\delta
m_{1}$ to the gradient. Comparing the network’s gradient to the FWI gradient
we can say that the network has learned an approximation of the inverse
Hessian. Thus, if we apply the network to the FWI gradient, we get a Hessian
corrected update as shown in Figure 3d. For certain, we do not expect the
network to learn the exact inverse Hessian due to the linear link between the
gradient and $\delta m_{1}$ whilst the FWI inverse Hessian is non-linear. Not
surprisingly, the update presents more balanced contributions throughout the
entire depth axis and at the edges of the horizontal axis compared to the
gradient; in other words, the effect of the acquisition geometry and
geometrical spreading has been compensated for by the network.

Figure 4: The normalized training loss for the first three iterations of FWI.

To demonstrate the effectiveness of our approach, we compare its speed of
convergence with that of a gradient-based method with Barzilai-Borwein (BB)
update Barzilai and Borwein (1988) and a quasi-Newton method, namely L-BFGS
Liu and Nocedal (1989). The BB method is a two-point step size gradient method
used to approximate the quasi-Newton method. More precisely, the Hessian is
approximated with a scalar using the difference between two evaluations of the
gradient and two solutions

$\alpha_{k}^{long}=\frac{\Delta x\cdot\Delta x}{\Delta x\cdot\Delta g}$ (17)
$\alpha_{k}^{short}=\frac{\Delta x\cdot\Delta g}{\Delta g\cdot\Delta g}$ (18)

where $\Delta x$ and $\Delta g$ are the differences between two solutions and
two gradients, respectively. In our implementation we used $\alpha^{short}$.
For nonlinear problems, the BB method tends to converge faster than the
conjugate-gradient methods as demonstrated in Fletcher (2005).
In order to provide a fair comparison in terms of computational cost for the
different algorithms, the data misfit loss is plotted as function of the the
number of wave equation solves (instead of the iteration itself). In our
approach, in fact, every FWI iteration requires solving the wave equation
twice more than in the conventional approaches. Alongside inspecting the data
misfit loss, we also assess the quality of the updated model with respect to
the true model (whenever possible) using the structural similarity index
measure (SSIM) Wang et al. (2004); this choice is dictated by the fact that a
small vertical (or horizontal) shift in the predicted model could lead to a
small element-wise loss (e.g., mean squared error), even if the model closely
resembles the true one.
As shown in Figure 5 for the data loss and the SSIM index, the rate of
convergence as function of the number of wave equations solved is higher for
the proposed approach than for the BB and L-BFGS methods. Looking at the
models after the first iteration of FWI, we observe that the predicted model
from our method (Figure 6b) is already of higher quality than that obtained
from the BB method (Figure 6a). The shallow layers are better defined, and
even the deeper parts of the model (at a depth of 3000 m or greater) are
already updated. Moreover, looking at the final predicted models for the BB,
L-BFGS, and $\delta m$ methods, respectively (Figure 7), it is clear that our
method managed to better update the deeper part of the subsurface. This is
expected when the Hessian is included in the model update as of the case with
L-BFGS and $\delta m$, but $\delta m$ resolved the deeper part better when
comparing the two approaches to the true model. Moreover, the update in the
deeper part with $\delta m$ did not overestimate the width of the high
velocity layers looking at the two ridges in the right and left of the deeper
part of the model (depicted in dark red), for example.

Figure 5: Data loss and SSIM index for FWI on Marmousi using the BB method,
L-BFGS, and $\delta m$. Figure 6: First model update using a) the BB method,
and b) $\delta m$.

Figure 7: Final updated models: a) using the BB method, b) using L-BFGS, c)
using $\delta m$

, and d) the true model for comparison.

###  3.2 Volve synthetic

Figure 8: a) The true Volve synthetic model and b) its smoothed version. The
arrays of sources and receivers are shown in white and red, respectively.
Figure 9: Data loss and SSIM index for FWI on the Volve synthetic dataset
using the BB method, L-BFGS, and $\delta m$. Figure 10: Third model update
using a) the BB method, and b) $\delta m$ for Volve synthetic. Figure 11:
Final model update: a) using the BB method, b) using L-BFGS, and c) using
$\delta m$.

In our second example, we consider a synthetic dataset created to closely
resemble the openly available Volve field data. This dataset has been created
in Ravasi et al. (2022) as follows: first, a high-resolution acoustic
impedance model is obtained by applying post-stack seismic inversion to a 2D
section of the migrated Volve field data. Next, the inverted acoustic
impedance model is converted into a velocity model with the aid of well logs
and used to simulate a synthetic dataset with a 20 Hz Ricker wavelet. An array
of 110 sources, spaced 50 m apart and positioned at a depth of 6 m, along with
an array of 180 receivers, spaced 25 m apart at depths ranging from 86 m to 99
m, are used to model the data. Figure 8 shows the inverted velocity model,
which we refer to here as the true Volve synthetic velocity model and its
smoothed version used as an initial model for FWI.

For the Volve synthetic, which is more challenging given its more realistic
data and model, and larger model size, we utilized the data up to 7 Hz and we
estimated the wavelet from it. Inspecting the data loss curves in Figure 9a,
our approach converges at a similar rate as the BB method; however the SSIM
curves show that our method’s estimated model is of a superior quality
compared to that of the BB method (Figure 9b); this is especially the case for
later iterations where the BB model seems to be stagnant, whilst our model
continues to improve. Even though L-BFGS reached a similar data residual to
the BB and $\delta m$ methods at the end, the SSIM shows that the inverted
model started to degrade half-way through the iterations. Also, looking at the
updated models at the third (we chose the third iteration rather than the
first because the difference is easier to spot in the third iteration.) and
last iterations (Figures 10 and 11) using both approaches confirms that our
approach resulted in a better inversion since in our approach the effect of
the limited aperture is less pronounced (if any), and between 2200 to 2700 m
depth the velocity estimation is better in general. In addition, our approach
is not just a scalar approximation for the Hessian (like the BB step size)
which is evident by the correction for limited aperture in the acquisition as
shown in Figure 11c.

###  3.3 Volve field data

Figure 12: a) the velocity model obtained for Volve using tomography, and b)
its smoothed version. The arrays of sources and receivers are shown in red.
Figure 13: Data loss for FWI on the Volve dataset a) using the BB method and
$\delta m$, and b) using L-BFGS. Figure 14: First model update: a) using the
BB method, and b) using $\delta m$ for Volve. Figure 15: Final model update:
a) using the BB method, b) using L-BFGS, and c) using $\delta m$ for Volve.
Figure 16: RTM images computed using: a) the initial model, b) the updated
model with the BB method, c) the updated model with L-BFGS, and d) the updated
model with $\delta m$ for Volve. Figure 17: ADCIGs computed using the initial
model at 5.0, 5.5, 6.0, 6.5, and 7.0 km in the $x$ direction for Volve. Figure
18: ADCIGs computed using the updated model with the BB method at 5.0, 5.5,
6.0, 6.5, and 7.0 km in the $x$ directionfor Volve. Figure 19: ADCIGs computed
using the updated model with L-BFGS at 5.0, 5.5, 6.0, 6.5, and 7.0 km in the
$x$ direction for Volve. Figure 20: ADCIGs computed using the updated model
with $\delta m$ at 5.0, 5.5, 6.0, 6.5, and 7.0 km in the $x$ direction for
Volve. Figure 21: Observed data compared to simulated data computed with the
model inverted a) using the BB method, b) using L-BFGS, and c) using $\delta
m$. First and third quarters of each panel are the simulated data while the
second and fourth quarters are the observed data.

In 2008, the production of oil commenced in the Volve oil field located in the
Norwegian North Sea and continued until its shutdown in 2016. From the 3D OBC
survey acquired by Statoil in 2010 for Volve, we selected a 2D line to test
our algorithm. Since we removed multiples by multi-dimensional deconvolution
(MDD), the processed data has now co-located sources and receivers at the
seafloor. Thus, this 2D line now contains 180 sources and receivers,
positioned at depths ranging from 86.8 m to 99.2 m, with a spacing of 25 m. In
addition, we applied other processing methods such as denoising, vector-
fidelity corrections, anti-aliasing filtering, and a scaling by $\sqrt{t}$ to
manage the conversion from 3D to 2D due to geometrical spreading. For a
detailed description of the processing sequence, the reader is referred to
Ravasi et al. (2015, 2016). Figure 12 shows a 2D section of the tomographic
velocity model provided as part of the open Volve dataset along the line of
sources and receivers selected for this study in red, and its smoothed version
used as starting model for FWI.

Working with field data is more challenging due to factors such as unknown
wavelet, insufficient frequency content in the data, different physics from
the one used in the numerical modeling which leads to events in the data that
cannot be matched/reproduced, etc. Thus, after pre-processing the data to be
suitable for seismic processing and imaging as described earlier, we
approximated the shape and scale of the wavelet from the data and ran multi-
scale FWI for 4, 7, and 10 Hz maximum frequency as commonly done in practical
applications of FWI. Figure 13a shows the normalized data loss for the three
frequencies from lowest to highest for the BB method and $\delta m$ while
Figure 13b is for L-BFGS. Since we are plotting three different runs of FWI
and we have less control on the number of iteration when running the L-BFGS
optimization, we plotted the curves for L-BFGS in a separate plot. For 4Hz and
10 Hz our approach converges faster than the BB method, whilst for 7Hz it
seems to be slightly slower. In contrast, L-BFGS takes more time to catch up
with the other two methods for 4 Hz and does not perform as well in the other
two frequency scales. Note that large increases in data loss when the
frequency range is changed are expected because there is a jump in the
frequency content of the data which is not represented in the updated velocity
model. Figure 14 shows the first updated model for the BB and $\delta m$
methods: it is clear that our approach did not only provide a better update
for the shallow part of the model, but it also was able to update the deep
layers. Similarly, when looking at the final updated models in Figure 15 for
the three approaches, we can conclude that our approach has led to more
significant changes in the deeper layers, whilst at the same time the
velocities of the shallower layers are not overestimated as in the BB method
and L-BFGS.

To validate our observations, we finally performed imaging using both the
initial and the three estimated velocity models (Figure 16). The image
produced using the velocity model from the BB method presents an artifact
between around $2-2.5$ km depth (red arrow in Figure 16b), which does not
appear in the other three images. Also, the image produced from L-BFGS (Figure
16c) is of lower quality as the reflectors in the deeper part are less
focused. Most likely, this is because L-BFGS was more aggressive in its
velocity update. Lastly, comparing the image produced using our velocity model
with the one obtained from the initial velocity model, we notice that the
former presents more continuous reflectors (red arrow in Figure 16d);
moreover, the reflector at around $2.5$ km has been slightly shifted downward.
In addition, Figures 17, 18, 19, and 20 show the angle-domain common image
gathers (ADCIGs) produced using the same three velocity models used for
imaging. The image gathers produced from the initial model are in general
flat; this suggests that the initial model is already good in terms of
kinematic which is because the tomographic model is highly accurate J. Szydlik
et al. (2006), but we smoothed it to mimic the scenario where such a model is
not available. However, the angle gathers for the BB method are curving down
as pointed out by the red arrows in Figure 18 (i.e., the updated velocity
model is not better). On the other hand, the angle gathers produced with the
velocity model estimated by our approach show similar flatness as the ones
from the initial model but with more focused reflectors and better continuity,
indicated by red arrows in Figure 20. As for L-BFGS, like the RTM image, the
gathers are not focused and not flat, see the red arrow in Figure 19 for
example. Finally, Figures 21a-c compares the observed data with the simulated
data using the inverted models with the BB methos, L-BFGS, and $\delta m$,
respectively. In each panel, the first and third quarters are the simulated
data, and the second and the fourth quarters are the observed data. In
general, the BB method, and $\delta m$ data align more accurately and resemble
better the field data than L-BFGS data. However, $\delta m$ data has stronger
events in the deeper part.

##  4 Discussion

The novelty of this work stems from the fact that we embed the concept of
approximating the inverse of the Hessian into a neural network, which is
trained to map a Hessian-adjusted gradient to the FWI gradient. The effect of
the Gauss-Newton Hessian is added to the FWI gradient by means of Born
modeling and its adjoint, and the role of the network is that of approximating
the inverse of the Hessian, with the goal of producing an improved update
vector at each iteration of FWI. Since we expect the inverse of the Hessian to
slowly change from iteration to iteration, we leverage the ability of neural
networks to learn a task and adapt quickly to a similar, but slightly
different task via transfer learning. This would not be the case for standard
approaches with matching filters, as they would have to be performed from
scratch at every FWI iteration. Furthermore, since the inverse of the Hessian
can be considered as a deblurring operator that can increase the resolution of
images, the use of neural networks to approximate it is attractive due to
their success in super-resolution applications.

Despite the aforementioned advantages, our method comes also with some
shortcomings. First, the Hessian that links the gradient with the doubly
migrated image based on Born modeling and its adjoint is an approximation of
the FWI Hessian. Thus, the resulting approximation of the inverse of the
Hessian most likely does not provide as fast convergence (which is, also, the
case with most matching filters methods) as the actual inverse of the FWI
Hessian, even though we expect this to be better than having no Hessian at
all. Second, the bandwidth of the signal we are attempting to recover (the
full FWI update) is larger than what the network is provided within training.
That is due to the fact that the FWI gradient and the doubly migrated image
are ’convolved‘ by one and two wavelets, respectively. Thus, the doubly
migrated image has the same frequency support of the FWI gradient or even
less, and both have less support than the sought after scattering potential.
Thus, even though the network has learned to produce an enhanced version of
the gradient, there is no guarantee that this can extend the frequency content
beyond that present in the training data. Third, our approach presents a
computational bottleneck as it requires to perform an additional step of
modeling-migration at each FWI iteration, as well as training or fine-tuning
the network. Applying the network to the gradient once it is trained is fast,
however the training cost is not negligible.

We have applied our method to three examples, two with synthetic datasets and
one with the Volve field dataset. In general, our approach converges faster
(or has a similar convergence rate) than the quasi-Newton methods (BB and
L-BFGS). However, as discussed in the two Volve examples, our approach
resulted in a better inverted model as assessed by the SSIM index in the
synthetic case or by the different imaging products in the field case. Field
data are contaminated by noise and unwanted events, and thus, we do not wish
to perfectly fit the data. Thus, having a neural network embedded within the
FWI iteration can have a regularization effect, and thus, we might avoid
fitting the noise in the data as supported by our Volve example.

Finally, we have observed experimentally that our approach works best when
inverting low frequency data; however, when we start FWI from higher
frequencies, our method does not perform as well as when starting from low
frequencies. Nevertheless, our method is still applicable to the multi-scale
FWI, provided that we train the network from scratch at the start of each new
frequency band. On the other hand, even though our approximation of the
Hessian entails a cost equals to one conventional FWI iteration to obtain the
remigrated image, we have a better model update from the start as illustrated
in the previous examples. Thus, we can use our approximation of the Hessian in
the early stage of FWI to steer the optimization problem into a better
direction, and then, continue with the conventional FWI. This approach will
not only reduce the cost, but it is also more practical taking into account
that the role of inverse of the Hessian is generally more important in the
early iterations of FWI.

##  5 Conclusions

We have proposed a deep learning-based approach to approximate the inverse of
the Hessian in the context of FWI. Our method is based on obtaining a
remigrated image that is related to the gradient through the inverse of the
Hessian by means of Born modeling and its adjoint; subsequently, a neural
network is trained to learn the mapping from the remigrated image to the
gradient; the trained network is then applied to the gradient to enhance its
resolution and compensate for illumination effects. Numerical results on three
different datasets show that the enhanced gradients can lead to an FWI
algorithm with faster convergence and improved inverted models.

## Acknowledgments

The author thanks KAUST and the DeepWave Consortium sponsors for supporting
this research, as well as Equinor and the Volve license partners for releasing
the field dataset.

## References

  * Lailly and Bednar [1983] Patrick Lailly and J Bednar.  The seismic inverse problem as a sequence of before stack migrations.  1983\.
  * Tarantola [1984] Albert Tarantola.  Inversion of seismic reflection data in the acoustic approximation.  _Geophysics_ , 49(8):1259–1266, 1984.
  * Bording et al. [1987] R Phillip Bording, Adam Gersztenkorn, Larry R Lines, John A Scales, and Sven Treitel.  Applications of seismic travel-time tomography.  _Geophysical Journal International_ , 90(2):285–303, 1987.
  * Symes [2008] William W Symes.  Migration velocity analysis and waveform inversion.  _Geophysical prospecting_ , 56(6):765–790, 2008.
  * Bunks et al. [1995] Carey Bunks, Fatimetou M Saleck, S Zaleski, and G Chavent.  Multiscale seismic waveform inversion.  _Geophysics_ , 60(5):1457–1473, 1995.
  * Miller and Willsky [1996] Eric Lawrence Miller and Alan S Willsky.  A multiscale, statistically based inversion scheme for linearized inverse scattering problems.  _IEEE transactions on geoscience and remote sensing_ , 34(2):346–357, 1996.
  * Fichtner et al. [2013] Andreas Fichtner, Jeannot Trampert, Paul Cupillard, Erdinc Saygin, Tuncay Taymaz, Yann Capdeville, and Antonio Villasenor.  Multiscale full waveform inversion.  _Geophysical Journal International_ , 194(1):534–556, 2013.
  * Alkhalifah [2016] Tariq Alkhalifah.  Full-model wavenumber inversion: An emphasis on the appropriate wavenumber continuation.  _Geophysics_ , 81(3):R89–R98, 2016.
  * Zhang et al. [2019] Pan Zhang, Ru-Shan Wu, and Liguo Han.  Seismic envelope inversion based on hybrid scale separation for data with strong noises.  _Pure and Applied Geophysics_ , 176:165–188, 2019.
  * Jin and Madariaga [1994] Side Jin and Raul Madariaga.  Nonlinear velocity inversion by a two-step monte carlo method.  _Geophysics_ , 59(4):577–590, 1994.
  * Pratt et al. [1998] R Gerhard Pratt, Changsoo Shin, and GJ Hick.  Gauss–newton and full newton methods in frequency–space seismic waveform inversion.  _Geophysical journal international_ , 133(2):341–362, 1998.
  * Shin et al. [2001] Changsoo Shin, Seonghyung Jang, and Dong-Joo Min.  Improved amplitude preservation for prestack depth migration by inverse scattering theory.  _Geophysical prospecting_ , 49(5):592–606, 2001.
  * Liu and Nocedal [1989] Dong C Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  _Mathematical programming_ , 45(1-3):503–528, 1989.
  * Métivier et al. [2014] Ludovic Métivier, François Bretaudeau, Romain Brossier, Stéphane Operto, and Jean Virieux.  Full waveform inversion and the truncated newton method: quantitative imaging of complex subsurface structures.  _Geophysical Prospecting_ , 62(6):1353–1375, 2014.
  * AlTheyab et al. [2013] Abdullah AlTheyab, Xin Wang, and Gerard T Schuster.  Time-domain incomplete gauss-newton full-waveform inversion of gulf of mexico data.  In _2013 SEG Annual Meeting_. OnePetro, 2013.
  * Wu and Alkhalifah [2015] Zedong Wu and Tariq Alkhalifah.  Simultaneous inversion of the background velocity and the perturbation in full-waveform inversion.  _GEOPHYSICS_ , 80(6):R317–R329, 2015.  doi:10.1190/geo2014-0365.1.  URL https://doi.org/10.1190/geo2014-0365.1.
  * Guitton [2004] Antoine Guitton.  Amplitude and kinematic corrections of migrated images for nonunitary imaging operators.  _Geophysics_ , 69(4):1017–1024, 2004.
  * Aoki and Schuster [2009] Naoshi Aoki and Gerard T Schuster.  Fast least-squares migration with a deblurring filter.  _Geophysics_ , 74(6):WCA83–WCA93, 2009.
  * Valenciano [2008] Alejandro A Valenciano.  _Imaging by wave-equation inversion_.  PhD thesis, Citeseer, 2008.
  * Liu et al. [2022] Cewen Liu, Mengyao Sun, Nanxun Dai, Wei Wu, Yanwen Wei, Mingjie Guo, and Haohuan Fu.  Deep learning-based point-spread function deconvolution for migration image deblurring.  _Geophysics_ , 87(4):S249–S265, 2022.
  * Kumar et al. [2022] A Kumar, G Hampson, T Rayment, and T Burgess.  A deep learning inverse hessian for least-squares migration.  In _83rd EAGE Annual Conference & Exhibition_, volume 2022, pages 1–5. EAGE Publications BV, 2022.
  * Peng et al. [2023] H. Peng, I. Vasconcelos, and M. Ravasi.  Image-domain seismic inversion by deblurring with invertible recurrent inference machines.  _Geophysics_ , 89(2):1–62, 2023.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18_ , pages 234–241. Springer, 2015.
  * Brougois et al. [1990] A Brougois, Marielle Bourget, Patriek Lailly, Michel Poulet, Patrice Ricarte, and Roelof Versteeg.  Marmousi, model and data.  In _EAEG workshop-practical aspects of seismic data inversion_ , pages cp–108. European Association of Geoscientists & Engineers, 1990.
  * Barzilai and Borwein [1988] Jonathan Barzilai and Jonathan M. Borwein.  Two-Point Step Size Gradient Methods.  _IMA Journal of Numerical Analysis_ , 8(1):141–148, 01 1988.  ISSN 0272-4979.  doi:10.1093/imanum/8.1.141.  URL https://doi.org/10.1093/imanum/8.1.141.
  * Fletcher [2005] Roger Fletcher.  On the barzilai-borwein method.  In Liqun Qi, Koklay Teo, and Xiaoqi Yang, editors, _Optimization and Control with Applications_ , pages 235–256, Boston, MA, 2005. Springer US.  ISBN 978-0-387-24255-2.
  * Wang et al. [2004] Zhou Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE Transactions on Image Processing_ , 13(4):600–612, 2004.  doi:10.1109/TIP.2003.819861.
  * Ravasi et al. [2022] Matteo Ravasi, Tamin Selvan, and Nick Luiken.  Stochastic multi-dimensional deconvolution.  _IEEE Transactions on Geoscience and Remote Sensing_ , 60:1–14, 2022.  doi:10.1109/TGRS.2022.3179626.
  * Ravasi et al. [2015] M. Ravasi, I. Vasconcelos, A. Curtis, and A. Kritski.  Multi-dimensional free-surface multiple elimination and source deblending of volve obc data.  2015(1):cp-451-00802, 2015.  ISSN 2214-4609.  doi:https://doi.org/10.3997/2214-4609.201413355.  URL https://www.earthdoc.org/content/papers/10.3997/2214-4609.201413355.
  * Ravasi et al. [2016] Matteo Ravasi, Ivan Vasconcelos, Alexander Kritski, Andrew Curtis, Carlos Alberto da Costa Filho, and Giovanni Angelo Meles.  Target-oriented Marchenko imaging of a North Sea field.  _Geophysical Journal International_ , 205(1):99–104, 02 2016.  ISSN 0956-540X.  doi:10.1093/gji/ggv528.  URL https://doi.org/10.1093/gji/ggv528.
  * J. Szydlik et al. [2006] T. J. Szydlik, S. Way, P. Smith, L. Aamodt, and C. Friedrich.  3d pp/ps prestack depth migration on the volve field.  art. cp-2-00185, 2006.  ISSN 2214-4609.  doi:https://doi.org/10.3997/2214-4609.201402177.  URL https://www.earthdoc.org/content/papers/10.3997/2214-4609.201402177.
  * Ravasi and Vasconcelos [2020] Matteo Ravasi and Ivan Vasconcelos.  Pylops—a linear-operator python library for scalable algebra and optimization.  _SoftwareX_ , 11:100361, 2020.
  * Richardson [2023] Alan Richardson.  _Deepwave_ , September 2023.  URL https://doi.org/10.5281/zenodo.8381177.

## Appendix: Network architecture

Table 1: Summary of the network architecture Layer | Operation | Channels
---|---|---
Input Conv | 2 $\times$ Conv2d + LeakyReLU | 128
Downsampling 1 | AvgPool2d + 2 $\times$ Conv2d + LeakyReLU | 256
Downsampling 2 | AvgPool2d + 2 $\times$ Conv2d + LeakyReLU | 512
Downsampling 3 | AvgPool2d + 2 $\times$ Conv2d + LeakyReLU | 1024
Downsampling 4 | AvgPool2d + 2 $\times$ Conv2d + LeakyReLU | 2048
Upsampling 1 | ConvTranspose2d + 2 $\times$ Conv2d + LeakyReLU | 1024
Upsampling 2 | ConvTranspose2d + 2 $\times$ Conv2d + LeakyReLU | 512
Upsampling 3 | ConvTranspose2d + 2 $\times$ Conv2d + LeakyReLU | 256
Upsampling 4 | ConvTranspose2d + 2 $\times$ Conv2d + LeakyReLU | 128
Output Conv | Conv2d | 1