Seismic imaging is an ill-posed inverse problem that is challenged by noisy data and modeling inaccuracies—due to errors in the background squared-slowness model. Uncertainty quantification is essential for determining how variability in the background models affects seismic imaging. Due to the costs associated with the forward Born modeling operator as well as the high dimensionality of seismic images, quantification of uncertainty is computationally expensive. As such, the main contribution of this work is a survey-specific Fourier neural operator surrogate to velocity continuation that maps seismic images associated with one background model to another virtually for free. While being trained with only 200 background and seismic image pairs, this surrogate is able to accurately predict seismic images associated with new background models, thus accelerating seismic imaging uncertainty quantification. We support our method with a realistic data example in which we quantify seismic imaging uncertainties using a Fourier neural operator surrogate, illustrating how variations in background models affect the position of reflectors in a seismic image.

Seismic imaging involves estimating the short-wavelength component of the Earth’s subsurface squared-slowness model—known as the seismic image—given shot records and an estimation of the smooth background squared-slowness model. This linearized imaging problem is challenged by the computationally expensive forward operator as well as presence of measurements noise, linearization errors, modeling errors, and the nontrivial nullspace of the linearized forward Born modeling operator (Lambaré et al., 1992; Schuster, 1993; Nemeth et al., 1999). These challenges highlight the importance of uncertainty quantification (UQ) in seismic imaging, where instead of finding one seismic image estimate, a distribution of seismic images is obtained that explains the observed data (Tarantola, 2005), consequently reducing the risk of data overfit and enabling UQ (Malinverno and Parker, 2006; Martin et al., 2012; Ray et al., 2017; Zhu et al., 2016; Fang et al., 2018; G. K. Stuart et al., 2019; Herrmann et al., 2019; Siahkoohi et al., 2020a, 2020b, 2021; Siahkoohi and Herrmann, 2021).

The seismic imaging uncertainty can be attributed to two main sources (Thore et al., 2002; Osypov et al., 2013; Ely et al., 2018): (1) errors in the data, which include measurement and linearization errors; and (2) modeling errors, which include errors in the estimation of the background squared-slowness model. In this paper, we focus on uncertainties with respect to the background model as it has the main contributing factor to imaging uncertainty due to its effect on reflector positioning (Fomel and Landa, 2014; O. V. Poliannikov and Malcolm, 2016; Ely et al., 2018). To this end, we assume there are no measurement or linearization errors and consider the map from shot records to the seismic image to be represented by the deterministic reverse-time migration (RTM) algorithm. In this setup, quantifying the uncertainty in seismic imaging—due to errors in the background model—involves computing numerous RTMs with all the background model posterior samples (Zhu et al., 2016; Fang et al., 2018). Due to the high-dimensionality of seismic images, the number of seismic imaging posterior samples required to obtain accurate estimations of the imaging posterior moments is large, making UQ with respect to the background model computationally expensive for 2D and unfeasible for large 3D problems. These costs can be alleviated via velocity continuation methods (Fomel, 2003a, 2003b; Duchkov and De Hoop, 2009; Leeuwen and Herrmann, 2012; M. Yang et al., 2021) that map seismic images associated with one background model to another, without directly solving the imaging problem. Velocity continuation methods are designed to be frugal compared to computing a new RTM for each new background model, which makes them suitable for large-scale seismic imaging uncertainty quantification (Fomel and Landa, 2014; O. V. Poliannikov and Malcolm, 2016).

While several velocity continuation methods have been proposed (Fomel, 2003a, 2003b; Duchkov and De Hoop, 2009; Leeuwen and Herrmann, 2012; M. Yang et al., 2021), they typically involve solving partial differential equations (PDEs), which can be still costly in the context of seismic imaging UQ. To address this challenge, we propose a neural network surrogate for velocity continuation that is capable of mapping seismic images associated with one background model to another with negligible computational cost. Motivated by the success of Fourier neural operators (FNOs, Z. Li et al., 2021) in approximating the solution operator of PDEs (Toledo-Marín et al., 2021; Konuk and Shragge, 2021; N. Kovachki et al., 2021), we chose them as the architecture for our neural network surrogate. Due to our main interest in accelerating velocity continuation in the context of UQ, we train a survey-specific FNO that acts as a surrogate for velocity continuation for the specific survey at hand. This choice, while not offering generalizations across different survey areas in the Earth, can speed up seismic imaging UQ for the survey at hand, which involves computing seismic images associated with many background model posterior samples. We show that the FNO can be trained using \(200\) pairs of background and seismic image pairs, making the survey-specific training procedure computationally viable. To scale this method to industry size problems, transfer learning (Yosinski et al., 2014) can further reduce the upfront costs of training the FNO. After training, the FNO can be used to obtain samples from the imaging posterior almost free of cost.

In the next section, we describe seismic imaging by introducing the forward Born modeling operator, through which the seismic image relates to the background model. Next, we define velocity continuation, followed by describing our proposed FNO-based approach for accelerating seismic imaging UQ. Finally, we evaluate the performance of the trained neural operator on a realistic dataset, and we demonstrate how the imaging uncertainty affects the positioning of the reflectors in the seismic image.

We introduce a deep-network surrogate for velocity continuation in order to enable faster quantification of uncertainty—due to errors in the background velocity model—in seismic imaging. We begin with an introduction to seismic imaging and the linearized forward model associated with it.

The inverse problem that we tackle involves the process of estimating the short-wavelength component of the Earth’s unknown subsurface squared-slowness model given measurements recorded at the surface. This problem, also known as seismic imaging, can be formulated as a linear inverse problem by linearizing the nonlinear relationship between shot records and the squared-slowness model, governed by the wave-equation. In its simplest acoustic form, the linearization with respect to the slowness model—around a background squared slowness model \(\B{m}_0\)—leads to a linear inverse problem for estimating the ground truth seismic image \(\delta \B{m}^{\ast}\) with the following forward model, \[ \begin{equation} \delta \B{d}_i = \B{J}(\B{m}_0, \B{q}_i) \delta \B{m}^{\ast}. \label{linear-fwd-op} \end{equation} \] In the above expression, \(\delta \B{d}=\left\{\delta \B{d}_{i}\right\}_{i=1}^{n_s}\) are \(n_s\) linearized shot records and \(\B{J}(\B{m}_0, \B{q}_i)\) represents the linearized Born scattering operator. This operator is parameterized by the source signature \(\B{q}_{i}\) and the background squared-slowness model \(\B{m}_0\), which is typically estimated in the previous inversion steps (Zhu et al., 2016; Fang et al., 2018). In this work, we focus on potential inaccuracies in the background model, the main source of variability and uncertainty in seismic imaging (Fomel and Landa, 2014; O. V. Poliannikov and Malcolm, 2016; Ely et al., 2018). Therefore, we assume there are not measurement and linearization errors, which leads to a deterministic mapping from \(\left\{\delta \B{d}_{i}\right\}_{i=1}^{n_s}\) to \(\delta \B{m}\) for a given background model. We use RTM for this deterministic map, which involves applying the adjoint linearized Born scattering operator to the linearized shot records for all source experiments, \[ \begin{equation} \delta \B{m}_{\text{RTM}} = \sum_{i=1}^{n_s} \B{J}(\B{m}_0, \B{q}_i)^{\top} \delta\B{d}_i. \label{rtm} \end{equation} \] Since this map is deterministic, the uncertainty in the background model can be translated into uncertainty in seismic imaging by evaluating Equation \(\ref{rtm}\) for the background models sampled from the posterior distribution \(p(\B{m}_0\mid \B{d})\), where \(\B{d}\) represents the shot records with no linearization (the field data) (Zhu et al., 2016; Fang et al., 2018). Bayesian inference involving high-dimensional posterior distributions, for example seismic imaging and full-waveform inversion, requires many samples from the posterior distribution for accurate Monte Carlo approximation of high-dimensional integrals (Gelman and Rubin, 1992). As a result, the mapping in Equation \(\ref{rtm}\) must be evaluated over numerous background model samples from \(p(\B{m}_0\mid \B{d})\), which is computationally expensive due to costs of applying \(\B{J}(\B{m}_0,\B{q}_i)^{\top},\,i=1,\ldots,n_s\). To address this computational challenge, our proposed method involves a FNO-based velocity continuation approach that is capable of mapping seismic images associated with one background model to another, effectively replacing a costly demingration-migration. We describe this approach in the next section.

At its core, velocity continuation is a process which alters a seismic image based on changes in the background model (Fomel, 2003a, 2003b; Leeuwen and Herrmann, 2012; M. Yang et al., 2021). Using this technique, an initial seismic image associated with an initial background model \(\B{m}_{\text{init}}\) is altered to approximate the seismic image associated with a target background model \(\B{m}_{\text{target}}\) without computing Equation \(\ref{rtm}\). This process can be interpreted as a map (Duchkov and De Hoop, 2009): \[ \begin{equation} \mathcal{T}_{\left(\B{m}_{\text{init}}, \B{m}_{\text{target}}\right)}: \delta \mathcal{M} \to \delta \mathcal{M}, \label{velocity-continuation} \end{equation} \] where \(\delta \mathcal{M}\) denotes the space of seismic images, and the velocity continuation map is parameterized by \(\B{m}_{\text{init}}\) and \(\B{m}_{\text{target}}\). We take advantage of recent advances in deep learning to train a surrogate model for velocity continuation that approximates \(\mathcal{T}_{\left(\B{m}_{\text{init}},\B{m}_{\text{target}}\right)} = \B{J}(\B{m}_{\text{target}},\B{q}_i)^{\top}(\B{J}(\B{m}_{\text{init}}, \B{q}_i)^{\top})^{\dagger}\), which can be evaluated virtually for free instead of solving four PDEs. Through this approach, we are able to significantly accelerate velocity continuation, trading wave-equation solves for a simple neural network inference. This is of considerable importance in the context of seismic imaging UQ. Before demonstrating the benefits of our approach, we introduce FNOs in the context of velocity continuation.

In light of their success in learning mesh-free solution operators to PDEs (Toledo-Marín et al., 2021; Konuk and Shragge, 2021; N. Kovachki et al., 2021), we choose FNOs (Z. Li et al., 2021) as a surrogate model for velocity continuation—a process that can be interpreted as a double linearized PDE solve. The main components of FNOs are the Fourier layers, which involve a Fourier transform over the spatial dimensions of their input, followed by a learned pointwise multiplication and an inverse Fourier transform. These layers act as long-kernel convolutional layer, akin to pseudo-spectral methods, which explains the representation power of FNOs (N. Kovachki et al., 2021). To train a FNO as a surrogate model for velocity continuation, we define it as a map \[ \begin{equation} \mathcal{G}_{\B{w}} : \mathcal{M} \times \delta \mathcal{M} \to \delta \mathcal{M}, \label{fno-map} \end{equation} \] where \(\mathcal{G}_{\B{w}}\) denotes the FNO with weights \(\B{w}\), and \(\mathcal{M}\) is the space of background models. In words, we design \(\mathcal{G}_{\B{w}}\) as a neural network that takes as input the target background model and the initial seismic image, and outputs the target seismic image. This choice is analogous to the structure of the velocity continuation map \(\mathcal{T}_{\left(\B{m}_{\text{init}},\B{m}_{\text{target}}\right)}\) (cf. Equation \(\ref{velocity-continuation}\)), except that we fix the initial background model \(\B{m}_{\text{init}}\) to an arbitrary background model posterior sample, hence making the dependence of \(\mathcal{G}_{\B{w}}\) to \(\B{m}_{\text{init}}\) implicit. With this choice of inputs and output for \(\mathcal{G}_{\B{w}}\), the FNO’s task is to perturb the input initial seismic image according to the provided target background model in order to predict the target seismic image.

Due to our interest in accelerating velocity continuation in the context of UQ, we train a survey-specific FNO that acts as a surrogate for velocity continuation for the specific survey at hand. This choice, while not offering generalizations across different survey areas in the Earth, can speed up seismic imaging UQ for the survey at hand, which involves computing seismic images associated with many posterior background model samples \(\B{m}_0 \sim p(\B{m}_0\mid \B{d})\). By training the FNO on a small number of these background model and seismic image pairs, i.e., approximately \(200\), we can accelerate the velocity continuation process for the rest of the background model posterior samples while limiting the risk of introducing generalization errors due to the strong heterogeneity of Earth. To achieve this, we construct a set of \(N\) training input-output pairs in the form of \[ \begin{equation} \bigg\{ \Big ( \big(\B{m}_{0}^{(i)},\delta \B{m}_{\text{init}} \big),\delta \B{m}_{\text{RTM}}^{(i)} \Big) \ \Big | \ i=1,\ldots,N \bigg \}, \label{fno-training-pairs} \end{equation} \] where \((\B{m}_{0}^{(i)},\delta \B{m}_{\text{init}})\) is the input target background and initial seismic image training pair, and \(\delta \B{m}_{\text{RTM}}^{(i)}\) is the associated target seismic image. Training involves minimizing the squared \(\ell_2\)-norm of the difference between the FNO output and the target seismic image with respect to FNO weights, \[ \begin{equation} \B{w}^{\ast} = \argmin_{\B{w}}\, \frac{1}{N} \sum_{i=1} \| \mathcal{G}_{\B{w}}(\B{m}_{0}^{(i)}, \delta \B{m}_{\text{init}}) - \delta \B{m}_{\text{RTM}}^{(i)} \|_2^2. \label{fno-train} \end{equation} \] We solve optimization problem \(\ref{fno-train}\) with the Adam stochastic optimization algorithm (Kingma and Ba, 2014). After training, the trained FNO approximates the velocity continuation map for a fixed homogeneous initial background model, i.e., \[ \begin{equation} \mathcal{G}_{\B{w}^{\ast}}(\B{m}_{\text{target}}, \delta \B{m}_{\text{init}}) \approx \mathcal{T}_{\left(\B{m}_{\text{init}}, \B{m}_{\text{target}}\right)} (\delta \B{m}_{\text{init}}) = \delta \B{m}_{\text{target}}, \label{fno-approx} \end{equation} \] where \(\B{m}_{\text{init}}\) is fixed to an arbitrary background model posterior sample and \(\delta \B{m}_{\text{target}}\) denotes the seismic image associated with \(\delta \B{m}_{\text{init}}\). Given background model samples \(\B{m}_0 \sim p(\B{m}_0\mid \B{d})\) as input, the FNO outputs samples from the imaging posterior, \(p(\B{m}_0\mid \delta \B{d})\). This process accelerates seismic imaging uncertainty quantification as no further RTMs (Equation \(\ref{rtm}\)) need to be computed. Training the FNO with as little as \(200\) training pairs makes the survey-specific training procedure computationally viable. To further accelerate the process, training of the FNO can be started as the same time as the background model posterior sampling phase, using the already collected posterior samples as training data. In the next section, we show the results of approximating the velocity continuation map with a FNO via a quasi-real seismic experiment.

The purpose of the presented numerical experiments here is to demonstrate the ability to approximate the velocity continuation map (Equation \(\ref{velocity-continuation}\)) using a FNO. We further show how this trained FNO can be used for UQ by showing the effect of imaging uncertainty on the positioning of the reflectors. We begin by describing the training setup, including the seismic acquisition geometry.

Our examples involve imaging a 2D subset of the Parihaka (Veritas, 2005; WesternGeco., 2012) prestrack Kirchhoff migration field dataset. We use this 2D section to create linearized data according to the linear forward model in Equation \(\ref{linear-fwd-op}\), where we consider no measurement and linearization errors. The model is discretized on a \(12.5\,\mathrm{m}\) vertical and \(20\,\mathrm{m}\) horizontal grid, and the data is acquired with \(102\) equally spaced sources and \(204\) fixed receivers. We use a Ricker wavelet with a central frequency of \(30\, \mathrm{Hz}\) as the source signature. For presentation purposes, we augment a \(125\, \mathrm{m}\) water column on top of these models to limit the near source imaging artifacts. We create background models by assuming access to an oracle, which provides geologically consistent background models with the 2D section. Given \(200\) background models, we migrate the simulated seismic data via Equation \(\ref{rtm}\) to obtain the corresponding seismic images. Figures 1a and 1b show the vertical profile of five randomly selected background and seismic images at \(2.50\,\mathrm{km}\), respectively. These images, displaying strong amplitude and phase differences, highlight the importance of quantifying the uncertainty in imaging when dealing with errors in the background model. To reduce the costs associated with UQ, we train the FNO surrogate using the \(200\) training pairs (cf. Equation \(\ref{fno-training-pairs}\)) for \(500\) epochs with the Adam optimizer (Kingma and Ba, 2014).

To evaluate the accuracy of the trained FNO in predicting seismic images, we compare its output to the target seismic image, i.e., the RTM image obtained via Equation \(\ref{rtm}\) when using the target background model (Figure 2a) to parameterize the Born scattering forward operator. This is conducted over the testing dataset, which is derived using the same procedure as the training dataset using the background model creating oracle. Figures 2b and 2c show the target and predicted seismic images, respectively, and Figures 2d includes the difference between them. We observe that the network has accurately predicted the target image, where errors are mostly due to amplitude differences. The accuracy of the prediction in terms of phase and amplitude can be further confirmed by focusing on two vertical profiles, at \(2.50\,\mathrm{km}\) (Figure 2e) and \(4.125\,\mathrm{km}\) (Figure 2f) horizontal locations, which include regions with fault and torturous reflectors. We use the trained FNO for UQ by providing testing background models and predicting the associated seismic images. We visualize the obtained uncertainties by showing the variability in the location of reflectors, determined via an automatic horizon tracking software (X. Wu and Fomel, 2018). We pass the seismic images predicted by the FNO to the horizon tracker for \(25\) selected horizons. As a result, we obtain multiple instances of each horizons, from which we compute pointwise mean and standard deviations. Figure 3 indicates the result where the solid lines correspond to the mean among different instances of each horizons and shaded areas indicate the mean plus and minus the pointwise standard deviation. These shaded areas indicate uncertainties in the location of the reflectors, which are due to the variability in the background models. As expected, we find a general increase of uncertainty with depth. We also observe that the areas of high uncertainty are correlated with areas of poor illumination, faults and tortuous reflectors.

In this example, we use Devito (Luporini et al., 2018; Louboutin et al., 2019) for the wave-equation based simulations. We based our PyTorch FNO implementation on the original implementation. The code to reproduce our results are made available on GitHub.

Quantifying the uncertainty in seismic imaging due to errors in the background model involves solving many seismic imaging problems that vary in the parameterization of the background model of the forward operator. To reduce the computational cost of this process—mainly due to the computational costs of the forward operator—we proposed to train a survey-specific Fourier neural operator surrogate that mimics velocity continuation. This surrogate model maps seismic images associated with one background model to another virtually for free, which has the benefit of accelerating uncertainty quantification. We showed that this surrogate model can be trained with as few as \(200\) training pairs while still providing a good seismic image prediction accuracy. Further research is required in training a reliable global surrogate, being able to generalize across other survey areas and more realistic physics.

This research was carried out with the support of Georgia Research Alliance and partners of the ML4Seismic Center.

Duchkov, A. A., and De Hoop, M. V., 2009, Velocity continuation in the downward continuation approach to seismic imaging: Geophysical Journal International, **176**, 909–924.

Ely, G., Malcolm, A., and Poliannikov, O. V., 2018, Assessing uncertainties in velocity models and images with a fast nonlinear uncertainty quantification method: GEOPHYSICS, **83**, R63–R75. doi:10.1190/geo2017-0321.1

Fang, Z., Silva, C. D., Kuske, R., and Herrmann, F. J., 2018, Uncertainty quantification for inverse problems with weak partial-differential-equation constraints: GEOPHYSICS, **83**, R629–R647. doi:10.1190/geo2017-0824.1

Fomel, S., 2003a, Time-migration velocity analysis by velocity continuation: Geophysics, **68**, 1662–1672.

Fomel, S., 2003b, Velocity continuation and the anatomy of residual prestack time migration: Geophysics, **68**, 1650–1661.

Fomel, S., and Landa, E., 2014, Structural uncertainty of time-migrated seismic images: Journal of Applied Geophysics, **101**, 27–30.

Gelman, A., and Rubin, D. B., 1992, Inference from Iterative Simulation Using Multiple Sequences: Statistical Science, **7**, 457–472. Retrieved from http://www.jstor.org/stable/2246093

Herrmann, F. J., Siahkoohi, A., and Rizzuti, G., 2019, Learned imaging with constraints and uncertainty quantification: Neural Information Processing Systems (NeurIPS) 2019 Deep Inverse Workshop. Retrieved from https://arxiv.org/pdf/1909.06473.pdf

Kingma, D. P., and Ba, J., 2014, Adam: A method for stochastic optimization: ArXiv Preprint ArXiv:1412.6980. Retrieved from https://arxiv.org/pdf/1412.6980.pdf

Konuk, T., and Shragge, J., 2021, Physics-guided deep learning using fourier neural operators for solving the acoustic vTI wave equation: In 82nd eAGE annual conference & exhibition (Vol. 2021, pp. 1–5). European Association of Geoscientists & Engineers.

Kovachki, N., Lanthaler, S., and Mishra, S., 2021, On universal approximation and error bounds for fourier neural operators: Journal of Machine Learning Research, **22**, Art–No.

Lambaré, G., Virieux, J., Madariaga, R., and Jin, S., 1992, Iterative asymptotic inversion in the acoustic approximation: Geophysics, **57**, 1138–1154.

Leeuwen, T. van, and Herrmann, F. J., 2012, Wave-equation Extended Images-Computation and Velocity Continuation: In 74th eAGE conference and exhibition incorporating eUROPEC 2012 (pp. cp–293). European Association of Geoscientists & Engineers.

Li, Z., Kovachki, N. B., Azizzadenesheli, K., liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A., 2021, Fourier Neural Operator for Parametric Partial Differential Equations: In International conference on learning representations. Retrieved from https://openreview.net/forum?id=c8P9NQVtmnO

Louboutin, M., Lange, M., Luporini, F., Kukreja, N., Witte, P. A., Herrmann, F. J., … Gorman, G. J., 2019, Devito (v3.1.0): An embedded domain-specific language for finite differences and geophysical exploration: Geoscientific Model Development, **12**, 1165–1187. doi:10.5194/gmd-12-1165-2019

Luporini, F., Lange, M., Louboutin, M., Kukreja, N., Hückelheim, J., Yount, C., … Gorman, G. J., 2018, Architecture and performance of devito, a system for automated stencil computation: CoRR, **abs/1807.03032**. Retrieved from http://arxiv.org/abs/1807.03032

Malinverno, A., and Parker, R. L., 2006, Two ways to quantify uncertainty in geophysical inverse problems: GEOPHYSICS, **71**, W15–W27. doi:10.1190/1.2194516

Martin, J., Wilcox, L. C., Burstedde, C., and Ghattas, O., 2012, A Stochastic Newton MCMC Method for Large-scale Statistical Inverse Problems with Application to Seismic Inversion: SIAM Journal on Scientific Computing, **34**, A1460–A1487. Retrieved from http://epubs.siam.org/doi/abs/10.1137/110845598

Nemeth, T., Wu, C., and Schuster, G. T., 1999, Least-squares migration of incomplete reflection data: GEOPHYSICS, **64**, 208–221. doi:10.1190/1.1444517

Osypov, K., Yang, Y., Fournier, A., Ivanova, N., Bachrach, R., Yarman, C. E., … Woodward, M., 2013, Model-uncertainty quantification in seismic tomography: Method and applications: Geophysical Prospecting, **61**, 1114–1134.

Poliannikov, O. V., and Malcolm, A. E., 2016, The effect of velocity uncertainty on migrated reflectors: Improvements from relative-depth imaging: Geophysics, **81**, S21–S29.

Ray, A., Kaplan, S., Washbourne, J., and Albertin, U., 2017, Low frequency full waveform seismic inversion within a tree based Bayesian framework: Geophysical Journal International, **212**, 522–542. doi:10.1093/gji/ggx428

Schuster, G. T., 1993, Least-squares cross-well migration: In SEG technical program expanded abstracts 1993 (pp. 110–113). Society of Exploration Geophysicists.

Siahkoohi, A., and Herrmann, F. J., 2021, Learning by example: fast reliability-aware seismic imaging with normalizing flows: In First international meeting for applied geoscience & energy (pp. 1580–1585). Society of Exploration Geophysicists.

Siahkoohi, A., Rizzuti, G., and Herrmann, F. J., 2020a, A deep-learning based bayesian approach to seismic imaging and uncertainty quantification: 82nd EAGE Conference and Exhibition 2020. Retrieved from https://slim.gatech.edu/Publications/Public/Submitted/2020/siahkoohi2020EAGEdlb/siahkoohi2020EAGEdlb.html

Siahkoohi, A., Rizzuti, G., and Herrmann, F. J., 2020b, Uncertainty quantification in imaging and automatic horizon tracking—a Bayesian deep-prior based approach: SEG technical program expanded abstracts 2020. doi:10.1190/segam2020-3417560.1

Siahkoohi, A., Rizzuti, G., and Herrmann, F. J., 2021, Deep bayesian inference for seismic imaging with tasks: ArXiv Preprint ArXiv:2110.04825.

Stuart, G. K., Minkoff, S. E., and Pereira, F., 2019, A two-stage Markov chain Monte Carlo method for seismic inversion and uncertainty quantification: GEOPHYSICS, **84**, R1003–R1020. doi:10.1190/geo2018-0893.1

Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM. doi:10.1137/1.9780898717921

Thore, P., Shtuka, A., Lecour, M., Ait-Ettajer, T., and Cognot, R., 2002, Structural uncertainties: Determination, management, and applications: Geophysics, **67**, 840–852.

Toledo-Marín, J. Q., Fox, G., Sluka, J. P., and Glazier, J. A., 2021, Deep learning approaches to surrogates for solving the diffusion equation for mechanistic real-world simulations: Frontiers in Physiology, **12**.

Veritas, 2005, Parihaka 3D Marine Seismic Survey - Acquisition and Processing Report: No. New Zealand Petroleum Report 3460. New Zealand Petroleum & Minerals, Wellington.

WesternGeco., 2012, Parihaka 3D PSTM Final Processing Report: No. New Zealand Petroleum Report 4582. New Zealand Petroleum & Minerals, Wellington.

Wu, X., and Fomel, S., 2018, Least-squares horizons with local slopes and multi-grid correlations: GEOPHYSICS, **83**, IM29–IM40. doi:10.1190/geo2017-0830.1

Yang, M., Graff, M., Kumar, R., and Herrmann, F. J., 2021, Low-rank representation of omnidirectional subsurface extended image volumes: Geophysics, **86**, S165–S183.

Yosinski, J., Clune, J., Bengio, Y., and Lipson, H., 2014, How Transferable Are Features in Deep Neural Networks? In Proceedings of the 27th international conference on neural information processing systems (pp. 3320–3328). Retrieved from http://dl.acm.org/citation.cfm?id=2969033.2969197

Zhu, H., Li, S., Fomel, S., Stadler, G., and Ghattas, O., 2016, A bayesian approach to estimate uncertainty for full-waveform inversion using a priori information from depth migration: Geophysics, **81**, R307–R323.