# IntraSeismic: a coordinate-based learning approach to seismic inversion

Juan Romero
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Wolfgang Heidrich
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Nick Luiken
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Matteo Ravasi
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>

###### Abstract

Seismic imaging is the numerical process of creating a volumetric
representation of the subsurface geological structures from elastic waves
recorded at the surface of the Earth. As such, it is widely utilized in the
energy and construction sectors for applications ranging from oil and gas
prospection, to geothermal production and carbon capture and storage
monitoring, to geotechnical assessment of infrastructures. Extracting
quantitative information from seismic recordings, such as an acoustic
impedance model, is however a highly ill-posed inverse problem, due to the
band-limited and noisy nature of the data. This paper introduces IntraSeismic,
a novel hybrid seismic inversion method that seamlessly combines coordinate-
based learning with the physics of the post-stack modeling operator. Key
features of IntraSeismic are i) unparalleled performance in 2D and 3D post-
stack seismic inversion, ii) rapid convergence rates, iii) ability to
seamlessly include hard constraints (i.e., well data) and perform uncertainty
quantification, and iv) potential data compression and fast randomized access
to portions of the inverted model. Synthetic and field data applications of
IntraSeismic are presented to validate the effectiveness of the proposed
method.

Figure 1: Schematic representation of the IntraSeismic framework: The input to
IntraSeismic consists of spatial coordinates that are processed through a
nonlinear mapping module. Initially, this module maps the input into a higher-
dimensional space of trainable feature vectors using multiresolution hash
encoding. Subsequently, these vectors are input to a MLP that outputs the
impedance model. In the subsequent physics modeling module, the output from
the MLP is added to the background model. A modeling operator then computes
the predicted data from the impedance model, which is then incorporated into
the loss function.

##  1 Introduction

Reflection seismology is the dominant geophysical method used to study the
subsurface [1], making it an indispensable tool for industries seeking to
produce detailed insights into the Earth’s interior. Seismic imaging, similar
to non-light-of-sight imaging [2], relies on the principle of emitting waves
from a source and interpreting the back-scattered signals to gain information
about objects that are not directly visible. A variety of seismic imaging
techniques have been developed over the years and found applications in
several industrial sectors to aid oil and gas prospecting [3, 4, 5],
geotechnical assessment in civil engineering [6], construction projects [7],
and more recently, geothermal production and carbon capture and storage
monitoring [8]. Seismic inversion represents the foremost method for deriving
quantitative information from seismic data [9]. A seismic inverse problem
endeavors to approximate the elastic parameters of the subsurface, such as
wave propagation velocities and density, from either raw or processed seismic
data. Even when this process is carried out post-migration (i.e., after the
raw dataset has been imaged), the associated inverse problem, called post-
stack seismic inversion, is highly ill-posed due to the limited range of
frequencies excited in a seismic experiment and the inevitable presence of
noise in the data.

Historically, this problem has been tackled through three distinct families of
methods: 1. variational methods, 2. deep learning-based approaches, and  3.
hybrid methods. Variational methods aim to minimize a data-misfit term
accompanied by one or more regularization terms that guide the solution to a
physically plausible outcome. A standard regularization choice is Tikhonov
regularization [10], which favors the retrieval of a smooth subsurface model.
More recently, total-variation regularization [11] has been used to enforce a
layered structure of the subsurface, which is a more realistic representation
of the acoustic impedance [12, 13]. Deep learning-based approaches, on the
other hand, leverage the prowess of neural networks to learn a direct mapping
between the seismic data and acoustic impedance model, providing that a
suitable training dataset is available. Various architectures have been used
so far in the literature, including convolutional neural networks [14, 15],
temporal convolutional networks [16], and attention-based networks [17].
However, seismic applications deal with a shortage of labeled datasets because
the ground truth is inaccessible except for specific locations along drilled
wells (see Section 3 for further details). Consequently, the application of
supervised learning approaches is limited in this domain and is feasible only
when a dense well coverage is available [18]. Hybrid seismic inversion
methodologies harness the strengths of both paradigms by amalgamating the
foundational geological knowledge and physics-based modeling with the pattern
recognition capabilities of deep learning techniques. Examples of hybrid
methods are Deep Image Prior (DIP) [19], which uses a neural network as a
nonlinear preconditioner in the solution of an inverse problem and has been
successfully applied to the post-stack seismic inversion problem [20].
Similarly, the Plug-and-Play (PnP) framework [21] has recently been proposed
to solve inverse problems using a denoiser to act as regularizer; this concept
has led to the seamless inclusion of pre-trained deep denoisers [22] in the
solution of inverse problems of any kind: [23] showed that this is also
beneficial to 2D post-stack seismic inversion, achieving state-of-the-art
results. For an extensive comparison of some of the previously described
methods, we refer the reader to [18, 24].

These hybrid methods represent a forward-thinking direction for seismic
inversion applications. Model-driven components ensure inversion results to be
grounded in well-understood geophysical principles, providing interpretable
and physically plausible solutions. Simultaneously, deep-learning-based
components can capture complex, nonlinear relationships in the data and/or
model that might be overlooked or overly simplified in purely model-driven
techniques. In this spirit, we introduce IntraSeismic, a hybrid method that is
able to:

  * •

Leverage the power of coordinate-based learning to solve ill-posed seismic
inverse problems whilst retaining the physics of the seismic modeling operator
through a modular framework.

  * •

Seamless integrate hard constraints in the seismic inversion process.

  * •

Perform uncertainty quantification via Monte Carlo Dropout.

  * •

Provide a compact representation of subsurface models, serving as an efficient
data compression mechanism.

As shown by the numerical examples, IntraSeismic surpasses standard
variational methods and other hybrid approaches in the inversion of both 2D
and 3D post-stack seismic data in terms of reconstruction quality and
convergence. This paper begins with an overview of post-stack seismic
inversion, coordinate-based learning, and its application to inverse problems.
Subsequently, we detail the components of the IntraSeismic framework and
demonstrate the effectiveness of the proposed method on two synthetic
datasets, namely the 2D Marmousi model [25] and a subvolume of the 3D SEAM
Phase 1 model [26], and the Sleipner field dataset [27]. Finally, we
illustrate the ability of IntraSeismic to include hard constraints, perform
uncertainty quantification, and act as a flexible data compressor.

##  2 Background and related work

#### Post-stack seismic inversion

Post-migration seismic data (or seismic images) can be mathematically
represented via the so-called convolutional model [28]. This entails the
convolution of a source function or wavelet $w(t)$ with the earth reflectivity
series $r(t)$:

$d(t)=w(t)*r(t),$ (1)

where $t$ indicates the time axis, $r(t)=(I_{p}(t+\Delta
t)-I_{p}(t))/(I_{p}(t+\Delta t)+I_{p}(t))$, and $I_{p}$ corresponds to the
acoustic or $\mathrm{P}$-wave impedance. As shown in [29], a logarithmic
approximation holds for small to medium contrasts between geological layers
and renders Equation 1 linear in $I_{p}(t)$:

$r(t)=\frac{1}{2}\frac{d}{dt}\log(I_{p}(t)).$ (2)

Hence, post-stack seismic data can be modeled as

$d(t)=\frac{1}{2}w(t)*\frac{d}{dt}\log(I_{p}(t)),$ (3)

which can be written in compact matrix-vector notation:

$\mathbf{d}=\mathbf{WDm},$ (4)

where $\mathbf{W}$ represents a block-Toeplitz convolution matrix, which
encapsulates the seismic wavelet, while $\mathbf{D}$ is the first-order
derivative operator. Post-stack seismic inversion transforms seismic images
$\mathbf{d}\in\mathbb{R}^{N_{y}N_{x}N_{t}}$ into quantitative estimates of the
subsurface’s acoustic impedance $\mathbf{m}\in\mathbb{R}^{N_{y}N_{x}N_{t}}$
[30, 31], and it is expressed as:

$\mathbf{m^{*}}=\operatorname*{arg\,min}_{\mathbf{m}}\frac{1}{2}\left\|\mathbf{Gm}-\mathbf{d}\right\|_{2}^{2}+\mathcal{R}(\mathbf{m}),$
(5)

where $\mathbf{G}=\mathbf{WD}$ and $\mathcal{R}$ corresponds to the chosen
regularization term.

Coordinate-based learning for inverse problems Coordinate-based learning,
wherein neural networks directly utilize spatial coordinates as input, has
emerged as a focal point in recent deep learning advancements [32]. This
approach has achieved state-of-the-art results in various domains, including
image and shape representation [33, 34] and volume density and view synthesis
[35], to name a few. Neural Radiance Fields (NeRFs) [35] is an example of
pioneering work in this domain, which has showcased the potent application of
spatial coordinates in rendering photo-realistic 3D scenes.

A recent notable contribution to the inverse problems field is the IntraTomo
framework [36]. IntraTomo synergizes a coordinate-based multi-layer perceptron
(MLP) to estimate a density field with a subsequent geometry refinement module
that solves an optimization problem to refine the details in such a density
field. The data- and model-driven combination of IntraTomo has achieved state-
of-the-art results in computer tomography. In the domain of seismic inversion,
a parallel methodology has been suggested, involving the inversion of full-
waveform seismic data using spatial coordinates as inputs to an MLP with
variable activation functions [37].

In the present study, we draw inspiration from the IntraTomo and further
develop this concept to produce a new technique called "IntraSeismic". In
IntraSeimic, we take advantage of multi-resolution hash encoding with
trainable feature vectors as a way to enrich the expressiveness of the
network’s input [38]. As such, we are able to employ a very small, streamlined
MLP architecture with two layers of 64 neurons each. Finally, we assimilate
the ’geometry refinement module’ of IntraTomo, by directly embedding model-
based regularization terms in the neural network’s loss function.
Consequently, IntraSeismic emerges as a compact and computationally efficient
framework that is tailor-made for solving highly ill-posed inverse problems in
geophysics.

##  3 The IntraSeismic framework

IntraSeismic integrates a self-supervised, coordinate-based learning module
with a modeling operator to invert seismic data for acoustic impedance
subsurface models (Figure 1). The proposed framework primarily uses
coordinates as input, which first undergo a multi-resolution hash encoding
process [38]. Next, the encoded coordinates are passed through an MLP, whose
output is subsequently fed to the seismic modeling operator to allow
evaluation of the loss function:

$\displaystyle\operatorname*{arg\,min}_{\Theta}\frac{1}{2}\left\|\mathbf{G}\boldsymbol{m}_{\Theta}-\mathbf{d}\right\|_{2}^{2}+\mathcal{R}(\boldsymbol{m}_{\Theta})$
(6)
$\displaystyle\boldsymbol{m}_{\Theta}=\mathbf{F}_{\Theta}(\boldsymbol{x})+\boldsymbol{m}_{\text{back}}\hskip
28.45274pt,$

where $\boldsymbol{m}_{\Theta}$ denotes the subsurface acoustic impedance
volume obtained from the inversion process.
$\mathbf{F}_{\Theta}(\boldsymbol{x})$ is the network with $\Theta$
representing the learnable parameters from both the hash table and the MLP
network, and $\mathbf{x}\in\mathbb{R}^{N_{y},N_{x},N_{t}\times 3}$ corresponds
to the spatial coordinates of the volume of interest. These coordinates can be
processed in IntraSeismic either as a single batch or in mini-batches. The
term $\boldsymbol{m}_{\text{back}}$ refers to the background model. The
regularization term, $\mathcal{R}(\boldsymbol{m}_{\Theta})$, employs
anisotropic TV regularization on the model
$\left(\mathcal{R}(\boldsymbol{m}_{\Theta})=\left\|\boldsymbol{m}_{\Theta}\right\|_{TV}\right)$
to achieve a high-resolution layered representation of the subsurface. Our
experiments demonstrate improved performance when incorporating additional
regularization on the network’s output, such as sparsity-promoting
regularization
$\mathcal{R}(\boldsymbol{m}_{\Theta})=\left\|\boldsymbol{m}_{\Theta}\right\|_{1}$,
or dynamically scheduled $L_{2}$ regularization
$\mathcal{R}(\boldsymbol{m}_{\Theta})=\beta_{k}\left\|\boldsymbol{m}_{\Theta}\right\|_{2}^{2}$,
where the regularization weight $\beta$ diminishes across iterations $k$ (for
further details refer to the supplementary material). In the following text,
we elucidate the different components of IntraSeismic.

Input coordinates: A post-stack seismic dataset is a 3D volume ($V$) of
spatial size $N_{y},N_{x},N_{t}$ that can be represented via voxel
discretization. Each voxel $p_{i}$ within the volume $V$ is uniquely defined
by its spatial coordinates $\boldsymbol{x_{i}}=(y,x,t)$. As shown in equation
6, the collection of all coordinates of the volume of interest represents the
input of the IntraSeismic framework. Note the third dimension $N_{t}$ is time
samples, as seismic data is recorded in time.

Background Model: A background or initial prior model plays a pivotal role in
seismic inverse problems. This model must contain the low-frequency components
of the subsurface impedance model, which are not captured in the data. In
variational methods, the purpose of the background model is to steer the
inverse process towards geologically consistent solutions. When using
synthetic data, the background model is typically a smoothed version of the
ground truth. In contrast, for field data, the background model is derived
from an interpolated, smoothed version of well-log data. Well-log data refer
to a suite of physical measurements acquired inside wells using various
logging tools. These logs provide a detailed record of the geological
formations the well has penetrated, offering crucial information such as the
rock’s seismic velocity, density, resistivity, and natural radiation, among
others. Well-log data and seismic data differ significantly in scale and
resolution: well-logs are high-frequency records with data points captured at
intervals as small as a few centimeters, whereas seismic data have lower
frequencies, typically ranging from 10 Hz to 100 Hz, corresponding to a
resolution measured in meters.

Multiresolution hash encoding: Encoding input coordinates into a higher-
dimensional space via pre-defined Fourier feature embeddings (or Gaussian-
based alternatives) has been shown to significantly increase the ability of
coordinate-based learning to capture high-frequency information in the sought-
after solution. Results are reported across a variety of computer vision tasks
[32, 35, 39, 40]. Notably, [38] introduced multiresolution hash encoding, a
sparse grid-based parametric embedding that outperformed state-of-the-art non-
parametric frequency-based encodings. This method utilizes a multi-resolution
hash table with trainable feature vectors that are optimized through
backpropagation, alongside the subsequent neural network parameters. This
method leverages spatial hashing, where coordinates serve as the keys for a
hash function, facilitating the retrieval of hash table indices. The hash
function $h$ for a coordinate point $\mathbf{x}$ is defined as:

$h(\mathbf{x})=\left(\bigoplus_{i=1}^{3}x_{i}\pi_{i}\right)\bmod T,$ (7)

where $\bigoplus$ denotes the binary XOR operation, and $\pi_{i}$ represents a
large prime number chosen by the user. This function computes the indices for
a hash table of size $T$ that contains feature vectors of dimension $F$ at
every entry. Given its multi-resolution nature, this operation is performed at
multiple predefined resolutions. Each input coordinate is scaled to a specific
resolution, and through floor and ceiling operations, the surrounding points’
corners at a given resolution are ascertained—yielding four corner points for
a 2D problem and eight for a 3D problem. These corners are fed into the hash
function, with the resultant feature vectors being interpolated to obtain the
corresponding feature vector at the input point.

Neural Network Architecture: The neural network component of IntraSeismic is a
two-layer MLP network with ReLU activation function. Each layer contains 64
neurons.

Modeling Operator: In the applications presented in subsequent sections, we
utilize the post-stack seismic modeling operator, as defined in Equation 4.

Optimizer: We use the Adam optimizer [41] with a learning rate of $10^{-3}$ to
solve for the IntraSeismic parameters from the MLP network and hash encoding
tables.

Implementation details: The IntraSeismic framework is fully implemented in
PyTorch 2.0.1, and all experiments are conducted on an AMD EPYC 7713 64-Core
Processor equipped with a single NVIDIA TESLA A100. The modeling operator is
implemented and wrapped into PyTorch using the PyLops library [42].

##  4 Experiments and results

###  4.1 Synthetic data

To assess the effectiveness and performance of IntraSeismic across different
scenarios, we apply the proposed algorithm to both 2D and 3D synthetic post-
stack seismic datasets, incorporating band-pass Gaussian noise at varying
standard deviations: noise-free ($\sigma=0$), $\sigma=0.1$, and $\sigma=0.2$.
It is worth noting that the seismic dataset is scaled from $-1$ to $1$. In
Table 1, we report both the Signal-to-Noise Ratio (SNR) of the obtained
solution and the number of iterations needed to achieve such a result. In the
subsequent sections, we illustrate the results when using a moderate noise
level ($\sigma=0.1$), which is reflective of the usual noise levels
encountered in field data. Figures 2a and 4a display the corresponding data
for the Marmousi and SEAM models. In both cases, the seismic data is modeled
using a Ricker wavelet [43] with a peak frequency of 15 Hz, and the background
models are smoothed versions of the original acoustic impedance models
(Figures 2b, and 4b).

Table 1: Inversion results for IntraSeismic and benchmark methods for the
Marmousi and SEAM datasets and different noise levels. T SNR /Iterations are
shown in each entry, and bold is used to indicate the best performing
algorithm for any combination of model and noise.

.  Marmousi SEAM Method $\sigma=0$ $\sigma=0.1$ $\sigma=0.2$ $\sigma=0$
$\sigma=0.1$ $\sigma=0.2$ Tikhonov [23] 21.50 /71 21.17 /69 20.35 /69 34.11
/352 33.32 /400 31.53 /332 PD-TV [23] 24.74 /4000 23.92 /38.70 22.37 /3580
36.9 /4000 36.01 /3000 34.50 /2000 DIP [18] 22.19 /1475 20.21 /590 17.61 /306
36.1 /7000 33.4 /230 28.00 /7000 DIP enhanced 24.64 /1238 24.00 /407 23.71
/349 38.44 /4100 34.6 /4600 32.66 /2600 PnP [23] 25.25 /1280 24.54 /1660 23.37
/1230 N/A N/A N/A IntraSeismic 27.34 /1500 25.00 /457 22.93 /437 40.84 /393
36.92 /217 34.5 /307

The IntraSeismic result on synthetic data is benchmarked with two variational
approaches, Tikhonov regularization (the standard method used in industry) and
total-variation regularization (TV) using the Primal-Dual solver [44].
Additionally, we compare IntraSeismic with the state-of-the-art hybrid methods
proposed in the literature for post-stack seismic inversion, namely Deep Image
Prior [20] and Plug-and-Play priors with deep denoisers [23].

For the PnP results, we utilized the published code [23]. This code employs
the DRUNet denoiser [22] within a Primal-Dual framework, functioning as a
regularizer which achieved state-of-the-art outcomes in 2D post-stack seismic
inversion. As of the publication date of our paper, we are not aware of any
published implementation of PnP for 3D data. Consequently, we do not include a
comparison with PnP in our 3D synthetic and field data experiments. In the
case of DIP, [18] showed that the originally proposed method becomes sub-
optimal in the presence of noise in post-stack seismic inversion. This
limitation may stem from the absence of prior information that guides the
inversion toward geologically consistent solutions. In our study, we further
modify DIP by incorporating the same regularization terms as those used in
IntraSeismic (see equation 6), and we refer to it as DIP enhanced. We believe
that this facilitates a more equitable comparison between the DIP and
IntraSeismic methods, where their main difference lies in the parameterization
of the acoustic impedance model (i.e., a CNN that warps a noise vector into
the sought-after solution in DIP, and an MLP that transforms coordinates into
values of the solution in IntraSeismic). Finally, we decided not to provide
any comparison with state-of-the-art supervised learning approaches as that
would introduce additional factors (e.g., size of the training dataset) that
greatly affect the quality of the final solution.

Marmousi model
The Marmousi model [25] is a synthetically created, realistic 2D subsurface
velocity model that emulates complex geological structures with sharp velocity
contrasts (see Figure 2i). The Marmousi model is the most widely used model in
geophysics, thereby serving as a benchmark for various seismic imaging
algorithms. Figures 2(c to f) illustrate the inversion results for the
Marmousi dataset at a noise level of $\sigma=0.1$, comparing standard
variational methods, hybrid algorithms, and IntraSeismic. IntraSeismic
retrieves a high-resolution acoustic impedance image of the Marmousi model,
where discontinuities and layer contrasts are clearly defined. As indicated in
Table 1, IntraSeismic achieves the highest SNR in the Marmousi model for both
noise-free conditions and noisy data with $\sigma=0.1$, and it is slightly
surpassed by PnP for the highest noise level ($\sigma=0.2$). Moreover, Table 1
and Figure 3 highlight the fast convergence of IntraSeismic, reaching the
highest SNR in less than $500$ iterations for noisy data.

Figure 2: Inversion of the synthetic 2D Marmousi dataset with a noise level of
$\sigma=0.1$: a) seismic data, b) background model, c) Tikhonov
regularization, d) TV-regularization using the Primal-Dual solver, e) DIP, f)
DIP with additional priors, g) PnP, h) IntraSeismic, and i) ground truth.
Figure 3: SNR evolution through iterations for IntraSeismic and benchmark
methods in the inversion of the Marmousi data.

SEAM model
The model employed in this study is a volume of size $180\times 180\times
180$, extracted from the SEAM Phase I velocity model created by the SEG
Advanced Modeling Consortium [26]. The SEAM model showcases realistic
geological heterogeneities and discontinuities across all three spatial
dimensions and exhibits high contrasts between some of its layers (Figure 4g).
Evidently, the IntraSeismic method (Figure 4f) produces the solution with the
highest SNR, leading to a more precise 3D reconstruction of the SEAM model. As
demonstrated in Table 1, IntraSeismic achieves the highest SNR across all
noise levels for the SEAM dataset. Moreover, it converges more rapidly than
existing benchmark methods, requiring fewer than 400 iterations in every
evaluated noise scenario. This is illustrated in Figure 5, where the SNR
evolution over iterations is shown for some of the benchmark methods and
IntraSeismic for the case of noisy data with $\sigma=0.1$.

Figure 4: Inversion of the synthetic 3D SEAM dataset: a) seismic data, b)
background model, c) Tikhonov regularization, d) TV-regularization using the
Primal-Dual solver, e) DIP, f) DIP with additional priors, g) IntraSeismic and
h) ground truth. Figure 5: SNR evolution through iterations for IntraSeismic
and benchmark methods in the inversion of the SEAM data.

###  4.2 Field data: Sleipner dataset

The Sleipner seismic dataset was acquired in 2001 in the North Sea as part of
a seismic monitoring program of an ongoing carbon capture and storage project
in the Sleipner field [27]. We use a subvolume of size $300\times 300\times
250$ (Figure 6a), which exhibits a predominantly horizontal geology disrupted
by strong contrast values in the CO2 reservoir. The statistical wavelet and
the background model (Figure 6b) are taken from [45]. Figure 6 shows the
inversion results for the benchmark methods and IntraSeismic. In this case, we
do not observe any major scale difference between the IntraSeismic and PD-TV
results; however, when looking at the inversion results in detail,
IntraSeismic shows smoother transitions in the horizontal direction (a sought-
after characteristic in geology) and heightened resolution in the vertical
direction that accentuates the layer contrasts. After 200 iterations (spanning
40 seconds), we notice minimal changes in the model update, which we interpret
as the convergence point.

Figure 6: Inversion of Sleipner dataset using IntraSeismic and benchmark
methods: a) seismic data, b) background model, c) Tikhonov regularization, d)
TV-regularization using the Primal-Dual solver, e) DIP, f) DIP with additional
priors, and g) IntraSeismic.

###  4.3 Uncertainty quantification

Uncertainty quantification (UQ) is pivotal in any seismic inversion process
given the non-uniqueness of the seismic inverse problem, the inherent
complexity and variability of the Earth’s subsurface, and the previously
discussed limitations of the observed data. As such, UQ offers insight into
the range of possible subsurface models that can equally fit the observed
data. This knowledge is instrumental for informed decision-making in scenarios
with significant financial and environmental implications [46]. Furthermore,
understanding the uncertainties associated with a seismic inversion product
can guide subsequent data acquisition campaigns to refine results and provide
a clearer picture of result reliability and associated risks. When a neural
network is involved directly (via supervised learning) or indirectly (as a
regularizer or preconditioner) in the solution of the seismic inverse problem,
additional uncertainty is introduced in the process. Among other techniques to
quantify such an uncertainty component, Monte Carlo Dropout (MCD) has been
shown to be an effective and easy-to-implement method to quantify the so-
called epistemic uncertainty [47]. By retaining the dropout layers during the
inference phase, rather than just during training, and performing multiple
evaluations of the network with different neurons "dropped out" each time, one
can obtain samples of an approximate posterior distribution [48]. We propose
to add dropout layers to IntraSeismic’s MLP network and to perform MCD
sampling once the optimization process is completed. The proposed technique is
applied to the Marmousi model. More specifically, in this case, we apply
dropout to both of the layers of the MLP network (Figure 1) with a dropout
probability of $0.2$, and we compute a total of $100$ realizations, which are
used to evaluate sample statistics (e.g., sample mean and standard deviation).
Figure 7 corresponds to the obtained standard deviation divided by the mean of
the Marmousi model parameters that shows the IntraSeismic’s confidence levels
and highlights areas in the model with potential ambiguity in the solution.
The layers with higher impedance values exhibit the highest uncertainty,
possibly due to TV regularization, which might lead to a diminished contrast.

Figure 7: Monte Carlo Dropout uncertainty quantification result for Marmousi
model: a) std/mean, the red line indicates the extracted column shown in b).

###  4.4 Adding hard-constraints

Well-log data provide the closest approximation to the ground truth in seismic
inversion, as they are actual recordings of subsurface physical properties. As
such, they can serve as critical anchor points, ensuring that the seismic
inversion results are consistent with real-world measurements at specific
locations. One way to incorporate well-log data in the seismic inversion
process is by integrating it directly into the original optimization
framework. This can be achieved by introducing a soft constraint to the
optimization problem as formulated in equation 5. Specifically, the term
$\gamma\left\|\boldsymbol{m}_{i}-\boldsymbol{m}_{\text{well}}\right\|_{2}^{2}$
can be added, where $\gamma$ is a weighting factor that balances the influence
of the well-log data, $\mathbf{m_{i}}$ represents the model parameters at the
well position, and $\boldsymbol{m}_{\text{well}}$ denotes the measurements
from the well-log. This addition ensures that the inversion solution aligns
more closely with the well-log observations. An alternative for incorporating
well-log data into the seismic inversion process involves the application of a
hard constraint, which enforces an exact match between the impedance model and
the well-log data ($\boldsymbol{m}_{i}=\boldsymbol{m}_{\text{well}}$). This
approach is particularly advantageous in instances of highly noisy data where
the integrity of well data is prioritized or in scenarios with numerous wells
within a development field, where the objective is to construct a robust model
that accurately reflects all available well data. In IntraSeismic, the
integration of well data can be naturally achieved by redefining the network’s
output as follows [49]:

$\boldsymbol{m}=(\boldsymbol{x}-\boldsymbol{x}_{\text{well}})\Delta\boldsymbol{m}_{\Theta}+\boldsymbol{m}_{\text{well}}$
(8)

where $\boldsymbol{x}$ represents the coordinate vector, while
$\boldsymbol{x}_{\text{well}}$ and $\boldsymbol{m}_{\text{well}}$ denote the
coordinates and values of the well-log data, respectively. The term
$\Delta\boldsymbol{m}_{\Theta}$ corresponds to the output of the neural
network. Within this framework, the neural network’s output ensures compliance
with the hard constraint and acts as an update to the well data away from its
position. This approach is implemented for the Sleipner dataset, forcing the
solution to match the available acoustic impedance well-log in the middle of
the model (Figure 8).

Figure 8: Hard-constrained inversion of Sleipner dataset: a) Output of the
network $\Delta\boldsymbol{m}_{\Theta}$, b) final estimated acoustic impedance
model and c) comparison at the well location of the inversion with and without
hard-contraint.

###  4.5 Irregular data access and data compression

Geoscientists are constantly tasked to analyze vast amounts of seismic data, a
process usually referred to as seismic interpretation [50]. As such, the
property models produced by different seismic inversion algorithms are
commonly visualized along different axes (i.e., inline, crossline, time/depth)
or along irregular surfaces corresponding to key geological horizons.
Conventional seismic interpretation software packages perform such extraction
procedures operating directly with data stored on fast-access file storage
systems or on local memory when the data size allows it.
However, a user must select upfront the layout of the data, placing the most
commonly accessed direction over the fast axis and the least commonly accessed
direction over the slow one or replicating the same information over multiple
files with different layouts. In order to reduce the amount of I/O, the
seismic industry has recently placed tremendous focus on compression
algorithms, with the popular ZFP data format [51] being widely adopted in both
proprietary and open-source data formats (e.g., MDIO [52]). A second by-
product of the IntraSeismic framework is represented by the coordinate-based
parameterization of the sought-after acoustic impedance model. As such, the
cost of accessing the model values at $N$ (possibly random) coordinates is
independent of the location and order of such points. We believe that this
unique feature of IntraSeismic could revolutionize how seismic interpretation
software packages store and access seismic volumes.
Moreover, given that the complete knowledge of the acoustic impedance model is
encapsulated within the learned hash tables and the parameters of the MLP, the
IntraSeismic framework can be considered an effective data compression tool.
The SEAM model, illustrated in Figure 4, comprises of $5.83$ million data
points, yet the IntraSeismic framework can represent it with only $1.16$
million parameters. This yields a compression ratio of approximately $5$. In
contrast, the ZFP compression algorithm achieves a compression ratio of
approximately $15$ for the SEAM model, with an error margin of $0.1\%$. A
significant factor in the high compression efficiency of ZFP can be attributed
to its quantization step performed prior to storing the compressed data. With
quantization emerging as a focal point of interest within the field of Deep
Learning —for its utility in representing neural network weights in a reduced
and efficient format— it is plausible to anticipate an improvement in
IntraSeismic’s compression ratio with the application of quantization to the
hash table and MLP parameters.

##  5 Conclusions and future directions

In this study, we have introduced the IntraSeimic framework, a hybrid data-
driven and model-driven strategy to solve seismic inverse problems. We
presented and critically evaluated the IntraSeismic framework by applying it
to both realistic synthetic and field post-stack seismic datasets, and
benchmarked it against prevalent variational approaches and deep-learning-
based methods in the seismic inversion domain. Our findings underscore
IntraSeismic’s superior capability in producing high-resolution impedance
models and faster convergence. Additionally, the IntraSeismic framework’s
inherent coordinate-based parametrization offers revolutionary prospects for
the seismic interpretation domain, especially in terms of data access and
compression. The introduction of Monte Carlo Dropout to IntraSeismic further
empowers it to quantify uncertainties, a critical aspect in seismic inversion.

The visual similarity of the IntraSeismic results to those obtained using the
PD-TV method stems mainly from the shared use of the TV regularization
strategy. Both IntraSeismic and PD-TV exhibit challenges in recovering high
impedance layers—a feat where the PnP method showed superiority in the
Marmousi example. Such limitations might arise from selecting the
regularization term itself, which can inadvertently constrain the framework’s
modeling capacity, leading to inherent trade-offs like the loss of contrast
commonly associated with TV regularization. Therefore, future research should
aim to explore learned regularization strategies (like PnP) within
IntraSeismic-type frameworks. Finally, while this study has focused on the
application of IntraSeismic to post-stack seismic inversion, the framework’s
design allows for straightforward adaptation to other inversion scenarios,
such as pre-stack seismic inversion and image domain least-squares migration
(additional details provided in the supplementary material).

## Acknowledgments

The authors thank King Abdullah University of Science & Technology (KAUST) for
supporting this research as well as Equinor and partners for releasing the 4D
Sleipner dataset (available at
https://co2datashare.org/dataset/sleipner-2019-benchmark-model). For computer
time, this research used the resources of the Supercomputing Laboratory at
KAUST in Thuwal, Saudi Arabia.

## References

  * [1] Philip Kearey, Michael Brooks, and Ian Hill.  An Introduction to Geophysical Exploration, 3rd Edition | Wiley, 2002.
  * [2] Daniele Faccio, Andreas Velten, and Gordon Wetzstein.  Non-line-of-sight imaging.  Nature Reviews Physics, 2:318––327, 2020.
  * [3] William Murray Telford, Lloyd P. Geldart, and Robert E. Sheriff.  Applied Geophysics, October 1990.  ISBN: 9781139167932 Publisher: Cambridge University Press.
  * [4] Robert E. Sheriff and Lloyd P. Geldart.  Exploration Seismology.  Cambridge University Press, Cambridge, 2 edition, 1995.
  * [5] Öz Yilmaz.  Seismic data analysis.  Society of Exploration Geophysicists, 2001.
  * [6] Öz Yilmaz and Richard D. Miller.  Engineering Seismology: With Applications to Geotechnical Engineering.  Society of Exploration Geophysicists, January 2015.
  * [7] Don W. Steeples and Richard D. Miller.  Seismic Reflection Methods Applied To Engineering, Environmental, And Ground-Water Problems.  page cp. European Association of Geoscientists & Engineers, March 1988.  ISSN: 2214-4609.
  * [8] David Lumley.  4D seismic monitoring of CO2 sequestration.  The Leading Edge, 29(2):150–155, February 2010.  Publisher: Society of Exploration Geophysicists.
  * [9] Albert Tarantola.  Inversion of seismic reflection data in the acoustic approximation.  GEOPHYSICS, 49(8):1259–1266, August 1984.  Publisher: Society of Exploration Geophysicists.
  * [10] Andrey N. Tikhonov and Vasiliy .Y. Arsenin.  Solution of Ill-posed Problems.  Washington: Winston & Sons, 1977.
  * [11] Leonid I. Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena, 60(1):259–268, November 1992.
  * [12] Ali Gholami.  Nonlinear multichannel impedance inversion by total-variation regularization.  Geophysics, 80:R217–R224, 2015.
  * [13] Dehua Wang, Jinghuai Gao, and Hongan Zhou.  Data-driven multichannel seismic impedance inversion with anisotropic total variation regularization.  Journal of Inverse and Ill-Posed Problems, 2018(2):229–241, 2019.
  * [14] Reetam Biswas, Mrinal K. Sen, Vishal Das, and Tapan Mukerji.  Prestack and poststack inversion using a physics-guided convolutional neural network.  Interpretation, 7(3), August 2019.  Publisher: Society of Exploration Geophysicists.
  * [15] Vishal Das, Ahinoam Pollack, Uri Wollner, and Tapan Mukerji.  Convolutional neural network for seismic impedance inversion.  GEOPHYSICS, 84(6):R869–R880, November 2019.  Publisher: Society of Exploration Geophysicists.
  * [16] Robert Smith, Philippe Nivlet, Hussain Alfayez, and Nasher AlBinHassan.  Robust deep learning-based seismic inversion workflow using temporal convolutional networks.  Interpretation, 10(2):SC41–SC55, May 2022.  Publisher: Society of Exploration Geophysicists.
  * [17] Bangyu Wu, Qiao Xie, and Baohai Wu.  Seismic Impedance Inversion Based on Residual Attention Network.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–17, 2022.  Conference Name: IEEE Transactions on Geoscience and Remote Sensing.
  * [18] Matteo Ravasi, Nick Luiken, Juan Romero, and Miguel Corrales.  Deep Learning to replace or augment model-based seismic inversion?  volume 2023, pages 1–5. European Association of Geoscientists & Engineers, June 2023.  ISSN: 2214-4609 Issue: 1.
  * [19] Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Deep Image Prior.  International Journal of Computer Vision, 128(7):1867–1888, July 2020.  arXiv:1711.10925 [cs, stat].
  * [20] Vicenzo Lipari, Francesco Picetti, P. Bestagini, and Stefano Tubaro.  Post-Stack Inversion with Uncertainty Estimation through Bayesian Deep Image Prior.  volume 2021, pages 1–5. European Association of Geoscientists & Engineers, October 2021.  ISSN: 2214-4609 Issue: 1.
  * [21] Singanallur V. Venkatakrishnan, Charles A. Bouman, and Brendt Wohlberg.  Plug-and-Play priors for model based reconstruction.  In 2013 IEEE Global Conference on Signal and Information Processing, pages 945–948, December 2013.
  * [22] Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte.  Plug-and-Play Image Restoration With Deep Denoiser Prior.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(10):6360–6376, October 2022.  Conference Name: IEEE Transactions on Pattern Analysis and Machine Intelligence.
  * [23] Juan Romero, Miguel Corrales, Nick Luiken, and Matteo Ravasi.  Plug and Play Post-Stack Seismic Inversion with CNN-Based Denoisers.  volume 2022, pages 1–5. European Association of Geoscientists & Engineers, October 2022.  ISSN: 2214-4609 Issue: 1.
  * [24] Matteo Ravasi, Juan Romero, Miguel Corrales, Nick Luiken, and Claire Birnie.  Chapter Six - Striking a balance: Seismic inversion with model- and data-driven priors.  In Shib Sankar Ganguli and Vijay Prasad Dimri, editors, Developments in Structural Geology and Tectonics, volume 6 of Reservoir Characterization, Modeling, and Quantitative Interpretation, pages 153–200. Elsevier, January 2023.
  * [25] Aline Brougois, Marielle Bourget, Patrick Lailly, Michel Poulet, Patrice Ricarte, and Roelof Versteeg.  Marmousi, model and data.  page cp. European Association of Geoscientists & Engineers, May 1990.  ISSN: 2214-4609.
  * [26] 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.
  * [27] Anne-Kari Furre, Ola Eiken, Håvard Alnes, Jonas Nesland Vevatne, and Anders Fredrik Kiær.  20 years of monitoring co2-injection at sleipner.  Energy Procedia, 114:3916–3926, 2017.  13th International Conference on Greenhouse Gas Control Technologies, GHGT-13, 14-18 November 2016, Lausanne, Switzerland.
  * [28] Pierre L. Goupillaud.  An approach to inverse filtering of near-surface layer effects from seismic records.  GEOPHYSICS, 26(6):754–760, December 1961.  Publisher: Society of Exploration Geophysicists.
  * [29] R.H. Stolt and A.B. Weglein.  Migration and inversion of seismic data.  Geophysics, 50:2458–2472, 1985.
  * [30] D.W. Oldebnbur, T. Scheuer, and S. Levy.  Recovery of the acoustic impedance from reflection seismograms.  Geophysics, 48:1318–1337, 1983.
  * [31] Brian Russell and Dan Hampson.  Comparison of poststack seismic inversion methods.  61st Annual International Meeting, SEG, pages 876–878, 1991.
  * [32] Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng.  Fourier Features Let Networks Learn High Frequency Functions in Low Dimensional Domains, June 2020.  arXiv:2006.10739 [cs].
  * [33] Kenneth O. Stanley.  Compositional pattern producing networks: A novel abstraction of development.  Genetic Programming and Evolvable Machines, 8(2):131–162, June 2007.
  * [34] Zhiqin Chen and Hao Zhang.  Learning Implicit Fields for Generative Shape Modeling, September 2019.  arXiv:1812.02822 [cs].
  * [35] Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng.  NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis, August 2020.  arXiv:2003.08934 [cs].
  * [36] Guangming Zang, Ramzi Idoughi, Rui Li, Peter Wonka, and Wolfgang Heidrich.  IntraTomo: Self-supervised Learning-based Tomography via Sinogram Synthesis and Prediction.  In 2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 1940–1950, October 2021.  ISSN: 2380-7504.
  * [37] Jian Sun, Kristopher Innanen, Tianze Zhang, and Daniel Trad.  Implicit seismic full waveform inversion with deep neural representation.  Journal of Geophysical Research: Solid Earth, 128(3):e2022JB025964, 2023.
  * [38] Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller.  Instant Neural Graphics Primitives with a Multiresolution Hash Encoding.  ACM Transactions on Graphics, 41(4):1–15, July 2022.  arXiv:2201.05989 [cs].
  * [39] Jianqiao Zheng, Sameera Ramasinghe, and Simon Lucey.  Rethinking Positional Encoding, October 2021.  arXiv:2107.02561 [cs].
  * [40] Yu Sun, Jiaming Liu, Mingyang Xie, Brendt Wohlberg, and Ulugbek S. Kamilov.  CoIL: Coordinate-based Internal Learning for Imaging Inverse Problems, February 2021.  arXiv:2102.05181 [eess].
  * [41] Diederik P. Kingma and Jimmy Ba.  Adam: A Method for Stochastic Optimization, January 2014.  arXiv:1412.6980 [cs].
  * [42] Matteo Ravasi and Ivan Vasconcelos.  Pylops—a linear-operator python library for scalable algebra and optimization.  SoftwareX, 11:100361, 2020.
  * [43] Norman Ricker.  Further developments in the wavelet theory of seismogram structure*.  Bulletin of the Seismological Society of America, 33(3):197–228, July 1943.
  * [44] Matteo Ravasi and Claire Birnie.  A joint inversion-segmentation approach to assisted seismic interpretation.  Geophysical Journal International, 28:893–912, 2022.
  * [45] Juan Romero, Nick Luiken, and Matteo Ravasi.  Seeing through the co2 plume: Joint inversion-segmentation of the sleipner 4d seismic data set.  The Leading Edge, 42:446–516, 2023.
  * [46] Ralph C. Smith.  Uncertainty Quantification: Theory, Implementation, and Applications.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2013.
  * [47] Yarin Gal and Zoubin Ghahramani.  Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.  In Proceedings of The 33rd International Conference on Machine Learning, pages 1050–1059. PMLR, June 2016.  ISSN: 1938-7228.
  * [48] Alex Kendall and Yarin Gal.  What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision?, October 2017.  arXiv:1703.04977 [cs].
  * [49] Enrico Schiassi, Carl Leake, Mario De Florio, Hunter Johnston, Roberto Furfaro, and Daniele Mortari.  Extreme Theory of Functional Connections: A Physics-Informed Neural Network Method for Solving Parametric Differential Equations, May 2020.  arXiv:2005.10632 [physics, stat].
  * [50] Donald A. Herron.  Seismic Intepretation.  Society of Exploration Geophysicists, 2011.
  * [51] Peter Lindstrom.  Fixed-rate compressed floating-point arrays.  IEEE Transactions on Visualization and Computer Graphics, 20(12):2674–2683, 2014.
  * [52] Altay Sansal, Sribharath Kainkaryam, Ben Lasscock, and Alejandro Valenciano.  MDIO: Open-source format for multidimensional energy data.  The Leading Edge, 42(7):465–473, July 2023.  Publisher: Society of Exploration Geophysicists.
  * [53] Keiti Aki and Paul G. Richards.  Quantitative Seismology.  W. H. Freeman and Company., 2002.
  * [54] K. Zoeppritz.  Uber reflexion und durchgang seismischer wellen durch unstetigkeitsflächen.  Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-physikalische Klasse, page 66–84, 1919.
  * [55] Gary S. Martin, Robert Wiley, and Kurt J. Marfurt.  Marmousi2 : An elastic upgrade for Marmousi.  The Leading Edge, 25(2):156–166, 02 2006.
  * [56] J E Downton.  Seismic parameter estimation from AVO inversion.  PhD thesis, University of Calgary, 2005.

##  Appendix A Samples from UQ

Figure 9 displays 12 of the 100 realizations produced as part of the
uncertainty quantification experiment conducted using IntraSeismic with Monte
Carlo dropout (MCD) on the Marmousi model [25]. Overall, the differences
between the different models are not readily apparent. However, a detailed
examination of each image reveals small variations at pixel resolution within
the high-impedance structures, which are the areas of the largest uncertainty.

##  Appendix B Ablation study of IntraSeismic loss function components

A quantitative analysis has been conducted to investigate the impact that
every loss term in the IntraSeismic loss function has on the final inverted
acoustic impedance model. To summarize, the loss function is composed of three
main terms as follows (for more details, see Section 3):

$\displaystyle\operatorname*{arg\,min}_{\Theta}\frac{1}{2}\underbrace{\left\|\mathbf{G}\boldsymbol{m}_{\Theta}-\mathbf{d}\right\|_{2}^{2}}_{\text{$\mathcal{L}_{d}$}}+\underbrace{\alpha\left\|\boldsymbol{m}_{\Theta}\right\|_{TV}}_{\text{$\mathcal{L}_{TV}$}}+\underbrace{\beta\left\|\mathbf{F}_{\Theta}(\boldsymbol{x})\right\|_{1}}_{\text{$\mathcal{L}_{m}$}}$
(A.9)

The impact of each loss term and various combinations of loss terms on the
final inverted model is shown in Figure 10. In the first scenario (Figure
10b), employing only the data loss term ($\mathcal{L}_{d}$), resulted in the
lowest SNR of 16.28. The model reconstructed in this case was of lower
resolution compared to the ground truth. In the second scenario (Figure 10c),
which combined the data loss with the MLP outputs L1 loss
($\mathcal{L}_{d}+\mathcal{L}_{m}$), an improvement in SNR was observed
compared to using the data loss alone. However, the resulting reconstructed
model appeared to contain a clear imprint of the noise in the data. The third
scenario (Figure 10d), involving the sum of the data loss and the TV loss on
the final model output ($\mathcal{L}_{d}+\mathcal{L}_{TV}$), yielded a high-
resolution model estimate. This estimate accentuated discontinuities and sharp
contrasts, but as observed in the profile view, it uniformly deviated from the
background model, thus not honoring the prior knowledge. In the fourth
scenario (Figure 10e), where all loss terms were combined
($\mathcal{L}_{d}+\mathcal{L}_{m}+\mathcal{L}_{TV}$), the reconstruction
achieved the highest SNR among all methods tested. This was reflected in a
high-resolution model estimate.

Additionally, we explored the possibility to use an $L_{2}$ norm with a
dynamically varying weight ($\beta$) across iterations ($k$) to regularize the
MLP output, as opposed to the previously described $L_{1}$ norm:
$\mathcal{L}_{m}=\beta_{k}\left\|\mathbf{F}_{\Theta}(\boldsymbol{x})\right\|_{2}^{2}$.
This regularization strategy aims to prioritize the fitting of the MLP output
to the background model
($\boldsymbol{m}_{\Theta}=\boldsymbol{m}_{\text{back}}$) in the initial
iterations. As the weight $\beta$ decreases over time, the network
progressively adapts to capture and fit the high-resolution features of the
model. The effectiveness of this approach is illustrated in Figure 10 f. These
results closely mirror those achieved with the $L_{1}$ norm, as depicted in
Figure 10 e, although with lower SNR.

##  Appendix C IntraSeismic for pre-stack seismic data

Angle-dependent seismic data can be modeled using a generalized form of the
convolutional model described in Section 2. Here, the data $d(t,\theta)$ can
be represented as the convolution of a source function or wavelet $w(t)$ and
the angle-dependent earth P-wave reflectivity series $r_{pp}(t,\theta)$:

$d(t,\theta)=w(t)*r_{pp}(t,\theta)$ (A.10)

The Aki-Richards approximation, which builds on the foundational Zoeppritz
equation [53, 54], expresses $r_{pp}$ as follows:

$r_{pp}(t,\theta)=\sum_{i=1}^{3}c_{i}(t,\theta)\frac{d}{dt}\log\left(m_{i}(t)\right),$
(A.11)

where the mixing coefficients $c_{i}$ are given by

$\displaystyle c_{1}(t,\theta)=\frac{1}{2}\left(1+\tan^{2}(\theta)\right)$
(A.12) $\displaystyle
c_{2}(t,\theta)=-4\left(\frac{\bar{V}_{s}(t)}{\bar{V}_{p}(t)}\right)^{2}\sin^{2}(\theta)$
$\displaystyle
c_{3}(t,\theta)=\frac{1}{2}-2\left(\frac{\bar{V}_{s}(t)}{\bar{V}_{p}(t)}\right)^{2}\sin^{2}(\theta)$

and

$\mathbf{m}(t)=\left(V_{p}(t),V_{s}(t),\rho(t)\right),$ (A.13)

where $\theta$ represents the incidence angle of the seismic ray-path, $V_{p}$
is the P-wave velocity, $V_{s}$ the S-wave velocity, and $\rho$ is the density
of the medium of interest. Note that equation 3, also known as the post-stack
seismic modeling equation, pertains to the special scenario where $\theta=0$,
a condition termed as zero-offset data. Combining equations A.10 and A.11, the
pre-stack seismic inverse problem can be expressed as:

$\mathbf{m^{*}}=\operatorname*{arg\,min}_{\mathbf{m}}\frac{1}{2}\left\|\mathbf{WCDm}-\mathbf{d}\right\|_{2}^{2}+\mathcal{R}(\mathbf{m}),$
(A.14)

where C is a matrix that contains the mixing coefficients.

IntraSeismic can be seamlessly extended to directly invert pre-stack seismic
data, simultaneously retrieving elastic properties such as $V_{p}$, $V_{s}$,
and $\rho$. To showcase the applicability of IntraSeismic to pre-stack
inversion, we use a portion of the synthetic elastic Marmousi model [55],
shown in Figures 11 a-c. Figures 11 d-f illustrate the background models used
as starting guesses for the inversion, which are obtained by smoothing the
corresponding ground truth models. Figures 11 g-i present three of the twenty-
one computed common-angle gathers (i.e., pre-stack data at a fixed angle).

When using IntraSeismic to invert the three subsurface models simultaneously,
as is the case for pre-stack seismic inversion, two approaches can be
considered: 1) employing the standard IntraSeismic architecture to produce
three network outputs, which are then processed by the physical modeling
module (as depicted in Figure 1); 2) enhancing IntraSeismic with three
networks (consisting of a hash table and MLP each), which individually output
a different elastic property. Their outputs are subsequently concatenated and
integrated into the physics modeling stage.

Figure 12 presents the outcomes of the first approach, showcasing a reasonable
reconstruction of $V_{p}$ and $V_{s}$, although the latter displays a lower
SNR. The estimation of $\rho$, on the other hand, although it appears notably
noisier than the other two, achieves a higher SNR; this is attributed to the
limited value range (2000 to 2500 $km/m^{3}$), resulting in a naturally lower
Total Variation (TV) term compared to other inverted properties. Conversely,
Figure 13 depicts the results from the second approach, yielding a more
effective reconstruction of $V_{p}$ and $V_{s}$. This is evident from their
visual alignment with the ground truth models (Figures 11 a-c) and an improved
SNR compared to the first approach. However, similar to the first method, the
$\rho$ estimate is less accurate compared to $V_{p}$ and $V_{s}$, resembling a
lower-resolution version of the ground truth. This issue is well known pre-
stack seismic inversion and is not specific to IntraSeismic, owing to the
limited sensitivity of pre-stack seismic data to density variations [56].
Consequently, modifying the regularization weights for each parameter
inevitably improves one model at the cost of others. Figure 14 illustrates the
SNR evolution over iterations for each model in both approaches. Both
approaches demonstrate swift convergence for the three models, with the SNR
peaks occurring nearly simultaneously and exhibiting stable behavior after
convergence.

##  Appendix D IntraSeismic sensitivity to network’s initialization

Finally, we evaluated the impact of the initialization of the parameters in
the hash tables and MLP layers of IntraSeismic. We conducted ten separate runs
for the Marmousi model, each with a different random seed. The resulting loss
functions, displayed in Figure 15a, demonstrate a high degree of alignment
across all runs, indicating a consistent behavior regardless of the
initialization used. The evolution of the SNR for the different runs exhibits
more variability (Figure 15b). Despite this variability, the final SNR for all
runs consistently exceeded a value of 24dB, which is a notable improvement
over standard methods, as detailed in Section 4.

Figure 9: 12 realizations of the Marmousi model obtained as part of the Monte
Carlo Dropout uncertainty quantification process. Figure 10: IntraSeismic
results on Marmousi model inversion using different loss functions: a) Ground
truth Marmousi model, b) data loss alone $\mathcal{L}_{d}$, c) data loss and
L1 loss of the IntraSeismic MLP output $\mathcal{L}_{d}+\mathcal{L}_{m}$, d)
data loss and TV regularization of the final model output
$\mathcal{L}_{d}+\mathcal{L}_{TV}$, e) full IntraSeismic loss
$\mathcal{L}_{d}+\mathcal{L}_{m1}+\mathcal{L}_{TV}$ being $\mathcal{L}_{m}$
the L1 norm of the MLP output, and f) full IntraSeismic loss
$\mathcal{L}_{d}+\mathcal{L}_{m2}+\mathcal{L}_{TV}$ being $\mathcal{L}_{m2}$
the L2 norm of the MLP output. Figure 11: Pre-stack Marmousi dataset: a-c
ground truth parameters $V_{p},V_{s}$ and $\rho$, d-f background models, g-i
partial stacks for angles $\theta=0$, $\theta=20$ and $\theta=40$. Figure 12:
Pre-stack inversion results of a portion of the elastic Marmousi model using
IntraSeismic: a) P-wave velocity $V_{p}$, b) S-wave velocity $V_{s}$, and c)
density $\rho$. Figure 13: Pre-stack inversion results of a portion of the
elastic Marmousi model using IntraSeismic with three separate networks: a)
P-wave velocity $V_{p}$, b) S-wave velocity $V_{s}$, and c) density $\rho$.
Figure 14: SNR evolution through iterations for each elastic parameter
($V_{p},V_{s},\rho$) using a) standard IntraSeismic, b) IntraSeismic with
three networks. Figure 15: Results of 10 random initializations of
IntraSeismic: a) loss function curves and b) SNR curves.