# Bi-directional Loop Closure for Visual SLAM

Ihtisham Ali, Sari Peltonen, and Atanas Gotchev Manuscript received xxx ;
revised xxx.Ihtisham Ali, Sari Peltonen, and Atanas Gotchev are affiliated
with Faculty of Information Technology and Communication Sciences, Tampere
University, Finland.Corresponding author: Ihtisham Ali (e-mail:
<EMAIL_ADDRESS>ihtishamaliktk@gmail.com).

###### Abstract

A key functional block of visual navigation system for intelligent autonomous
vehicles is Loop Closure detection and subsequent relocalisation. State-of-
the-Art methods still approach the problem as uni-directional along the
direction of the previous motion. As a result, most of the methods fail in the
absence of a significantly similar overlap of perspectives. In this study, we
propose an approach for bi-directional loop closure. This will, for the first
time, provide us with the capability to relocalize to a location even when
traveling in the opposite direction, thus significantly reducing long-term
odometry drift in the absence of a direct loop. We present a technique to
select training data from large datasets in order to make them usable for the
bi-directional problem. The data is used to train and validate two different
CNN architectures for loop closure detection and subsequent regression of
6-DOF camera pose between the views in an end-to-end manner. The outcome packs
a considerable impact and aids significantly to real-world scenarios that do
not offer direct loop closure opportunities. We provide a rigorous empirical
comparison against other established approaches and evaluate our method on
both outdoor and indoor data from the FinnForest dataset and PennCOSYVIO
dataset.

###### Index Terms:

Bi-directional, loop closure, pose estimation, siamese network, deep learning

##  I Introduction

Inffering where you are on a map of your local world is a core problem of
mobile robotics, navigation, and augmented reality. It refers to the process
of recognizing a previously visited place and determining your current pose
w.r.t the previous pose from the visual scene representation. This problem is
widely known as the lost or kidnapped robot. Loop closure is essential to
attain robust navigation in any intelligent transportation system as it aids
in significantly reducing the accumulated errors during visual navigation [1].
Traditionally, loop closure is detected in an environment that has been
previously viewed from a similar perspective e.g., a vehicle traveling toward
the north passes by the same location moving in the same direction. Such a
configuration maximizes the chances of place recognition. Both classical
methods, such as Binary feature descriptors in conjunction with Bag-of-Words
(BoWs), and deep learning-based approaches can be used to tackle this problem,
with the latter approach exhibiting better performance after the advent of
modern Convolutional Neural Network (CNN) architectures. In this paper, we
propose an approach to expand the capability of loop closure detection methods
towards bi-directional problems. Our proposed approach is able to correctly
recognize previously visited places and find the relative pose irrespective of
the direction of motion of the vehicle. Due to the recent success of CNNs in
the fields of image classiﬁcation [2, 3], semantic segmentation [4, 5], and
image retrieval [6, 7], many researchers opt to use them to predict camera
pose from visual data [8, 9]. In some studies, the learned deep feature
vectors are directly used for scene recognition following the concepts of
image retrieval. However, we show that decoupling the two tasks can result in
higher accuracy and better generalization.

Firstly, we introduce a novel automated technique to leverage the use of
existing large datasets for training CNNs towards the task of bi-directional
loop closure. This is essential since acquiring new data every time requires
considerable time and resources. Then, we present two CNN based networks that
are assigned to first identify potential candidates for loop closure between
the query and database and subsequently regress the pose between the matched
candidates. Our tests show that the models exhibit a successful learning
pattern, and the feature vectors can be used for scene recognition and pose
regression upon slight alterations between the networks. Moreover, the
generalization is successfully tested on unseen data thus exhibiting strong
comprehension of the visual cues by the model and not just scene memorization.

The article is organized as follows: In Section 2, we discuss the relevant
studies from the literature, their impact, and their limitations. In Section
3, we give out the system overview and formulate the proposed approaches for
place recognition and pose regression tasks. Moreover, we also discuss in
detail, the data preparation steps involved in this study for the task at
hand. Section 4 provides the experimental results and a comprehensive
comparison with other well-established methods on two challenging datasets.
Finally, Section 5 concludes the article.

##  II Related Works

Traditionally, the place recognition problem has been approached in a similar
manner as the Image retrieval problem. In general, a query image, whose
location needs to be estimated, is compared against a large geo-tagged
database of images from previous visits. Each image is represented as an
aggregate of numerous local invariant features. The state-of-the-art method
still relies on feature detectors and descriptors such as SIFT [10], ORB [11],
SURF [12], etc., that are used to extract local information from an image and
accumulated into a single feature vector for an entire image using encoded
representation through methods such as bag-of-visual-words [13], vector of
locally aggregated descriptors (VLAD) [14] or Fisher vector [15]. Fisher
Vector adopts the Gaussian mixture model (GMM) to build a visual word
dictionary matching. As a result, Fisher Vector encodes more image information
than BoW. and at times outperforms BoW in some computer vision tasks. In
contrast, VLAD is a simplification of Fisher Vector and provides a trade-off
between performance and computational efficiency. In most cases, VLAD performs
similarly to Fisher with better efficiency. These methods aid in compressing
the image representation and subsequent efficient retrieval of a match from
the database [16]. A popular approach based on the aforementioned concepts is
FAB-MAP [17] which learns a generative model for the BoW data. The model
observes and learns the co-occurrence of appearance words from common objects
that are likely to appear or disappear together thus providing valuable
probabilistic information. However, FAB-MAP proves to be computationally
expensive due to its complex methodology for image description and matching.

On the other hand, recent studies have shown that the deployment of CNNs
results in significant improvement in accuracy and reduction in complexity.
The models trained on very large datasets significantly outperform the local
descriptors such as SIFT in a variety of applications such as object and scene
recognition [18]. McManus et al. [19] proposed to learn features from image
patches and called them scene signatures. These scene signatures were for
matching and retrieving scenes under varying appearance changes. However, the
approach required a considerably more careful training phase with data of the
test environment under all possible environmental conditions. Some studies
directly opt for using the intermediate representations which are learned
using object recognition dataset and use them for scene identification [20,
21]. Sunderhauf et al. [21] propose the use of features from intermediate
layers to form a descriptor for matching. Features from higher layers of a CNN
encode semantic information about the place while features from the lower
layer encode more descriptive information about the geometry of the scene. The
authors experiment with varying combinations of these descriptions and attempt
to find the nearest neighbor based on the cosine distance between the feature
vectors of the query and the database.

It is noteworthy that all the aforementioned studies targeted uni-directional
or traditional loop closure cases. The closest work done to bi-directional
loop closure is [22]. The authors claim to target the problem of bi-
directional loop closure in panoramic images. In our humble opinion, the use
of panoramic images diminished the complications of the problem by providing
roughly similar views to a uni-directional case. The panoramas are captured in
an enclosed structural environment with a circular trajectory. As a result,
the motion in the reverse direction captures a considerable overlap of the
forward motion scenes with some spatial offset in images and only marginal
difference in perspectives. This can also be observed from the illustrations
in the study which exhibit only spatial changes in the scene. Moreover, the
study also reports that traditional methods such as FAB-MAP fail to close the
loop in practice even in these panoramic images. To the best of our knowledge,
this is the only study that targets place recognition and loop closure while
moving in the opposite direction and high perspective change. In our work, we
propose the use of a CNN topped with a neural network implementation of VLAD
known as NetVLAD [23]. The model learns to recognize places in an end-to-end
manner on the training data specifically prepared for the bi-directional loop
closure task. Moreover, we also present an approach to estimate the relative
pose between the query and the target image from the database, thus completing
the chain of tasks required for comprehensive relocalisation. The 6-DoF pose
regression in the proposed bi-directional case is significantly more
challenging compared to traditional cases and is performed with an independent
model. We employ two public datasets namely FinnForest Dataset [24] and
PennCOSYVIO dataset [25] to conduct our tests for place recognition and pose
regression.

##  III System Overview

In this section, we present the overall pipeline of the proposed approach for
localisation. The approach constitutes of two modules: a siamese CNN network
[26] with triplet structure for maximizing similarity learning and a bi-input
siamese model for 6 DoF relative camera pose regression. The overall pipeline
is shown in Figure 1. Initially, images are used from the database to train
the siamese with a triplet structure for place recognition tasks by learning
to identify maximum similarity. Each trained branch of the network is
essentially a feature encoder and the extracted feature vectors can be
employed to identify matches from the database of images that are nearest
neighbors (NN) to the query image in the feature space. Afterward, samples
with true positive matches are fed into the independently trained network for
pose regression to estimate the relative pose between true matches. The
processes are comprehensively explained in the following sections.

Figure 1:  Illustration of the overall system pipeline. A siamese network
constituting a VGG-16 base model topped with a VLAD pooling layer is used to
learn similarity in the scenes using a triplet loss. Once the training process
is completed, we employ the branch as a descriptor to compute feature
representations of database and query images for image retrieval towards
localisation. The pose regression network (lower) is independently trained to
directly regress the 6-DoF relative camera poses between the query and the
retrieved match from the database.

###  III-A Place Recognition

Feature Encoding As stated earlier, we adopt a siamese approach for the task
of place recognition whose illustration can be found in Figure 1. The network
constitutes of a base CNN model that takes three inputs namely query image
sample $I_{q}$, positive image sample $I_{p}$, and negative image sample
$I_{n}$ from the database $I_{D}$. These input images are pre-processed based
on the prescribed pre-processing technique adopted for the base model. Here,
we take VGG-16 as an example which will be mainly used in our work; however,
we also provide results with other base models for comparison later in the
study. VGG-16 takes an input image of 224x224 pixels and propagates it through
five sets of convolution and pooling layers, where the layers are connected
through Rectified Linear Unit (ReLU) as an activation function. Each layer in
the network learns a further abstraction of the input data with the highest-
level abstraction residing towards the last layers. The structure is essential
since features from the higher layers of the CNN hierarchy encode abstract
semantic information about the scene, while features from intermediate and
lower layers encode finer details from the image such as a change in
appearance and structure [21].

The outputs of the base model are normalized and fed into the neural network
form of VLAD descriptor known as NetVLAD pooling layer [23]. Essentially, VLAD
encodes information about the statistics of local descriptors aggregated over
an image in the form of feature distance from a cluster centre. For $N$
$D$-dimensional local image descriptors $\vec{x_{i}}$ as input, and $K$
cluster centres (visual words) $c_{k}$ as VLAD parameter, the output VLAD
image representation $V$ is $K\times D$-dimensional. The L2-Normalized vector
form of $V$ with elements $(j,k)$ is

$V(j,k)=\sum_{i=1}^{N}a_{k}(\vec{x_{i}})(x_{i}(j)-c_{k}(j)).$ (1)

Here, $x_{i}(j)$ and $c_{k}(j)$ are the j-th dimensions of the i-th descriptor
and k-th cluster centre, respectively. $a_{k}(\vec{x_{i}})$ indicates whether
the descriptor $\vec{x_{i}}$ belongs to the k-th visual word, i.e. it is 1 if
$\vec{x_{i}}$ belongs to the cluster $c_{k}$ and 0 otherwise.

To develop a layer reactive to training via backpropagation, it is required
that the layer’s operation is differentiable with respect to all its
parameters and the input. The original relation lacks this differentiation due
to the binary nature of $a_{k}(\vec{x_{i}})$. To overcome this issue NetVLAD
re-writes the original relation as:

$V(j,k)=\sum_{i=1}^{N}\frac{e^{w_{k}^{T}\vec{x_{i}}+b_{k}}}{\sum_{k{{}^{\prime}}}e^{w_{k}{{}^{\prime}}^{T}\vec{x_{i}}+b_{k}{{}^{\prime}}}}((x_{i}(j)-c_{k}(j)).$
(2)

where $w_{k}$ and $b_{k}$ are sets of trainable parameters for each cluster k
which are learned in an end-to-end manner during training. Conceptually, the
weight that the descriptor $\vec{x_{i}}$ is assigned to the cluster $c_{k}$
proportional to their proximity. Moreover, the relative proximity to other
cluster centres also plays a part in the relation.

For our study, we found emperically that 64 clusters and 512-dimensional VGG16
backbone work effectively for the localisation task. The NetVLAD feature
vector dimension becomes $512\times 64=32,768$. We further extract principle
components through a convolution block and retrieve the encoded description as
a normalized feature vector.

Loss Function The similarity in an image is learned by employing a triplet
loss over the output of each branch of the triplet siamese. For training, we
gather training sample set $S$ such that

$S=\\{(s_{i},{s_{i}}^{+},{s_{i}}^{-})|(s_{i},{s_{i}}^{+}\in
S^{+});(s_{i},{s_{i}}^{-}\in S^{-}),i=1,...,M\\}.$ (3)

Here, $S^{+}$ refers to the set of relevant image pairs, $S^{-}$ refers to
negative image pairs, and M indicates the span of the entire triplet set. The
triplet loss is then given as

$\ell(s_{i},{s_{i}}^{+},{s_{i}}^{-})=\\\
max(0,m+\left\|f(s_{i})-f({s_{i}}^{+})\right\|-\left\|f(s_{i})-f({s_{i}}^{-})\right\|)$

Here, margin $m$ is a scaler that defines an offset between positive and
negative pairs, and $f(.)$ is an embedding of an image sample. The global loss
over all triplet samples is given as

$L=\sum\limits_{(s_{i},{s_{i}}^{+},{s_{i}}^{-})\in\textit{{S}}}\ell(s_{i},{s_{i}}^{+},{s_{i}}^{-})$
(4)

Retrieving the nearest neighbours To retrieve a potential match for a query
image from the database of images, both images must have a suitable
representation before comparison. In the proposed case, we use one branch of
the fine-tuned network as a feature extractor to encode the query and database
images. This enables us to have the representation in the same embedding space
(i.e. 4096-dimensional feature vectors, see Figure 1). In the experimentation
section, we will use other methods of feature extraction to encode our images
for the sake of comparison. Finally, the top N-ranked database images,
$d=(d_{n}|d_{n}\in D,n=1...N)$ are selected as the nearest neighbors to the
query image based on the squared Euclidean distance in the embedding space. It
is important to note that a query image might have one, many, or no match in
the database as it depends on the number of keyframes generated during earlier
exploration of the environment.

Neighbour Confidence Sharing Place recognition is a critical task for loop
closure in visual SLAM. The problem can become significantly more challenging
when the environment contains repetitive textures even for distinct locations
such as in forests and large open areas. It is often the case that wrong
matches are generated due to similar semantics of different scenes. To
overcome this problem, we propose a confidence sharing scheme where the
confidence of the previously localized points is propagated to their neighbors
in a causal manner. In our case, we consider three neighbors for sharing
confidence. We incorporate the traveled distance between the neighbors in
order to ascertain the sanity of a potential match for a query point. A new
query point has a valid match if (1) it has a considerable match score (in
embedding space from the model) with an image from the database and (2) it has
nearby localized neighbors that agree in distance traveled with the estimates
from odometry. If a new neighbor is found far away from a nearby localized
neighbor and the odometry estimates run in favor of the previously localized
neighbor, then the new match is discarded as a possible wrong match.

###  III-B Pose Regression

The pose estimator block composes of a VGG-based siamese architecture that
takes two monocular images as raw input and predicts a 6-DoF relative
transformation between the poses for those specific inputs. The siamese
regression block is shown in Figure 1. The shared weights are initialized with
a network pre-trained for large-scale place classiﬁcation task [27] using the
Places 365 dataset, and later ﬁne-tuned for the relative pose estimation task
as described below. The output of each branch is vectorized and combined into
a single encoded description. The relative pose is regressed by passing the
feature vector through three fully connected (FC) layers activated through
Leaky ReLu functions and intermediate dropout layers for regularization, as
shown in Figure 1. The final FC layer gives out the predicted relative pose.
Different studies adopt different representations for the angles. In study
[28], the authors use deeply learned key points to estimate the Fundamental
matrix between the two images using two CNN modules. In a similar study titled
UnDeepVO [29], the authors opt for direct image alignment in conjunction with
the camera intrinsic to estimate the relative pose. The relative pose is
obtained as a decoupled combination of a translation vector and rotation
vector in Euler angles.

Although Euler angles carry briefer representation compared to the fundamental
matrix, however, it suffers from discontinuities in the form of gimbal lock.
On the other hand, rotation parameterization such as rotation matrices that
lie on a manifold, their distance computation requires ﬁnding a Euclidean
embedding. In our work, we represent the angles using quaternions similar to
the work [30]. It is important to note that quaternions lie on a unit sphere,
however, during optimization/training the difference becomes so small that the
distinction between spherical distance and Euclidean distance becomes
insigniﬁcant. Therefore, to avoid obstructing the optimization with
unnecessary additional constraints, we avoid the use of spherical constraint.
Hence, the distance between two quaternions can be measured by the Euclidean
l2 norm $\left\|q_{GT}-q\right\|$. The authors of popular study PoseNet [9]
and its derivative study [30] argue to use a decoupled approach with a
weighted parameterization of the angle, with a scale factor $\beta$, to
balance the loss function

$L=\left\|\Delta t_{GT}-\Delta t\right\|^{2}+\beta\left\|\Delta q_{GT}-\Delta
q\right\|^{2}.$ (5)

During experimentation, we observed that the approach seemed cumbersome as the
vale of $\beta$ has to be manually adjusted for each dataset. The authors [9]
remark that the value of $\beta$ can lie anywhere in the range of 120 to 2000
depending on the structure and semantics of the scene [9]. To avoid this issue
we propose to discard the scale factor $\beta$ and independently scale down
the entire translation vector and quaternions to the same range during the
preprocessing step of dataset preparation. Since we wish to use an adaptation
of ReLu activation for the FC layer, it is advisable to rescale both the
quaternions and the translation vectors between $[0,1]$. This makes the
relation invariant of any scale factor for the training phase. The scale
factors can be extracted from the range of the data in the dataset and applied
using relation

$d_{scaled}=\frac{(sc_{max}-sc_{min})*(d-d_{min})}{d_{max}-d_{min}}.$ (6)

where $d$ denotes the data array and $sc$ indicates the scaler values of the
desired range for scaling. The model is trained to predict an arbitrarily
scaled version of the pose where the scale is restored in a post-processing
step after the prediction. The MSE is then given as

$L_{MSE}=\frac{1}{n}\sum_{i=1}^{n}(T_{scaled\;GT}-T_{Est.})^{2}.$ (7)

Here, $T_{scaled\;GT}$ is the pose constituting scaled $[t_{x},t_{y},t_{z}]$
and scaled $[q_{w},q_{x},q_{y},q_{z}]$. $T_{Est.}$ is a similar vector to
$T_{scaledGT}$ which is predicted by the model. At test time, a pair of images
are fed into the regression model, consisting of two branches, which directly
estimates the relative camera pose vector. Finally, the estimated quaternion
and translation vectors are scaled up to retrieve real-world values.

###  III-C Data Selection and Dataset Preparation

Deep learning approaches require a large amount of data in order to train a
model that can generalize well for all cases. Producing new data every time is
often expensive and can drive the researcher away from the problem at hand.
For this reason, many public datasets have been introduced by researchers in
their respective fields. Visual SLAM is a widely researched problem and many
datasets exist for testing purposes [31, 32, 33, 34, 35]. However, most of the
datasets do not provide bi-directional motion since they are tailored for
handling the problem from a uni-directional perspective. For our work, we
found that FinnForest dataset [24] and parts of PennCOSYVIO dataset [25] can
be used for training and testing purposes.

As indicated by the name, the FinnForest dataset provides data, for visual
odometry and SLAM, in a forest landscape. The dataset provides recordings from
four RGB cameras that are synchronized with an Inertial Measurement Unit
(IMU), and a Global navigation satellite system (GNSS). The dataset contains
sequences for odometry that are well suited for this study. Each route is
traveled from both directions within the same sequence thus providing all the
relevant data for bi-directional loop closure. The dataset is challenging for
the problem since it contains repetitive texture, unlike an urban landscape
that provides more distinct landmarks over its trajectory. We will only
utilize the data recorded in the summer conditions in our study since it
offers slightly more landmarks than the winter condition for the localisation
block.

The second dataset that we use for training and testing offers indoor data.
This was specifically chosen to check the performance of our approach in both
indoor and outdoor environments. The dataset records similar sequences with
multiple configurations and types of sensors. We found four sequences,
recorded with GoPro Hero 4 Black, that are reasonably well suited for the task
of bi-directional loop closure. The sequences include slow and fast motions
that can represent vehicular motion using a forward-facing camera. Moreover,
the data includes ground truth poses that can be used for automatic extraction
of training and testing samples in our approach. Other sequences in this
dataset are wall-facing and more suited for Structure from Motion
applications.

Both of the datasets are passed through a data preparation phase in order to
generate sequences that can be used for training our localisation and pose
regression models. We generate sub-datasets out of the original datasets and
use them for training. Since the localisation and the pose regression are to
be performed on the same scenes we can use the data generated for localisation
in the pose regression block. For simplification, we split and consider two
cases of the bi-directional localisation problem. Assume that a route is
traversed in a straight line from both directions then we have images acquired
with a camera at somewhat regular intervals from both directions for roughly
the same location, given that the camera frame rate is high enough.

Considering the forward motion case, in Figure 2, an anchor sample is acquired
at query location (green). For these anchor samples, we can obtain positive
sample pairs from nearby locations that share the perspective view. Moreover,
it is fairly safe to assume that an image acquired further away or from far
back will provide a significantly different view and can be selected as
negative samples (red) for the localisation model training. For all the sample
pairs we compute the relative ground-truth poses in the form of translation
vector and quaternion angles which are later used in the pose regression
block. We generate 6 sample triplets for each query location, however, any
number of triplets can be generated according to the needs of the task. The
straight-line route shown in Figure 2 is a simplified case, however, we should
expect irregular movements and changes in angles in real cases. The FinnForest
dataset attempts to mimic the conditions that a heavy vehicle might face in a
real forest during its operation. we expect close-by samples from a query
point that might not share the same visual information due to sharp turns,
camera jitter, motion blur, sudden overexposure, or sun flare in the camera
view. These samples can deteriorate the learning performance of the model as
it expects them to be positive samples however they no longer share the
semantic and/or geometric information with the query sample.

Figure 2:  Illustration of training data selection based on distance and
direction. The cones represent camera body placed at various locations. The
shaded area in the backward case represents potentially overlapping regions in
the perspective view

To overcome this, we leverage structure from motion to autonomously generate
training triplets by validating the previously created triplet samples. We
generate 3D world points from a query stereo pair and track the corresponding
key points among all the positive samples for that specific query sample.
Similarly, the 3D points are also propagated between the camera frames using
the 6 DoF ground-truth poses. The transformed world points are then projected
into corresponding image space. Any sample with cumulative reprojection error
(for the tracked keypoints) higher than a threshold is discarded. The
reprojection error is computed using

$\begin{split}E_{px}=||P_{s_{i}^{+}}-\Pi(K,[q_{(s_{i},s_{i}^{+})}\>,\>_{s_{i}}t^{s_{i}^{+}}]_{HT},W_{s_{i}})||_{2}^{2}.\end{split}$
(8)

This validation step gets rid of the sample pairs with high angular changes in
perspective (such as in the case when the vehicle is turning). Here, $\Pi$ is
the perspective projection function that projects the 3D points
$W=(X,Y,Z,1)^{T}$ from world frame space to image space using the camera
intrinsic $K$. The superscript $T$ indicates the transpose of a vector. The
perspective projection yields $\tilde{x}=(\tilde{u},\tilde{v},1)^{T}$ in the
image space of the camera at the pose of interest. The reprojected points
$\tilde{x}$ are compared directly against the observed/tracked 2D points ($P$)
in the corresponding sample image. The symbol $[\>]_{HT}$ indicates the
conversion from quaternion $q$ and translation vector $t$ to the homogeneous
transformation matrix. We use quaternion angle representation for the sake of
coherence.

For the backwards case, We use the same query image samples initially selected
and filtered for the forward case and attempt to find pairs for it in the
backward motion part of the sequence. In contrast to the forward case, the
backward motion case can have positive samples only ahead of the query
location (see Figure 2). All the potential samples are expected to be oriented
in the direction opposite to the camera orientation at the query point.
Moreover, the assumption is that camera poses that are slightly ahead of the
query point would share potentially more of the same scene even if from the
opposite perspective. The similarity in the scene in this small range is what
we want our model to learn and discriminate. A camera pose that is too far
ahead or at the back of the query pose would have little and no match with the
query perspective, respectively. As before the positive samples are indicated
in blue while negative samples are shown in red. The reprojection-based
verification is not possible for the backward case since the traditional
feature detector cannot detect and track features with such high perspective
changes. As a result, there are no reference key point image positions for the
reprojected 3D points.

Figure 3: Visualization of the automatic sample selection scheme for the case
of backward motion. Valid training samples within the field of view of the
anchor/query pose (green) are shown in blue while the negative samples are
shown in red. The poses shown in black indicate rejected samples in the
vicinity of the query.

A visualization of the automatic sample selection, based on the camera poses,
from a sequence of FinnForest dataset is shown in Figure 3. Samples with in
the field overview of the query pose are considered as potential candidates
for triplet grouping. Similarly, we discard very close samples since target
samples that are too close will share very little view with the query
perspective. Here, the camera poses are shown that pass the constraints set on
the field of view, distance from the query pose, and the tolerance of
orientation difference from the query pose orientation. The candidates for
positive samples are shown in blue while the negative samples are shown in
red. Additionally, some camera poses that do not pass the constraints are
visualized in black for the sake of understanding.

It is difficult to conclude from merely visually observing the data as to what
should be the minimum and maximum distance between the query and the positive
samples. To understand the relationship we follow an empirical approach and
train the model with different data distributions. In Figure 4, we explain the
effect of data distribution. The figure shows an Upto scale visualization of
the training loss (purple) and generalization gap (indigo) for the backward
case. The horizontal axis indicates the minimum valid distance i.e, the
distance from the query to the nearest sample while the vertical indicates the
maximum valid distance i.e, the distance from the query to the furthest valid
sample. We observed that the training loss and the generalization gap were
minimum for the training data when the nearest positive sample was limited to
a distance of 2 meters and the farthest sample was kept to be at 11 meters
from the query pose. This indicates that the maximum overlap is found within
this range for the given data and that anything out of this range is a
potential outlier to what the model attempts to learn. This is obviously
specific to the datasets in this study and might vary slightly depending upon
the camera and optics used for recordings.

Figure 4: Upto scale visualization of the training loss (purple) and
generalization gap (indigo) when the training data is changed for the backward
case. The horizontal axis indicates the minimum valid distance i.e, the
distance from query to the nearest sample while the vertical axis indicates
the maximum valid distance i.e, the distance from query to the furthest valid
sample. We indicate the best result with a green overlay.

##  IV Results and Discussion

In this section, we provide our experimental results and quantitatively
demonstrate the effectiveness of the proposed system on the FinnForest and
PennCOSYVIO datasets. To ascertain the generalization capability of our
pipeline, on data previously unseen during training, we hold out one of the
scenes in the FinnForest dataset (S5) and PennCOSYVIO dataset (C2-bs) for
evaluation and train our model on the remaining scenes. We will discuss the
results in the forthcoming subsections and compare our results with other
well-established methods. The network models were implemented with the Tensor
Flow framework using Keras API. We employ the Adam optimizer to train the
network with an early stop. The network setup preferred a small learning rate
that started from .0000001 and decreased by one-tenth every epoch.

###  IV-A Place Recognition

For testing place recognition capability, we compare our proposed method with
three other approaches to gauge the relative quality of the results. Among
these methods, VGG-FC and VGG FC-norm are variants of deep learning approaches
where we deploy two fully connected layers with dropouts applied after the
VGG-16 network. VGG FC-norm has an additional normalization layer before the
feature encoded vectors are extracted from the network. The third approach
which we term here as ORB-VLAD uses ORB feature detector and descriptor to
encode keypoints from the images and uses VLAD to further re-encode and reduce
the dimensionality of the feature vectors. The use of VLAD helps us to have a
more direct comparison with our proposed approach since we employ a variant of
VLAD known as netVLAD.

(a)

(b)

Figure 5: Precision-recall curves for bi-directional loop closures in the (a)
FinnForest dataset and (b) PennCOSYVIO dataset.

For both datasets, the localisation performance is expressed in the form of a
Precision-Recall (PR) curve. The results are shown in Figure 5. It can be
observed from the results that the proposed approach VGG16-VLAD and
VGG16-FCnorm outperform VGG16-FC and ORB-VLAD. The difference between the area
under the curve (AUC) for VGG16-VLAD and VGG16-FCnorm in case of both datasets
is almost the same. Nonetheless, we remark that VGG16-VLAD is more suited for
the task at hand. Our proposition is based on the observation that VGG16-VLAD
performs better on the FinnForest dataset which is considerably more
challenging compared to the PennCOSYVIO dataset. FinnForest dataset is
recorded over a significantly larger spatial area which has repetitive
textures and fewer discriminative landmarks. On the other hand, the
PennCOSYVIO dataset offers the same indoor scene in all sequences where the
route is the same and motion speed is slightly varied. This means we can
expect a high correlation in the training and testing data in the case of
PennCOSYVIO dataset. In contrast, the route and scene are varied in the
FinnForest dataset which will result in a lower correlation between training
and testing data, and higher data center distribution (in space that houses
encoded data clusters). Hence, we can infer that VGG16-VLAD has better
generalization capability compared to the other methods. It is important to
mention that we also tested the SURF and SIFT features in combination with
VLAD. However, the results were not included as they were poor and inhibited
the readability of the PR curve. These classical feature descriptors work well
for the uni-directional cases where the perspective does not change a lot,
however, they fail to perform well in the bi-directional cases.

To gain a better understanding of what the network has learned and what it
sees in an image, we overlay the activation maps on their corresponding images
for visual observation. The activation maps for the forward motion case are
shown in Figures 7 and 9 for FinnForest and PennCOSVIO datasets, respectively.
Similarly, the activation maps for the motion in the opposite direction (bi-
directional case) are shown in Figures 7 and 9. For the forward motion case,
it can be observed that almost the same regions are activated in the query and
positive test image pairs, which is an intuitive conclusion. In contrast, the
activation maps are flipped left-right in the backward motion case for the
query image and the positive sample pair. This effect makes sense since the
images are acquired from the opposite directions for roughly the same
location. This flip effect is dominant and easily observed from image pairs in
Figure 7.a,d and 7.b,e. Moreover, the regions closer to the camera exhibit
stronger activation compared to the regions that are far away. All these
observations are in agreement with our hypothesis mentioned in Section III-C
that motivated the study.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Activation maps over sample image triplet used for testing from
FinnForest Dataset for forward/uni-directional case, where maps in (a-c) are
for query images, (d-f) are for corresponding (column wise) positive samples.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: Activation maps over sample image triplet used for testing from
FinnForest Dataset for bi-directional motion case, where maps in (a-c) are for
query images, (d-f) are for corresponding (column wise) positive samples

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: Activation maps over sample image triplet used for testing from
PennCOSYVIO Dataset for forward/uni-directional motion case, where maps in
(a-c) are for query images, (d-f) are for corresponding (column wise) positive
samples.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: Activation maps over sample image triplet used for testing from
PennCOSYVIO Dataset for bi-directional motion case, where maps in (a-c) are
for query images, (d-f) are for corresponding (column wise) positive samples.

Figure 10:  illustration of localisation results for FinnForest dataset. The
top row shows the results for detections from the forward pass while the lower
row shows the results for localisation from the backward pass.

Figure 11:  illustration of localisation results for PennCOSYVIO dataset. The
top row shows the results for detections from the forward pass while the lower
row shows the results for localisation from the backward pass.

PR curves are a good metric for binary classiﬁcation and to understand the
overall performance of a system. However, to fully prove that our method is
capable of accurately closing loops in practice, we perform loop closure
offline in a causal manner. For this experiment, we simulate keyframe
selection for loop closure at regular intervals based on the distance
traveled. A loop closure candidate is detected, if the score of the query and
a keyframe from the database is above a priori threshold $\tau$ and if it has
passed the confidence sharing criterion among its neighboring localized
keyframes. We eliminate the most recent images from the search space and wait
until the database is large enough to start loop detection. The values of
$\tau$ for the experiments shown in Figures 11 and 11 were selected from
Figure 5 such that the recall rate is maximized with good precision returns. A
slightly higher threshold was selected to illustrate all the possible outcomes
in the experiments. As illustrated, a search can result in no match, a single
unique match, multiple valid matches, a valid match with one or many false
positives, or an invalid match. Multiple matches can be observed when the
query images are acquired in open spaces and the scene does not change much
among subsequent keyframes. In some cases, a keyframe can find a true and
false positive at the same time due to visual similarity between multiple
keyframes. This is more apparent in the FinnForest case where we observe quite
frequent repetitive textures in the trees and on the road. Nonetheless, these
false positives and invalid matches can be removed with the confidence-sharing
approach between the neighboring keyframes.

It can be observed that the matches for the forward motion cases are
significantly higher than the reverse case for the FinnForest. This makes
sense since the forward motion provides more opportunities to observe similar
scenes. In contrast, it is difficult to localize in the case of motion in the
opposite direction for the FinnForest dataset as most of the scenes do not
provide distinct landmarks within that short spatial window of observation.
Nonetheless, the model was able to successfully identify the absence of
matches and found enough matches that could be used to significantly improve
odometry results. Even one valid match is enough to drastically improve
odometry results and mitigate the errors due to drift that are accumulated in
the odometry results.

On the other hand, the model works very well on the PennCOSYVIO dataset since
the indoor constrained environment provides distinct visual landmarks that are
specific to their corresponding locations. Extrapolating on this observation,
we can postulate that the approach would be effective in outdoor scene
recorded within a city environment that provides distinct landmarks for bi-
directional localisation. As mentioned earlier, we were not able to use the
existing visual odometry datasets that were recorded in an urban environment
as they are tailored for uni-directional odometry and SLAM purposes.

###  IV-B Pose Regression

Earlier in Section III-B, we discussed the pose regression block in detail and
stressed the necessity of pose estimation through end-to-end learning. Here,
we provide the experimental results of the proposed approach and compare the
results against alternative approaches. The experimental results are tabulated
in Tables I and II for FinnForest and PennCOSYVIO datasets, respectively. For
comparison, we test the network by replacing the base model (VGG16) with
Resnet50. Moreover, the relative impact of weight initialization is also
studied by initializing the base models with weights acquired from models
previously trained on ImageNet and Places 1365 (an extension of the Places 365
dataset) for classification tasks.

For both datasets, the pose regression performance is measured as the absolute
difference between the predicted and ground truth values for the location (in
meters) and orientation (in degrees). Similar to localisation, we test the
performance of the pose regression network and state the results on individual
sequences and the combined case. The combined testing results are effectively
the average of the individual results. The number of test image samples and
the spatial extent of the area where each sequence was recorded are also
mentioned.

It can be observed that the network that has VGG16 as the base model and
initialized with the weights of Places 1365 yields the best results followed
by VGG16 initialized with ImageNet. This improvement was observed since Places
1365 incorporates visual scenes for scene classification that are quite
relevant to our localisation task. In contrast, ImageNet is a more diverse
dataset that is tailored for object classification. As a result,
initialization with Places 1365 aids our network to generalize better to
landscapes. On the other hand, Resnet50 performed poorly for both datasets. It
is important to remember that the task at hand is localisation and not visual
odometry. The pose regression is aimed at finding the relative pose between
the query and a potential match for localisation. Traditional methods fail
when we consider bi-directional cases of localisation. The results obtained
for bi-directional pose regression in this study match the performance of
other state-of-the-art approaches that are reported in studies conducted for
uni-directional loop closure [9, 30].

TABLE I: Comparison of pose estimation results from the regressor model
trained on FinnForest dataset.

Sequence |  Test Samples |  Spatial Extent (m) |  Resnet50 Imagenet |  VGG-Imagenet |  VGG-Places 1365
---|---|---|---|---|---
S1 | 2044 | 47 x 193 | 5.38m, 1.02∘ | 2.42m, 0.352∘ | 2.26m, 0.3∘
S3 | 2706 | 800 x 190 | 5.26m, 1.00∘ | 2.38m, 0.33∘ | 2.31m, 0.29∘
S4 | 3566 | 812 x 568 | 5.68m, 1.06∘ | 2.54m, 0.39∘ | 2.36m, 0.32∘
S5 | 8866 | 1826 x 1883 | 7.35m, 0.84∘ | 3.35m, 0.44∘ | 3.23m, 0.53∘
Combined | 17182 | 2633 x 2014 | 5.92m, 0.98∘ | 2.67m, 0.38∘ | 2.54m, 0.36∘

TABLE II: Comparison of pose estimation results from the regressor model
trained on PennCOSYVIO dataset.

Sequence |  Test Samples |  Spatial Extent (m) |  Resnet50 Imagenet |  VGG-Imagenet |  VGG-Places 1365
---|---|---|---|---|---
C2-af | 3361 | 144 x 36 | 3.79m, 0.71∘ | 1.51m, 0.21∘ | 1.35m, 0.22∘
C2-bs | 3330 | 144 x 36 | 3.85m, 0.72∘ | 1.51m, 0.20∘ | 1.33m, 0.21∘
C2-bf | 3090 | 144 x 36 | 3.81m, 0.73∘ | 1.49m, 0.19∘ | 1.36m, 0.22∘
C2-bs | 3375 | 144 x 36 | 5.75m, 0.80∘ | 2.20m, 0.40∘ | 1.81m, 0.26∘
Combined | 13156 | 144 x 36 | 4.3m, 0.74∘ | 1.68m, 0.25∘ | 1.46m, 0.22∘

##  V Conclusion

The article presents for the first time, to our knowledge, a learning-based
approach to solve the problem of bi-directional loop closure in monocular
images. We segregate the tasks of localisation into place identification and
pose regression and solve them in two end-to-end deep learning steps. We
demonstrate that it is indeed possible to achieve bi-directional loop closure
on monocular images by carefully posing the problem and leveraging the
training data for the networks. Moreover, we demonstrate that the networks
generalize well and aim for learning the geometric and spatial relations in
images rather than memorize the scenes/locations. This is validated by the
performance of the model on unseen data. We compare the proposed approach with
other deep learning methods and classical approaches and demonstrate superior
performance for localisation. Moreover, we provide both qualitative and
quantitative results to corroborate the claim. A natural extension of the work
would be to extend the case scenarios and test the approach with more
datasets.

## References

  * [1] X. Gao, R. Wang, N. Demmel, and D. Cremers, “Ldso: Direct sparse odometry with loop closure,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 2198–2204.
  * [2] S. Anwar, M. U. Qadir, I. Ali, S. Razzaq, B. Ahmad, K. Shah, S. A. Shah, M. T. Khan _et al._ , “Computer aided diagnosis of acoustic neuroma: a neural network perspective,” _Journal of Medical Imaging and Health Informatics_ , vol. 7, no. 2, pp. 371–377, 2017.
  * [3] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” _Advances in neural information processing systems_ , vol. 25, 2012.
  * [4] S. Hong, H. Noh, and B. Han, “Decoupled deep neural network for semi-supervised semantic segmentation,” _Advances in neural information processing systems_ , vol. 28, 2015.
  * [5] H. Noh, S. Hong, and B. Han, “Learning deconvolution network for semantic segmentation,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1520–1528.
  * [6] A. Babenko, A. Slesarev, A. Chigorin, and V. Lempitsky, “Neural codes for image retrieval,” in _European conference on computer vision_. Springer, 2014, pp. 584–599.
  * [7] A. Gordo, J. Almazán, J. Revaud, and D. Larlus, “Deep image retrieval: Learning global representations for image search,” in _European conference on computer vision_. Springer, 2016, pp. 241–257.
  * [8] A. Kendall and R. Cipolla, “Geometric loss functions for camera pose regression with deep learning,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 5974–5983.
  * [9] A. Kendall, M. Grimes, and R. Cipolla, “Posenet: A convolutional network for real-time 6-dof camera relocalization,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 2938–2946.
  * [10] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International journal of computer vision_ , vol. 60, no. 2, pp. 91–110, 2004\.
  * [11] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: An efficient alternative to sift or surf,” in _2011 International conference on computer vision_. Ieee, 2011, pp. 2564–2571.
  * [12] H. Bay, A. Ess, T. Tuytelaars, and L. Van Gool, “Speeded-up robust features (surf),” _Computer vision and image understanding_ , vol. 110, no. 3, pp. 346–359, 2008.
  * [13] J. Philbin, O. Chum, M. Isard, J. Sivic, and A. Zisserman, “Object retrieval with large vocabularies and fast spatial matching,” in _2007 IEEE conference on computer vision and pattern recognition_. IEEE, 2007, pp. 1–8.
  * [14] R. Arandjelovic and A. Zisserman, “All about vlad,” in _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2013, pp. 1578–1585.
  * [15] H. Jégou, F. Perronnin, M. Douze, J. Sánchez, P. Pérez, and C. Schmid, “Aggregating local image descriptors into compact codes,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 34, no. 9, pp. 1704–1716, 2011.
  * [16] H. Jegou, M. Douze, and C. Schmid, “Product quantization for nearest neighbor search,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 33, no. 1, pp. 117–128, 2010.
  * [17] M. Cummins and P. Newman, “Fab-map: Probabilistic localization and mapping in the space of appearance,” _The International Journal of Robotics Research_ , vol. 27, no. 6, pp. 647–665, 2008.
  * [18] A. R. Memon, H. Wang, and A. Hussain, “Loop closure detection using supervised and unsupervised deep neural networks for monocular slam systems,” _Robotics and Autonomous Systems_ , vol. 126, p. 103470, 2020.
  * [19] C. McManus, B. Upcroft, and P. Newman, “Scene signatures: Localised and point-less features for localisation,” _Robotics: Science and Systems X_ , pp. 1–9, 2014.
  * [20] Z. Chen, O. Lam, A. Jacobson, and M. Milford, “Convolutional neural network-based place recognition,” _arXiv preprint arXiv:1411.1509_ , 2014\.
  * [21] N. Sünderhauf, S. Shirazi, F. Dayoub, B. Upcroft, and M. Milford, “On the performance of convnet features for place recognition,” in _2015 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2015, pp. 4297–4304.
  * [22] R. Arroyo, P. F. Alcantarilla, L. M. Bergasa, J. J. Yebes, and S. Gámez, “Bidirectional loop closure detection on panoramas for visual navigation,” in _2014 IEEE Intelligent Vehicles Symposium Proceedings_. IEEE, 2014, pp. 1378–1383.
  * [23] R. Arandjelovic, P. Gronat, A. Torii, T. Pajdla, and J. Sivic, “Netvlad: Cnn architecture for weakly supervised place recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 5297–5307.
  * [24] I. Ali, A. Durmush, O. Suominen, J. Yli-Hietanen, S. Peltonen, J. Collin, and A. Gotchev, “Finnforest dataset: A forest landscape for visual slam,” _Robotics and Autonomous Systems_ , vol. 132, p. 103610, 2020.
  * [25] B. Pfrommer, N. Sanket, K. Daniilidis, and J. Cleveland, “Penncosyvio: A challenging visual inertial odometry benchmark,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 3847–3854.
  * [26] S. K. Roy, M. Harandi, R. Nock, and R. Hartley, “Siamese networks: The tale of two manifolds,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 3046–3055.
  * [27] B. Zhou, A. Lapedriza, A. Khosla, A. Oliva, and A. Torralba, “Places: A 10 million image database for scene recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 40, no. 6, pp. 1452–1464, 2017\.
  * [28] Y.-Y. Jau, R. Zhu, H. Su, and M. Chandraker, “Deep keypoint-based camera pose estimation with geometric constraints,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 4950–4957.
  * [29] R. Li, S. Wang, Z. Long, and D. Gu, “Undeepvo: Monocular visual odometry through unsupervised deep learning,” in _2018 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2018, pp. 7286–7291.
  * [30] Z. Laskar, I. Melekhov, S. Kalia, and J. Kannala, “Camera relocalization by computing pairwise relative poses using convolutional neural network,” in _Proceedings of the IEEE International Conference on Computer Vision Workshops_ , 2017, pp. 929–938.
  * [31] J. Engel, V. Usenko, and D. Cremers, “A photometrically calibrated benchmark for monocular visual odometry,” _arXiv preprint arXiv:1607.02555_ , 2016\.
  * [32] N. Carlevaris-Bianco, A. K. Ushani, and R. M. Eustice, “University of michigan north campus long-term vision and lidar dataset,” _The International Journal of Robotics Research_ , vol. 35, no. 9, pp. 1023–1035, 2016.
  * [33] A. Z. Zhu, D. Thakur, T. Özaslan, B. Pfrommer, V. Kumar, and K. Daniilidis, “The multivehicle stereo event camera dataset: An event camera dataset for 3d perception,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 2032–2039, 2018.
  * [34] H. Jung, Y. Oto, O. M. Mozos, Y. Iwashita, and R. Kurazume, “Multi-modal panoramic 3d outdoor datasets for place categorization,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 4545–4550.
  * [35] W. Maddern, G. Pascoe, C. Linegar, and P. Newman, “1 year, 1000 km: The oxford robotcar dataset,” _The International Journal of Robotics Research_ , vol. 36, no. 1, pp. 3–15, 2017.

##  VI Biography Section

|  Ihtisham Ali received his B.Sc. in Mechatronics Engineering from the
University of Engineering and Technology, Pakistan (2014) and his M.Sc. in
Automation Engineering from Tampere University, Finland (2017). Currently, he
is a doctoral researcher in 3D Media Group at Tampere University. He has
worked on several industrial projects pertaining to machine automation using
visual cues. His research interest is focused on computer vision and robotics
specifically object pose estimation, 3D reconstruction, and visual SLAM.
---|---
|  Sari Peltonen received the M.Sc. degree in mathematics from the University
of Tampere, Finland, in 1996. She received the Ph.D. degree in signal
processing from the Tampere University of Technology in 2000. She is currently
University Lecturer in signal processing at the Unit of Computing Sciences in
Tampere University. Her research interests include robust estimation, image
processing and tomographic image reconstruction.
---|---
|  Atanas Gotchev received his M.Sc. degrees in radio and television
engineering (1990) and applied mathematics (1992), his Ph.D. degree in
telecommunications (1996) from the Technical University of Sofia, and the
D.Sc.(Tech.) degree in information technologies from the Tampere University of
Technology (2003). He is a Professor of Signal Processing and Director of the
Centre for Immersive Visual Technologies at Tampere University. His recent
work concentrates on the algorithms for multi-sensor 3-D scene capture,
transform-domain light-field reconstruction, and Fourier analysis of 3-D
displays.
---|---


# Seismic Full-Waveform Inversion Using Deep Learning Tools and Techniques

Alan Richardson (Ausar Geophysical)

###### Abstract

I demonstrate that the conventional seismic full-waveform inversion algorithm
can be constructed as a recurrent neural network and so implemented using deep
learning software such as TensorFlow. Applying another deep learning concept,
the Adam optimizer with minibatches of data, produces quicker convergence
toward the true wave speed model on a 2D dataset than Stochastic Gradient
Descent and than the L-BFGS-B optimizer with the cost function and gradient
computed using the entire training dataset. I also show that the cost function
gradient calculation using reverse-mode automatic differentiation is the same
as that used in the adjoint state method.

##  1 Introduction

Deep learning is a method for training deep (many layer) neural networks
(DNNs). It has recently attracted attention because of successful applications
in several fields, including computer vision [8]. DNNs consist of multiple
interconnected layers that each apply transformations to their inputs. When
the output of a layer is fed back into the same layer as an input, it is said
to be a recurrent neural network (RNN). Deep learning trains DNNs by inverting
for the variables that parameterize the layer transformations. Reverse-mode
automatic differentiation, including a backpropagation stage, is used to
calculate the gradient of the chosen cost function with respect to the
variables. The variables can then be iteratively updated so that the DNN
produces the desired output.

Full-waveform inversion (FWI, [12]) is a method for deriving a model of
subsurface properties such as wave speed from seismic data recorded on the
surface. To achieve this, it performs nonlinear inversion iteratively, using
the adjoint state method, including backpropagation, to determine the gradient
[10]. The cost function is generally the sum of squares of the difference
between the recorded data and data obtained by forward propagating a source
wave through the current model.

There is clearly similarity between deep learning and FWI, but what is the
benefit of combining them? Expressing an algorithm in the form of a DNN has
the advantage of facilitating the use of deep learning software, such as
TensorFlow [1]. These allow sophisticated implementations to be quickly
developed. Once an algorithm is expressed using TensorFlow, for example, it
can be executed on CPUs and GPUs, and run in parallel across a distributed
memory computer cluster. Furthermore, it makes it easier to apply new ideas
proposed in the active deep learning field to the algorithm, such as the Adam
optimizer [7].

My hypothesis is that the seismic algorithm FWI can be expressed using DNNs
and so implemented using deep learning software, and that it may benefit from
the deep learning practice of using minibatches of data with the Adam
optimizer.

##  2 Materials and Methods

FWI requires a method for forwarding modeling – numerically propagating a
source wave through the model and recording the wavefield amplitude at
receiver locations to create synthetic data. This allows the calculation of
the cost function, the sum of the squared differences between this synthetic
data and the true recorded data. FWI also requires a method for calculating
the gradient of the cost function with respect to the wave speed model and any
other model parameters used. This allows the model to be updated in a
direction that reduces the cost function value, resulting, after multiple
iterations, in a model that produces synthetic data that matches well with the
true data.

I use TensorFlow for my implementation, but it should also be possible with
other deep learning software.

In this section, I first describe the implementation of forward modeling of
the wave equation using an RNN, and then discuss two approaches to calculating
the gradient (automatic differentiation and the adjoint state method),
followed by a description of the SEAM model that I use for experiments.

###  2.1 Forwarding modeling

Forward modeling can be implemented by repeatedly applying a sequence of
arithmetic operations that correspond to one time step of a finite difference
wave equation discretization. Sequences of arithmetic operations are directly
representable using a neural network, where, in this case, the trainable
variables are the subsurface model parameters such as wave speed. The network
is therefore an RNN, and each repeating unit corresponds to one finite
difference time step. For simplicity, I describe the 1D, constant density,
acoustic (scalar wave equation) case. Extension to higher dimensions is
straightforward. The implementation of more sophisticated wave equations is
also, if numerical simulation of the chosen wave equation can be expressed
using a sequence of operations permissible in neural networks.

In an RNN, the same operations are applied in each unit, but the data that
these operations act upon changes. The data consist of state vectors, which
are passed from one unit to the next and signify the current state of the
system, and external inputs. In addition to passing state vectors to the next
RNN unit, each unit also produces an output vector that is not passed on.

When using a second-order finite difference approximation in time, two
adjacent time steps of the wavefield need to be stored simultaneously. To
prevent unwanted reflections from the edges of the simulation domain, I use a
perfectly matched layer (PML). This requires the storage of an additional
auxiliary wavefield. The state vector of the RNN therefore consists of two
time steps of the regular wavefield, which I refer to as $u(x,t)$ and
$u(x,t-1)$, and one auxiliary wavefield, which I call $\phi(x,t)$. The
external input to each time step is the source amplitude at that time step,
$f(x,t)$, and the output is either the full wavefield at time step $t$,
$u(x,t)$, or only the wavefield amplitude at the receiver locations. This
forward modeling neural network, using the conventional RNN structure, is
depicted in Figure 1.

Figure 1: Wave equation forward modeling using an RNN. Each unit applies the
operations to propagate forward one time step, taking the state from the
previous unit (the wavefield $\mathbf{u}$ at adjacent time steps and auxiliary
wavefield $\mathbf{\phi}$) and the source amplitude $\mathbf{f}$ as inputs,
and producing updated state vectors (advanced by one time step) and the
wavefield $\mathbf{u}$ as outputs.

In my implementation, I use convolution with a finite difference kernel to
implement the spatial derivatives in the wave equation simulation. The
convolution operation is optimized in DNN software due to its popularity in
computer vision applications.

As wavefield values are only required at the receiver locations to compute the
cost function, it is possible to only use those values as the output from the
RNN units instead of the whole wavefield $u(x,t)$. In my implementation I
output the whole wavefield as it is useful for the adjoint state method
gradient calculation.

###  2.2 Automatic differentiation

As all of the operations needed to numerically simulate wave propagation are
differentiable, automatic differentiation can be used to calculate the
gradient of the cost function. It is thus only necessary to describe the
forward modeling step of the algorithm and specify the cost function; the deep
learning software then has sufficient information to calculate the gradient.

I derive the equations for automatic differentiation of the FWI cost function
in Appendix A.

###  2.3 Adjoint state method

The traditional means of calculating the gradient in FWI is to use the adjoint
state method. This involves backpropagating the difference between the forward
modeled data and the recorded data, which can be achieved using the same
forward modeling RNN described above. The equations of the adjoint state
method and automatic differentiation are in fact the same, as I show in
Appendix A. The difference is that I use the adjoint state method to refer to
manually implementing the gradient computation rather than relying on the deep
learning software’s automatic differentiation.

Calculating the gradient of the cost function with respect to the wave speed
is achieved using

$\frac{\partial
J}{\partial\mathbf{c}}=\sum_{s=1}^{N_{s}}\sum_{t=1}^{N_{t}}-\frac{2}{\mathbf{c}^{3}}\frac{\partial^{2}F\left(\mathbf{f}_{s,t}\right)}{\partial
t^{2}}F^{-1}\left(2\mathbf{r}_{s,t}\right),$ (1)

where $J$ is the cost function value, $\mathbf{c}$ is the wave speed model,
$N_{s}$ and $N_{t}$ are the number of shots and time steps in the dataset,
respectively, $\mathbf{f}_{s,t}$ contains the source amplitudes for shot $s$
as a function of time $t$, $\mathbf{r}_{s,t}$ is the difference between the
synthetic and real receiver data for shot $s$ as a function of time $t$, and
$F(\mathbf{g}_{t})$ and $F^{-1}(\mathbf{g}_{t})$ are wave propagators of the
source amplitude $\mathbf{g}_{t}$ forward and backward in time, respectively.
As $F\left(\mathbf{f}_{s,t}\right)$ is $u(x,t)$, the output from the forward
modeling RNN, and $F^{-1}\left(2\mathbf{r}_{s,t}\right)$ uses the same forward
modeling RNN but runs it backward in time and uses $2\mathbf{r}_{s,t}$ as the
source term, we already have the necessary ingredients to implement the
adjoint state method.

###  2.4 SEAM dataset

Some of the experiments that I describe in the results section use a model
derived from the SEAM Phase I model [5]. I use a 2D section of the Vp model
extracted from the 23900 m North line, covering from 7000 m to 16000 m in the
horizontal direction and from 750 m to 6500 m in the depth direction. The
extracted model is interpolated onto a grid with 100 m cell spacing. The
dataset contains shots forward modeled on this model, using a 1 Hz Ricker
wavelet as the source with a source and receiver spacing of one grid cell
along the top surface.

In deep learning it is common to reserve part of the data as a “development”
or “validation” dataset. This allows hyperparameters to be tuned with reduced
risk of overfitting. In hyperparameter selection, models are trained using the
“training” dataset (the remainder of the data, sometimes after also extracting
a “test” dataset for final model testing) with various hyperparameter choices.
The cost function values of the resulting trained models using the development
dataset are used to select between the hyperparameter configurations, and,
because this dataset was not used in training, we have greater confidence that
this hyperparameter configuration will also result in a small cost function
value on data not included in the training set; it generalizes well by not
overfitting the data.

I reserve ten randomly selected shots for the development dataset, leaving the
remaining eighty shots for training. The training shots are also shuffled so
that minibatches of data are more likely to be representative of the whole
dataset.

##  3 Results

To compare FWI implemented using deep learning concepts with a more
conventional approach, I compare the gradients calculated using automatic
differentiation and the adjoint state method, and perform inversion using the
Adam and Stochastic Gradient Descent (SGD) optimizers (with minibatches) and
the L-BFGS-B optimizer (using the entire dataset to calculate the cost and
gradient), where the hyperparameters for the Adam and SGD optimizers were
chosen using the development dataset.

The code to reproduce these results is contained in the ancillary archive.

###  3.1 Gradient comparison

The gradient of the cost function with respect to the model parameters is
needed to update the model. In deep learning applications it is typically
calculated by using automatic differentiation included in most deep learning
software. In contrast, FWI applications generally use the adjoint state method
to calculate the gradients. It is also possible (but computationally
expensive) to use the finite difference method. This involves perturbing the
model at each point and forward modeling using this perturbed model to
determine the effect the perturbation has on the cost function. To compare
these approaches, I create a random 1D model as the true model, randomly
perturb it to produce an initial model, and compute the cost function
gradient.

In Figure 2 I compare the gradient at each point in the 1D model computed
using these three approaches. As expected, the results are similar. The
adjoint state method does appear to slightly overestimate the gradients, and
both it and automatic differentiation differ from the finite difference
results at the edges of the model. These errors are likely to be due to
numerical inaccuracies, and both are sufficiently small that they are unlikely
to have a significant effect on iterative optimization.

Figure 2: The gradient of the cost function with respect to the 1D wave speed
model, computed using three approaches. The results are similar.

###  3.2 Using a development dataset for hyperparameter selection

I wish to use minibatches of data with the Adam and SGD optimizers to perform
training with my SEAM dataset. Two important hyperparameters are thus the
batch size (the number of shots in each minibatch), and the learning rate (a
number that controls the size of model updates). Randomly choosing twenty
batch size and learning rate combinations, I train using each combination with
approximately forty shots from the training dataset, and evaluate the cost
function on the resulting models using the development dataset. The results
for the Adam optimizer are plotted in Figure 3. It appears that for Adam, a
learning rate between about 50 and 150 yields the lowest final cost function
value. From these results, it is not possible to determine whether batch size
has an effect. The lowest value results from a learning rate of 45.4 and a
batch size of 2. The best combination for SGD is 733.35 and 5.

Figure 3: The cost function value (loss) on the development dataset after
training using half of the training set with different combinations of batch
size and learning rate hyperparameters using the Adam optimizer. Of the
combinations considered, a learning rate of 45.4 and a batch size of 2 give
the lowest loss.

###  3.3 Adam and SGD minibatch updates compared to L-BFGS-B

Two particularly interesting techniques commonly used in deep learning are
minibatches and the Adam optimizer [7]. Minibatches attempt to significantly
reduce training time by updating the model after only computing the cost
function and gradient on a small number of data samples. Computing the cost
function and gradient using the entire dataset, a practice sometimes used in
FWI, can ensure that each model update always reduces the cost function value,
but this guarantee may not be worth the computational cost. The Adam optimizer
uses concepts such as momentum to improve convergence and has been used
successfully in many machine learning applications.

In this test I compare the Adam and SGD optimizers, using minibatches and the
hyperparameters found in the previous section, with the L-BFGS-B optimizer
[13]. For the latter, the cost and gradient are calculated using the entire
dataset, and I constrain the range of inverted values to be between 1490 and
5000 m/s.

Figure 4: The cost function value (loss) on the development dataset declines
more rapidly when using the Adam approach compared to SGD and L-BFGS-B. Adam
is also more stable than SGD.

Figure 4 plots the cost function value calculated using the reserved
development dataset on the current model of each optimizer as a function of
the number of shot evaluations (calculations of the cost and gradient of one
shot). It shows that Adam converges significantly more quickly than L-BFGS-B.
Adam’s convergence is also quicker and more stable than SGD. In fact, the cost
function when using Adam becomes flat after approximately one pass through the
training dataset. This means that even when the additional shot evaluations
used to determine the hyperparameters for Adam are considered, it still
converges more quickly than L-BFGS-B did in this experiment.

Figure 5: The model produced by Adam is close to the true model, while SGD has
only captured shallow features.

The final models produced by Adam and SGD after 400 shot evaluations (not
including those used for hyperparameter selection), are shown in Figure 5. As
expected, given its low final cost function value, the model produced by Adam
is quite similar to the true model: the sediment areas are mostly well
matched, and parts of the high velocity salt body on the right of the model
are approximately correct. Some shallow features, such as the top of the salt
body and the sea floor, are visible in the SGD model, but deeper parts of the
model are incorrect. The model produced by L-BFGS-B after 1200 shot
evaluations is not shown as it has not changed significantly from the initial
model. This is not surprising, as it has only evaluated 16 different models,
while the Adam model was produced after evaluating 201 different models (not
including the hyperparameter selection process).

##  4 Discussion

#### Automatic differentiation

There are several advantages of using automatic differentiation to calculate
the gradient of the cost function with respect to model parameters. Foremost
among these is that it greatly reduces development time as it avoids the need
to manually implement and verify the adjoint state method. Another advantage
is that it makes it easier to try new ideas, such as experimenting with other
deep learning concepts like dropout [11]. The major disadvantage, however, is
its high memory requirement, as it stores all of the wavefield time slices in
memory. In FWI, the wavefield at the final time step can certainly be affected
by interactions of the wavefield and the model during the initial time steps,
and so backpropagation must occur over the entire sequence of time steps. It
is not typical in deep learning applications to backpropagate through
thousands of layers. As a result, automatic differentiation in deep learning
software is not usually designed to efficiently handle such situations, and
does not employ strategies common in FWI software to reduce memory
requirements. These strategies include not saving the wavefield at every time
step, applying lossy compression to the wavefield time slices that are stored,
saving the wavefields to disk instead of keeping them in memory, and
regenerating wavefields during backpropagation rather than storing them. Using
the adjoint state method to compute the gradient instead of relying on
automatic differentiation provides enough control, even when using deep
learning software, to implement some of these strategies. This means that the
adjoint state method, rather than automatic differentiation, is likely to be
necessary when applying FWI to realistically sized models.

#### Adam

Instead of waiting to check every shot before updating the model, minibatch
updates allow more rapid exploration of the model space. One of the main
advantages of the Adam optimizer is its suitability for working with
minibatches. Although it is possible to use minibatches with some more
traditional optimizers, as in Stochastic Gradient Descent [3], Adam’s use of
concepts such as momentum make it particularly well adapted to minibatches.
Momentum not only helps to improve convergence by reducing the propensity of
optimizers to struggle in narrow valleys of the cost function, but it also
tackles the problem of small minibatches not being representative of the
entire dataset by smoothing gradients over several minibatches. This is likely
to be the cause of Adam’s greater stability compared to SGD in the SEAM
experiment. A further attraction of Adam is its wide availability in deep
learning software. In contrast, TensorFlow does not provide an L-BFGS
optimizer, and so I had to use L-BFGS-B from another software package, SciPy
[6]. An advantage of L-BFGS algorithms is that they have few or no
hyperparameters. Choosing hyperparameters for Adam and SGD required running
many trial inversions. In this example, even when these trial inversions are
included, Adam converges more quickly than L-BFGS-B, but this may not always
be the case. A recently proposed modification of L-BFGS allows it to also be
used with minibatches [4], which may make it more competitive with Adam.

#### Development datasets

Excluding part of the data from training for use as a development dataset is
standard practice in deep learning, but is likely to face resistance in the
FWI community. Part of the reason for the opposition is probably due to
seismic datasets generally having fewer data samples than those used in deep
learning applications, relative to the number of model parameters. The example
of selecting the Adam and SGD hyperparameters demonstrates that the practice
can result in fast convergence to a good model, but does not prove that a
similar result could not have been achieved using the entire dataset. The
necessity of a development dataset is more obvious when choosing
regularization-related hyperparameters, as in that situation the risk of
overfitting is clearer.

#### Other applications of deep learning to FWI

The deep learning approach to FWI I describe here is effectively the
conventional FWI algorithm implemented using deep learning software and using
several concepts popular in deep learning. There are other ways of applying
deep learning techniques to seismic inversion that have less resemblance to
traditional FWI. One of these is to use an approach more obviously inspired by
computer vision applications, where the DNN is trained to transform seismic
data into a model. With this approach, inversion only occurs during training.
Once the network has been trained, it can be applied to datasets without
performing wave propagation. Such a technique would be particularly useful for
specific tasks such as salt body reconstruction [9] and obtaining a low-
wavenumber starting model for traditional FWI [2].

##  5 Conclusion

The hypothesis that FWI can be implemented using an RNN with deep learning
software has been verified. Automatic differentiation calculates gradients
similar to those produced by the traditional adjoint state method approach,
but as the implementation of it in deep learning software is typically not
designed for RNNs as deep as those needed for FWI, the memory requirements are
likely to be prohibitive for realistic applications. The hyperparameters that
were found using a separate development dataset appear to have been good
choices for inverting the example dataset. Combined with the Adam optimizer
and the concept of minibatches, these resulted in significantly faster
convergence than when using the L-BFGS-B optimizer with the cost and gradient
calculated using the entire dataset, and even performed better than SGD.

##  6 Acknowledgments

I am grateful to Bram Willemsen for reviewing an earlier version of this
article.

## References

  * [1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] Mauricio Araya-Polo, Joseph Jennings, Amir Adler, and Taylor Dahlke.  Deep-learning tomography.  The Leading Edge, 37(1):58–66, 2018.
  * [3] Esteban Díaz and Antoine Guitton.  Fast full waveform inversion with random shot decimation.  In SEG Technical Program Expanded Abstracts 2011, pages 2804–2808. Society of Exploration Geophysicists, 2011.
  * [4] Gabriel Fabien-Ouellet, Erwan Gloaguen, and Bernard Giroux.  A stochastic L-BFGS approach for full-waveform inversion.  In SEG Technical Program Expanded Abstracts 2017, pages 1622–1627. Society of Exploration Geophysicists, 2017.
  * [5] Michael Fehler and P Joseph Keliher.  SEAM Phase 1: Challenges of subsalt imaging in tertiary basins, with emphasis on deepwater Gulf of Mexico.  Society of Exploration Geophysicists, 2011.
  * [6] Eric Jones, Travis Oliphant, Pearu Peterson, et al.  SciPy: Open source scientific tools for Python, 2001–.  http://www.scipy.org/; Accessed January 2018.
  * [7] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [8] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521(7553):436–444, 2015.
  * [9] Winston Lewis and Denes Vigh.  Deep learning prior models from seismic images for full-waveform inversion.  In SEG Technical Program Expanded Abstracts 2017, pages 1512–1517. Society of Exploration Geophysicists, 2017.
  * [10] René-Edouard Plessix.  A review of the adjoint-state method for computing the gradient of a functional with geophysical applications.  Geophysical Journal International, 167(2):495–503, 2006.
  * [11] Nitish Srivastava, Geoffrey E Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  Journal of machine learning research, 15(1):1929–1958, 2014.
  * [12] Jean Virieux, Amir Asnaashari, Romain Brossier, Ludovic Métivier, Alessandra Ribodetti, and Wei Zhou.  An introduction to full waveform inversion.  In Encyclopedia of Exploration Geophysics, pages R1–1. Society of Exploration Geophysicists, 2017.
  * [13] Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal.  Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization.  ACM Transactions on Mathematical Software (TOMS), 23(4):550–560, 1997.

##  Appendix A Derivation of automatic differentiation procedure

In this appendix I derive the automatic differentiation procedure for
calculating the gradient of the cost function with respect to the model
parameters. For simplicity, I use the 1D scalar wave equation, ignore the non-
reflecting PML regions, work in discretized space and time coordinates, and
assume that the data consist of only one shot with one receiver.

Using sum of squared errors as the cost function,

$J=\sum_{t=1}^{N_{t}}(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})^{2},$
(2)

where $J$ is the scalar cost function value, $N_{t}$ is the number of time
steps, $\mathbf{u}_{t}$ is a column vector containing the forward propagated
wavefield at time step $t$, $\mathbf{\delta_{x_{r}}}$ is a column vector that
is only non-zero at the location of the receiver (its purpose is to extract
the wavefield amplitude at the receiver), and $d_{t}$ is a scalar containing
the true receiver amplitude at time step $t$.

The cost function depends on the wave speed model implicitly through
$\mathbf{u}_{t}$. When using a second-order accurate central finite difference
scheme to approximate the second-order time derivative in the wave equation,
the wavefield depends on the wave speed, $\mathbf{c}$, and the two previous
time steps, $\mathbf{u}_{t-1}$, and $\mathbf{u}_{t-2}$,

$\mathbf{u}_{t}=\mathbf{c}^{2}\Delta_{t}^{2}\left(\mathbf{D_{x}}^{2}\mathbf{u}_{t-1}-\mathbf{f}_{t-1}\right)+2\mathbf{u}_{t-1}-\mathbf{u}_{t-2},$
(3)

where $\Delta_{t}$ is the time step size, $\mathbf{D_{x}}^{2}$ is a symmetric,
square matrix that applies the second spatial derivative, and
$\mathbf{f}_{t-1}$ is a column vector containing the shot’s source amplitude
at time step $t-1$ as a function of space.

Reverse-mode automatic differentiation uses the chain rule to calculate the
derivative of the cost function value with respect to variables as it proceeds
backward through the computation graph. Because the intermediate variables,
$\mathbf{u}_{t}$, have dependencies between themselves, it is important to
note which variables are being held constant in the partial derivative
notation. The gradient of the cost function with respect to one wavefield time
step is

$\displaystyle\left.\frac{\partial
J}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial
J}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
t}}}+\left.\frac{\partial
J}{\partial\mathbf{u}_{t+1}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{t+1}}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
t,t+1}}}$ (4) $\displaystyle+\left.\frac{\partial
J}{\partial\mathbf{u}_{t+2}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{t+2}}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
t,t+2}}},$

where $\left.\frac{\partial
J}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq t}}}$
indicates the row vector of partial derivatives of $J$ with respect to
elements of $\mathbf{u}_{t}$, while the wave speed, $\mathbf{c}$, and the
wavefields from other time steps, $\mathbf{u_{t^{\prime}\neq t}}$, are held
constant. Variables are only held constant for the partial derivative noted.
The partial derivative vector $\left.\frac{\partial
J}{\partial\mathbf{u}_{t+1}}\right|_{\mathbf{c}}$, for example, includes all
changes to $J$ caused by changes in $\mathbf{u}_{t+1}$, including those due to
dependencies of later wavefield time steps on $\mathbf{u}_{t+1}$.

In a similar way, automatic differentiation calculates the gradient of the
cost function with respect to the wave speed by accumulating terms of the sum
that can be evaluated using the chain rule as the method progresses backward
through time steps,

$\frac{\partial J}{\partial\mathbf{c}}=\sum_{t=1}^{N_{t}}\left.\frac{\partial
J}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{t}}{\partial\mathbf{c}}\right|_{\mathbf{u_{t^{\prime}\neq
t}}}.$ (5)

As an example, if $N_{t}=4$, the steps to calculate the gradient of $J$ would
be

$\left.\frac{\partial
J}{\partial\mathbf{u}_{4}}\right|_{\mathbf{c}}=\left.\frac{\partial
J}{\partial\mathbf{u}_{4}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq 4}}},$
(6) $\left.\frac{\partial
J}{\partial\mathbf{u}_{3}}\right|_{\mathbf{c}}=\left.\frac{\partial
J}{\partial\mathbf{u}_{3}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
3}}}+\left.\frac{\partial
J}{\partial\mathbf{u}_{4}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{4}}{\partial\mathbf{u}_{3}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
3,4}}},$ (7) $\left.\frac{\partial
J}{\partial\mathbf{u}_{2}}\right|_{\mathbf{c}}=\left.\frac{\partial
J}{\partial\mathbf{u}_{2}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
2}}}+\left.\frac{\partial
J}{\partial\mathbf{u}_{3}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{3}}{\partial\mathbf{u}_{2}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
2,3}}}+\left.\frac{\partial
J}{\partial\mathbf{u}_{4}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{4}}{\partial\mathbf{u}_{2}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
2,4}}},$ (8) $\left.\frac{\partial
J}{\partial\mathbf{u}_{1}}\right|_{\mathbf{c}}=\left.\frac{\partial
J}{\partial\mathbf{u}_{1}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
1}}}+\left.\frac{\partial
J}{\partial\mathbf{u}_{2}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{2}}{\partial\mathbf{u}_{1}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
1,2}}}+\left.\frac{\partial
J}{\partial\mathbf{u}_{3}}\right|_{\mathbf{c}}\left.\frac{\partial\mathbf{u}_{3}}{\partial\mathbf{u}_{1}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
1,3}}},$ (9)

where after each step we calculate the term of Equation 5 corresponding to the
current time step.

Calculating the required partial derivatives, we obtain

$\left.\frac{\partial
J}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
t}}}=2(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})\mathbf{\delta_{x_{r}}}^{T},$
(10)
$\left.\frac{\partial\mathbf{u}_{t+1}}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
t,t+1}}}=\mathbf{c}^{2}\Delta_{t}^{2}\mathbf{D_{x}}^{2}+2,$ (11)
$\left.\frac{\partial\mathbf{u}_{t+2}}{\partial\mathbf{u}_{t}}\right|_{\mathbf{c},\mathbf{u_{t^{\prime}\neq
t,t+2}}}=-1,$ (12)
$\left.\frac{\partial\mathbf{u}_{t}}{\partial\mathbf{c}}\right|_{\mathbf{u_{t^{\prime}\neq
t}}}=2\mathbf{c}\Delta_{t}^{2}\left(\mathbf{D_{x}}^{2}\mathbf{u}_{t-1}-\mathbf{f}_{t-1}\right).$
(13)

Substituting into Equation 5, we find

$\displaystyle\frac{\partial J}{\partial\mathbf{c}}$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{N_{t}}\left(2(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})\mathbf{\delta_{x_{r}}}^{T}+\left.\frac{\partial
J}{\partial\mathbf{u}_{t+1}}\right|_{\mathbf{c}}(\mathbf{c}^{2}\Delta_{t}^{2}\mathbf{D_{x}}^{2}+2)-\left.\frac{\partial
J}{\partial\mathbf{u}_{t+2}}\right|_{\mathbf{c}}\right)$ (15)
$\displaystyle\times
2\mathbf{c}\Delta_{t}^{2}\left(\mathbf{D_{x}}^{2}\mathbf{u}_{t-1}-\mathbf{f}_{t-1}\right)$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{N_{t}}\left(\mathbf{c}^{2}\Delta_{t}^{2}\left(\mathbf{D_{x}}^{2}\left.\frac{\partial
J}{\partial\mathbf{u}_{t+1}}\right|_{\mathbf{c}}+\frac{2}{\mathbf{c}^{2}\Delta_{t}^{2}}(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})\mathbf{\delta_{x_{r}}}^{T}\right)\right.$
$\displaystyle+2\left.\left.\frac{\partial
J}{\partial\mathbf{u}_{t+1}}\right|_{\mathbf{c}}-\left.\frac{\partial
J}{\partial\mathbf{u}_{t+2}}\right|_{\mathbf{c}}\right)2\mathbf{c}\Delta_{t}^{2}\left(\mathbf{D_{x}}^{2}\mathbf{u}_{t-1}-\mathbf{f}_{t-1}\right)$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{N_{t}}-F^{-1}\left(\frac{2}{\mathbf{c}^{2}\Delta_{t}^{2}}(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})\mathbf{\delta_{x_{r}}}^{T}\right)2\mathbf{c}\Delta_{t}^{2}\left(\mathbf{D_{x}}^{2}\mathbf{u}_{t-1}-\mathbf{f}_{t-1}\right),$
(16)

where $F^{-1}(\mathbf{g}_{t})$ is wave propagation backward in time of the
source amplitude $\mathbf{g}_{t}$. We also recognize that the wave equation
can be used to express the factor on the right as the second time derivative
of the forward propagated source wavefield (see Equation 3),

$\displaystyle\frac{\partial J}{\partial\mathbf{c}}$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{N_{t}}-F^{-1}\left(\frac{2}{\mathbf{c}^{2}\Delta_{t}^{2}}(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})\mathbf{\delta_{x_{r}}}^{T}\right)\frac{2\Delta_{t}^{2}}{\mathbf{c}}\frac{\partial^{2}\mathbf{u}_{t-1}}{\partial
t^{2}}$ (17) $\displaystyle=$
$\displaystyle\sum_{t=1}^{N_{t}}-F^{-1}\left(2(\mathbf{\delta_{x_{r}}}^{T}\mathbf{u}_{t}-d_{t})\mathbf{\delta_{x_{r}}}^{T}\right)\frac{2}{\mathbf{c}^{3}}\frac{\partial^{2}\mathbf{u}_{t-1}}{\partial
t^{2}},$ (18)

where I used the linearity of wave propagation to take factors out of the
source amplitude of the left term. The result is the same equation that is
used in the adjoint state method (Equation 1).