#  NuSPAN: A Proximal Average Network for Nonuniform Sparse Model —
Application to Seismic Reflectivity Inversion

Swapnil Mache†, Praveen Kumar Pokala†, Kusala Rajendran, and Chandra Sekhar
Seelamantula † Equal contribution.Swapnil Mache is with the Centre of
Excellence in Advanced Mechanics of Materials, Indian Institute of Science
(IISc), Bangalore, India, and with the Department of Electrical Engineering,
IISc. He was formerly with the Centre for Earth Sciences, IISc. Email:
<EMAIL_ADDRESS>Kumar Pokala and Chandra Sekhar Seelamantula
are with the Department of Electrical Engineering, IISc. E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Kusala Rajendran is
with the Centre of Excellence in Advanced Mechanics of Materials, IISc. She is
Retired Professor from the Centre for Earth Sciences, IISc. Email:
<EMAIL_ADDRESS>
###### Abstract

We solve the problem of sparse signal deconvolution in the context of seismic
reflectivity inversion, which pertains to high-resolution recovery of the
subsurface reflection coefficients. Our formulation employs a nonuniform, non-
convex synthesis sparse model comprising a combination of convex and non-
convex regularizers, which results in accurate approximations of the
$\ell_{0}$ pseudo-norm. The resulting iterative algorithm requires the
proximal average strategy. When unfolded, the iterations give rise to a
learnable proximal average network architecture that can be optimized in a
data-driven fashion. We demonstrate the efficacy of the proposed approach
through numerical experiments on synthetic $1$-D seismic traces and $2$-D
wedge models in comparison with the benchmark techniques. We also present
validations considering the simulated Marmousi$2$ model as well as real $3$-D
seismic volume data acquired from the Penobscot $3$D survey off the coast of
Nova Scotia, Canada.

###### Index Terms:

Geophysics, inverse problems, seismology, seismic reflectivity inversion,
geophysical signal processing, deep learning, neural networks, algorithm
unrolling, nonconvex optimization, sparse recovery.

##  I Introduction

Reflectivity inversion, which is a seismic deconvolution problem in reflection
seismology, is the means by which one can characterize and image the layered
subsurface structure. A recorded seismic trace
$\boldsymbol{y}\in\mathbb{R}^{n}$ is modeled as the convolution of the source
pulse $\boldsymbol{h}$, assumed to be a Ricker wavelet [1], and the subsurface
reflectivity $\boldsymbol{x}\in\mathbb{R}^{n}$, given by
$\boldsymbol{y}=\boldsymbol{h}*\boldsymbol{x}+\boldsymbol{n},$ where
$\boldsymbol{n}$ is the noise and $*$ denotes convolution. Reflectivity is
modeled as a sparse vector comprising subsurface layers, each with a constant
impedance, giving rise to a piecewise-constant impedance model [2, 3]. Figure
1 illustrates the model with two shale layers sandwiching a wet sand layer
(Figure 1(a)) [4]. The reflectivity (Figure 1(c)) convolved with the source
(Figure 1(d)) is recorded at the receiver as a seismic trace (Figure 1(e)).

Figure 1:  An ideal three-layer subsurface model. The operator $*$ denotes
convolution.

The subsurface geology (Figure 1(a)) is related to reflectivity (Figure 1(c))
through acoustic impedance (Figure 1(b)), which is defined as the product of
density ($\rho$) and P-wave velocity ($V$) [2], given by

$r_{i}=\frac{\rho_{i+1}V_{i+1}-\rho_{i}V_{i}}{\rho_{i+1}V_{i+1}+\rho_{i}V_{i}}.$
(1)

Solving the ill-posed linear inverse problem of estimating the subsurface
reflectivity through the classical least-squares formulation [5] leads to
nonuniqueness issues arising out of a convolution with a bandlimited wavelet
and loss of low and high-frequency information [2, 6, 7]. The nonuniqueness
aspect can be tackled through regularization [8], for example, by imposing a
sparsity prior on the solution through the $\ell_{1}$-norm [2, 5].

Zhang and Castagna [9] solved the $\ell_{1}$-norm constrained reflectivity
inversion problem through basis-pursuit inversion (BPI) [10], using a wavelet
dictionary of odd and even reflectivity pairs. The fast iterative shrinkage-
thresholding algorithm (FISTA) [11] has been employed for reflectivity
inversion [12] along with an amplitude recovery boost through debiasing steps
of least-squares inversion [13] and “adding back the residual” [14]. The
$\ell_{1}$-norm, although convex, is not the best sparsity constraint for
reflectivity inversion, and the accurate estimation of the sparsity of seismic
reflections is challenging [15, 16]. Further, $\ell_{1}$-norm regularization
underestimates the high-amplitude components and introduces a bias in the
estimate of the sparse code $\boldsymbol{x}$ [17, 18, 19]. Both [13] and [14]
observed the attenuation of reflectivity magnitudes due to the $\ell_{1}$-norm
regularization term and adopted a post-processing debiasing step [20].

Non-convex regularization strategies have been adopted to overcome the
shortcomings of $\ell_{1}$ regularization in sparse recovery problems [19, 21,
22]. Non-convex penalties such as the smoothly clipped absolute deviation
(SCAD) [23] and minimax concave penalty (MCP) [21] have been shown to be
superior to $\ell_{1}$ and non-convex regularization approaches in the context
of inverse problems. In the context of reflectivity inversion, [15] introduced
a data-driven $\ell_{p}$-loss-$\ell_{q}$-regularization ($p=2,0<q<1$), with an
adaptive approach for choosing the optimal $q$. Recently, [24] and [25]
developed proximal gradient-descent methods based on the proximal average
strategy [26, 27], by considering composite regularization, which is a convex
combination of several sparsity-inducing regularizers.

A new class of deep neural network (DNN) architectures was introduced by [28],
which are inspired by unrolling iterative algorithms into learnable networks
[29]. Gregor and LeCun [28] proposed the learned iterative shrinkage and
thresholding algorithm (LISTA), based on unfolding the update steps of ISTA
[30] into the layers of a neural network. This class of model-based
architectures has been demonstrated to be effective in solving sparse linear
inverse problems [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44].
Recently, Monga et al. [29] reviewed algorithm unrolling and its applications
to imaging, vision and recognition, speech processing, and other signal and
image processing problems.

Recently, learning-based approaches have also been employed for solving
inversion problems in seismology and geophysics [45, 16, 4, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67],
including seismic reflectivity inversion [45, 16, 4, 47, 53, 62]. Adler et al.
reviewed the development of deep learning solutions in the context of
seismology and geophysics from approaches that enhanced the performance of
full waveform inversion [48, 49, 50, 51] to end-to-end seismic inversion [52,
53, 54, 55, 56, 57, 58, 59, 60, 61] that employed physics-guided architectures
[62, 63] and deep generative modeling [64, 65, 66, 67].

Kim and Nakata [53] employed an elementary feedforward neural network and
observed superior support recovery compared to a least-squares approach but
poor amplitude recovery. Yuan and Su [16], and earlier, [45], employed sparse
Bayesian learning (SBL) for recovering the sparse reflection coefficients by
maximizing the marginal likelihood either by using the expectation-
maximization algorithm (SBL-EM) [68, 16], or by sequentially updating the
sparsity-controlling hyperparameters through a sequential algorithm-based
approach [45]. The former [68, 16] was demonstrated as being more robust to
noise, having higher accuracy, and preserving lateral continuity in the
seismic profile.

###  I-A Motivation and Contribution

Estimating the sparsity of the reflection data is difficult [15, 16] and the
recovery of support is prioritized [1]. One can resort to data-driven
approaches, which outperform conventional techniques, especially in support
recovery [53], when the knowledge about the underlying geology is limited,
with an added computational benefit [4]. Further, as opposed to an elementary
feedforward neural network approach for seismic reflectivity inversion [53],
one can employ model-based learning frameworks [39] such as deep-unrolled
architectures [29]. Such deep neural networks, where the architecture is
informed by the inverse problem itself, can provide insights into the problem
as well as model interpretability that is critical to gaining physical
insights into the system under consideration [46]. In our previous work [69],
we demonstrated the efficacy of FirmNet [37] and LISTA-like [28] formulations
in solving the seismic reflectivity inversion problem, comparing the
approaches with baselines such as BPI [10, 9], FISTA [11, 12], and SBL-EM [68,
16]. Here, we expand the work by constructing and learning from the data [46],
a composite nonuniform sparse prior from a convex combination of weighted
counterparts of three sparsity-promoting penalties (the $\ell_{1}$ norm, MCP
[21], and SCAD [23] penalties).

The contributions of this paper are stated below.

  1. 1.

We consider the problem of seismic reflectivity inversion based on a data-
driven prior, as opposed to a pre-designated prior, within the framework of
deep-unfolding. To the best of our knowledge, data-driven priors have not been
explored for solving the problem of seismic reflectivity inversion.

  2. 2.

We propose an optimization framework for seismic reflectivity inversion based
on a composite sparse-prior, which comprises multiple weighted regularizers.
The weights are allowed to be different for each component. Such a model is
commonly referred to as a nonuniform sparse model.

  3. 3.

We develop the nonuniform proximal-averaged thresholding algorithm (NuPATA)
and its deep-unfolded version, the nonuniform sparse proximal average network
(NuSPAN) to solve the problem under consideration.

  4. 4.

We demonstrate the efficacy of the proposed network, NuSPAN, over synthetic
$1$-D and $2$-D datasets, and simulated and real datasets.

###  I-B Organization of the Paper and Notations

This paper is structured as follows. In Section III, we introduce the
nonuniform sparse model for sparse seismic reflectivity inversion. We develop
an optimization algorithm and the corresponding unrolled network for solving
the problem under consideration. Section IV explains a generalized variant of
the problem formulation discussed in Section III. In Section V, we present
experimental results and demonstrate the efficacy of the proposed proximal
average network (NuSPAN) in comparison with baselines. Conclusions are
provided in Section VI.

The notational conventions used in this paper are listed in Table I.

TABLE I: Notations

Notation | Description
---|---
$x$ | Scalar (lower-case letter)
$\boldsymbol{x}$ | Vector (lower-case bold letter)
${\left\|\boldsymbol{x}\right\|}_{p}$ |  $\ell_{p}-\mathrm{norm}$ of $\boldsymbol{x}$
$\boldsymbol{x}^{(k)}$ |  $\boldsymbol{x}$ at the $k^{th}$ iteration of an algorithm
$x_{i}$ |  $i^{th}$ element of $\boldsymbol{x}$
$\boldsymbol{x}^{\prime}$ | Flipped version of $\boldsymbol{x}$
$\boldsymbol{1}$ | Vector of all ones
$\mathbf{H}$ | Matrix (upper-case bold letter)
$\mathbf{H}^{\dagger}$ | Pseudo-inverse of $\mathbf{H}$
$\mathbf{I}$ | Identity matrix
$\mathrm{sgn}$ | Signum function
$\nabla$ | Gradient operator
$*$ | Convolution operator
$\langle\cdot,\cdot\rangle$ | Inner-product
$\oplus$ | Element-wise sum
$\odot$ | Element-wise product
$\mathbb{R}^{n}_{>t}$ | Set of vectors in $\mathbb{R}^{n}$ with entries greater than $t$

TABLE II: Sparsity-promoting regularizers used in this study and corresponding
proximal operators.

Name | Penalty Function | Proximal Operator
---|---|---
$\ell_{1}-\mathrm{norm}$ $(\lambda>0)$ | $g_{1}(x)=\lambda\left|x\right|$ | $\mathcal{P}_{g_{1}}(x)=\mathrm{sgn}(x)~{}\mathrm{max}\left(\left|x\right|-\lambda,0\right)$
MCP [21] $(\mu>0,\gamma>1)$ | $g_{2}(x)=\begin{cases}{\mu\left(\left|x\right|-\cfrac{{\left|x\right|}^{2}}{2\mu\gamma}\right),}\quad{\mathrm{for}\left|x\right|\leq\gamma\mu,}\\\ {\cfrac{{\mu}^{2}\gamma}{2},}\qquad\qquad\qquad\,{\mathrm{for}\left|x\right|\geq\gamma\mu.}\end{cases}$ | $\mathcal{P}_{g_{2}}(x)=\begin{cases}{0,}\qquad\qquad\qquad\qquad\qquad{\mathrm{for}\;\left|x\right|\leq\mu,}\\\ {\mathrm{sgn}\left(x\right)\cfrac{\gamma}{\gamma-1}\left(\left|x\right|-\mu\right),}\enspace\,{\mathrm{for}\;\mu<\left|x\right|\leq\gamma\mu,}\\\ {x,}\qquad\qquad\qquad\qquad\qquad{\mathrm{for}\;\left|x\right|>\gamma\mu.}\end{cases}$
SCAD [23] $(\nu>0,a>2)$ | $g_{3}(x)=\begin{cases}{\nu\left|x\right|,}\qquad\qquad\qquad\quad\,{\mathrm{for}\left|x\right|\leq\nu,}\\\ {\cfrac{{\left|x\right|}^{2}-2a\nu\left|x\right|+\nu^{2}}{2(1-a)},}\enspace\,\,{\mathrm{for}\nu<\left|x\right|\leq a\nu,}\\\ {\cfrac{(a+1){\nu^{2}}}{2},}\qquad\qquad\enspace\,\,{\mathrm{for}\left|x\right|>a\nu.}\end{cases}$ | $\mathcal{P}_{g_{3}}(x)=\begin{cases}{\mathrm{sgn}\left(x\right)\mathrm{max}\left(\left|x\right|-\nu,0\right),}\,{\mathrm{for}\;\left|x\right|\leq 2\nu,}\\\ {\cfrac{\left(a-1\right)x-\mathrm{sgn}\left(x\right)a\nu}{a-2},}\enspace\,{\mathrm{for}\;2\nu<\left|x\right|\leq a\nu,}\\\ {x,}\qquad\qquad\qquad\qquad\qquad{\mathrm{for}\;\left|x\right|>a\nu.}\end{cases}$

##  II Non-convex penalties and their weighted counterparts

The $\ell_{1}$ penalty (Fig. 2), which is the convex relaxation of the
$\ell_{0}$ penalty, is the most popular one for promoting sparsity while
solving sparse linear inverse problems, owing to its convexity-preserving
attribute, making the optimization problem computationally more tractable than
$\ell_{0}$ or non-convex regularization. However, $\ell_{1}$-regularization
results in a biased estimate. In particular, it over-penalizes large-valued
coefficients due to the constant shrinkage parameter. Non-convex penalties
such as the minimax concave penalty (MCP) (Fig. 3) and smoothly clipped
absolute deviation (SCAD) (Fig. 4) result in nearly unbiased estimates
compared to the $\ell_{1}$ penalty. Non-convex regularization overcomes the
limitation of $\ell_{1}$, but the resulting optimization problem is quite
challenging. Proximal operators are central to developing proximal gradient-
descent methods for optimizing non-convex, nonsmooth sparse linear inverse
problems. Further, the resulting optimization problem is solved via proximal
methods since proximal operators of MCP and SCAD are computationally simple
[18, 23]. The proximal operator $\mathcal{P}_{g}:\mathbb{R}\to\mathbb{R}$
associated with the penalty $g:\mathbb{R}\to\mathbb{R}_{\geq 0}$ is defined
as:

$\displaystyle\mathcal{P}_{g}(y,\theta_{g})=\arg\,\min_{x}\frac{1}{2}(x-y)^{2}+\lambda
g(x),$ (2)

where $\theta_{g}$ denotes the parameters of $g$ and $\lambda>0$.

Figure 2:  The $\ell_{1}$ penalty and its corresponding proximal operator, the
soft-thresholding operator. Figure 3:  The minimax concave penalty (MCP) and
its corresponding proximal operator, the firm-thresholding operator. Figure 4:
The smoothly clipped absolute deviation (SCAD) penalty and its corresponding
proximal operator.

The scaled $\ell_{1}$ penalty is defined as

$g_{1}(x;\lambda)=\lambda\left|x\right|,$ (3)

where the weight factor $\lambda$ allows for adjusting the threshold value.
The MCP $g_{2}:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ is defined as

$\displaystyle
g_{2}(x;\mu,\gamma)=\left\\{\begin{array}[]{lr}\displaystyle\mu|x|-\frac{|x|^{2}}{2\gamma},&|x|\leq\mu\gamma,\\\
\displaystyle\frac{\mu^{2}\gamma}{2},&|x|\geq\mu\gamma,\end{array}\right.$ (4)

where $\mu>0,\gamma>1$. The SCAD penalty
$g_{3}:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ is defined as

$\displaystyle
g_{3}(x;\nu,a)=\left\\{\begin{array}[]{lr}\nu|x|,&|x|\leq\nu,\\\
\displaystyle\frac{2\nu a|x|-|x|^{2}-\nu^{2}}{2(a-1)},&\nu<|x|\leq\nu a,\\\
\displaystyle\frac{\nu^{2}(a+1)}{2},&|x|>\nu a,\end{array}\right.$ (5)

where $\nu>0,a>2$. The penalties used in this study and their proximal
operators are shown in Table II.

We propose weighted counterparts of the above regularizers as atoms to
construct a learnable sparsity-prior. Further, the chosen penalty functions
have closed-form proximal operators. The weighted counterparts of the penalty
functions and their proximal operators are given in Figure 5.

  1. 1.

$g_{1}(\boldsymbol{x})=\|\boldsymbol{\lambda}\odot\boldsymbol{x}\|_{1}$, where
$\boldsymbol{\lambda}\in\mathbb{R}^{n}_{>0}$, and $\odot$ denotes element-wise
product. $\mathbb{R}^{n}_{>0}$ denotes the set of vectors in $\mathbb{R}^{n}$
containing positive entries.

  2. 2.

$g_{2}(\boldsymbol{x})$ is weighted MCP defined as
$g_{2}(\boldsymbol{x})=\sum_{j=1}^{n}g_{2}(x_{j}).$

  3. 3.

$g_{3}(\boldsymbol{x})$ is weighted SCAD, which is defined as
$g_{3}(\boldsymbol{x})=\sum_{j=1}^{n}g_{3}(x_{j}).$

Figure 5:  The penalty functions and proximal operators corresponding to Table
II.

##  III Proximal Average Network for Reflectivity Inversion — Nonuniform
Sparse Model (Type-$1$)

The sparse reflection coefficients $\boldsymbol{x}\in\mathbb{R}^{n}$ are
recovered from the observed noisy seismic trace
$\boldsymbol{y}\in\mathbb{R}^{n}$ via optimization of the following composite-
regularized cost:

$\begin{split}\arg\min_{\boldsymbol{x},\,\\{\omega_{i}\\}}\,\Big{\\{}\mathcal{J(\boldsymbol{x})}=\underbrace{\frac{1}{2}{\left\|\boldsymbol{h}*\boldsymbol{x}-\boldsymbol{y}\right\|}_{2}^{2}}_{f(\boldsymbol{x})}+\underbrace{\sum\limits_{i=1}^{m}\omega_{i}{g_{i}(\boldsymbol{x})}}_{g(\boldsymbol{x})}\Big{\\}}\\\
\text{subject
to}\,\,\sum\limits_{i=1}^{m}{\omega_{i}=1},0<\omega_{i}<1\;\forall{i},\end{split}$
(6)

where $m$ is the number of regularizers, and $\\{\omega_{i}\\}$ are the
weights assigned to the regularizers.

The composite objective function $\mathcal{J}$ in (6), regularizer $g$ and
data-fidelity term $f$ satisfy the following properties, which are essential
for minimizing $\mathcal{J}$,

  1. P1.

$f:\mathbb{R}^{n}\rightarrow\left(-\infty,\infty\right]$ is proper, closed,
and $L$-smooth, i.e., $\|\nabla f(\boldsymbol{x})-\nabla
f(\boldsymbol{y})\|_{2}\leq L\|\boldsymbol{x}-\boldsymbol{y}\|_{2}$.

  2. P2.

$g$ is lower semi-continuous.

  3. P3.

$\mathcal{J}$ is bounded from below, i.e., inf $\mathcal{J}>-\infty$.

Given a set of training samples
$\\{\boldsymbol{y}_{p},\boldsymbol{x}_{p}\\},\,\,\forall
p\in\\{1,2,3,\ldots,N\\}$, our objective is to design a deep-unrolled network
that solves the optimization problem stated in Eq. (6).

Our algorithm is referred to as the nonuniform proximal-averaged thresholding
algorithm (NuPATA-$1$), which relies on Majorization-Minimization (MM) [70]
and the proximal average strategy [26, 27, 43]. Further, we unfold the
NuPATA-$1$ iterations into a learnable network called nonuniform sparse
proximal average network (NuSPAN-$1$).

###  III-A Optimization Algorithm – Type-$1$

Due to P1, there exists $\eta<1/L$ such that $f(\boldsymbol{x})$ is upper
bounded locally by a quadratic expansion about
$\boldsymbol{x}=\boldsymbol{x}^{(k)}$ as follows:

$f(\boldsymbol{x})\leq Q\left(\boldsymbol{x},\boldsymbol{x}^{(k)}\right),$ (7)

where
$Q(\boldsymbol{x},\boldsymbol{x}^{(k)})=f(\boldsymbol{x}^{(k)})+\langle\nabla
f(\boldsymbol{x}^{(k)}),\boldsymbol{x}-\boldsymbol{x}^{k}\rangle+\frac{1}{2\eta}\|\boldsymbol{x}-\boldsymbol{x}^{(k)}\|_{2}^{2}$.
The majorizer to $\mathcal{J}$ at $\boldsymbol{x}^{(k)}$ is

$\mathcal{J}(\boldsymbol{x})\leq\underbrace{Q(\boldsymbol{x},\boldsymbol{x}^{(k)})+g(\boldsymbol{x})}_{H(\boldsymbol{x},\boldsymbol{x}^{(k)})}.$
(8)

The update for $\boldsymbol{x}$ at the $(k+1)^{\text{th}}$ iteration is given
by

$\displaystyle\boldsymbol{x}^{(k+1)}$
$\displaystyle=\arg\,\min_{\boldsymbol{x}}H(\boldsymbol{x},\boldsymbol{x}^{(k)}),$
$\displaystyle=\arg\,\min_{\boldsymbol{x}}\frac{1}{2\eta}\left\|\boldsymbol{x}-\left(\boldsymbol{x}^{(k)}-\eta\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{(k)})\right)\right\|_{2}^{2}$
$\displaystyle\hskip 85.35826pt+g(\boldsymbol{x}).$ (9)

The above problem does not have a closed-form solution. Hence, we consider an
approximate variant of the problem corresponding to Eq. (9) based on the
proximal average strategy:

$\displaystyle\boldsymbol{x}^{(k+1)}$
$\displaystyle=\arg\min_{\boldsymbol{x}}\sum\limits_{i=1}^{m}\frac{\omega_{i}}{2\eta}{\left\|\boldsymbol{x}-\left(\boldsymbol{x}^{(k)}-\eta\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{(k)})\right)\right\|}_{2}^{2}$
$\displaystyle\hskip
56.9055pt+\sum\limits_{i=1}^{m}{\omega_{i}g_{i}(\boldsymbol{x})},$
$\displaystyle=\sum\limits_{i=1}^{m}\omega_{i}\mathcal{P}_{g_{i}}(\mathbf{u}^{(k)}),$
(10)

where
$\mathbf{u}^{(k)}=\boldsymbol{x}^{(k)}-\eta\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{(k)})$,
$\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{(k)})=-\boldsymbol{h}^{\prime}*\left(\boldsymbol{y}-\boldsymbol{h}*\boldsymbol{x}^{(k)}\right)$,
where $\boldsymbol{h}^{\prime}$ is the flipped version of $\boldsymbol{h}$,
and $\mathcal{P}_{g_{i}}$ represent the proximal operators corresponding to
$g_{i}$ (cf. Table II). The optimization procedure is listed in Algorithm 1.

Input: $\boldsymbol{y}$, $\boldsymbol{h}$,
$L={\left\|\boldsymbol{h}\right\|}_{2}^{2},k_{max}$,
$\\{\omega_{i}:0<\omega_{i}<1,\sum_{i=1}^{m}\omega_{i}=1\\}$,
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$,
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$

Initialize: $\boldsymbol{x}^{(0)}=\boldsymbol{0}$

while  _$k <k_{max}$_ do

$\boldsymbol{z}^{(k+1)}=\boldsymbol{x}^{(k)}+(1/2L)\left(\boldsymbol{h}^{\prime}*(\boldsymbol{y}-\boldsymbol{h}*{\boldsymbol{x}^{(k)}})\right)$

$\boldsymbol{x}^{(k+1)}=\sum\limits_{i=1}^{m}\omega_{i}\mathcal{P}_{g_{i}}\left(\boldsymbol{z}^{(k+1)}\right)$

$k=k+1$

Output: $\hat{\boldsymbol{x}}\leftarrow\boldsymbol{x}^{(k_{max})}$

Algorithm 1 Nonuniform Sparse Proximal-Averaged Thresholding Algorithm –
Type-$1$ (NuPATA-$1$)

Input: $\boldsymbol{y}$, $L={\left\|\mathbf{H}\right\|}_{2}^{2},k_{max}$,
$\\{\omega_{i}:0<\omega_{i}<1,\sum_{i=1}^{m}\omega_{i}=1\\}$,
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$,
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$

Initialize:
$\mathbf{W},\mathbf{S},\boldsymbol{x}^{(0)}=\sum\limits_{i=1}^{m}\omega_{i}\mathcal{P}_{g_{i}}(\mathbf{W}{\boldsymbol{y}})$

while  _$k <k_{max}$_ do

$\mathbf{c}^{(k+1)}=\mathbf{W}{\boldsymbol{y}}+\mathbf{S}{\boldsymbol{x}^{(k)}}$

$\boldsymbol{x}^{(k+1)}=\sum\limits_{i=1}^{m}\omega_{i}\mathcal{P}_{g_{i}}(\mathbf{c}^{(k+1)})$

$k=k+1$

$\hat{\boldsymbol{x}}_{\mathcal{S}(\hat{\boldsymbol{x}})}=\mathbf{H}_{\mathcal{S}(\hat{\boldsymbol{x}})}^{\dagger}\boldsymbol{y}$

Output: $\hat{\boldsymbol{x}}\leftarrow\boldsymbol{x}^{(k_{max})}$

Algorithm 2 Nonuniform Sparse Proximal Average Network – Type-$1$ (NuSPAN-$1$)

###  III-B Nonuniform Sparse Proximal Average Network – Type-$1$ (NuSPAN-$1$)

The update in (10) involves convolutions followed by nonlinear activation, and
can therefore be represented as a layer in a neural network. We unfold the
iterations of Algorithm 1 to obtain the deep-unrolled architecture, namely,
nonuniform sparse proximal average network (NuSPAN-$1$) for solving the
reflectivity inversion problem in Eq. (6). The structure for each layer in
NuSPAN-$1$ is given by

$\boldsymbol{x}^{(k+1)}=\sum_{i=1}^{m}\omega_{i}\mathcal{P}_{g_{i}}(\mathbf{W}\boldsymbol{y}+\mathbf{S}\boldsymbol{x}^{(k)}),$
(11)

where $\mathbf{W}=(1/L)\mathop{\boldsymbol{h}^{\prime}}$ and
$\mathbf{S}=\mathbf{I}-(1/L)\mathop{\boldsymbol{h}^{\prime}*\boldsymbol{h}}$
[28] are initialized as Toeplitz matrices. In the learning stage, they are
dense and unstructured. Given training data that consists of a large number of
independent and identically distributed examples
$\\{\boldsymbol{x}_{p},\boldsymbol{y}_{p}\\}_{p=1}^{N}$ and NuSPAN-$1$ with a
fixed number of layers, we optimize the smooth $\ell_{1}$ cost computed
between the true reflectivity $\boldsymbol{x}$ and the prediction
$\hat{\boldsymbol{x}}(\boldsymbol{\theta})$:

$\displaystyle\min_{\boldsymbol{\theta}}\frac{1}{N}\sum\limits_{i=1}^{N}\left(\beta\left\|{\boldsymbol{x}-\hat{\boldsymbol{x}}(\boldsymbol{\theta})}\right\|_{1}+(1-\beta)\left\|{\boldsymbol{x}-\hat{\boldsymbol{x}}(\boldsymbol{\theta})}\right\|_{2}^{2}\right),$

where $0<\beta<1$. We set $\beta=1$; it is not a trainable parameter. The
parameters $\boldsymbol{\theta}$ that need to be learned are the matrices
$\mathbf{W}$ and $\mathbf{S}$, the parameters of $g_{i}$:
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$ and
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$,
and the weights $\\{\omega_{i}:0<\omega_{i}<1,\sum_{i=1}^{m}\omega_{i}=1\\}$.
The NuSPAN-$1$ algorithm for $k_{max}$ layers is listed in Algorithm 2.

We enhance the amplitude recovery of NuSPAN-$1$ (and NuSPAN-$2$ in Section IV)
during the inference phase by re-estimating the amplitudes over the supports
given by NuSPAN-$1$ as:
$\hat{\boldsymbol{x}}_{\mathcal{S}(\hat{\boldsymbol{x}})}=\mathbf{H}_{\mathcal{S}(\hat{\boldsymbol{x}})}^{\dagger}\boldsymbol{y}$,
where $\hat{\boldsymbol{x}}$ is the sparse vector estimated by NuSPAN-$1$,
$\mathcal{S}(\hat{\boldsymbol{x}})$ is the support, i.e., non zero locations
of $\hat{\boldsymbol{x}}$ , and
$\mathbf{H}_{\mathcal{S}(\hat{\boldsymbol{x}})}^{\dagger}$ is the pseudo-
inverse of the Toeplitz matrix $\mathbf{H}$ of the kernel $\boldsymbol{h}$
over $\mathcal{S}(\hat{\boldsymbol{x}})$.

Figure 6:  Structure of layers in NuSPAN-$2$. The matrices $\mathbf{W}$,
$\mathbf{S}$, and
$\left[\boldsymbol{\omega}_{1}~{}\boldsymbol{\omega}_{2}~{}\boldsymbol{\omega}_{3}\right]$,
and $\boldsymbol{\Phi}_{i}$, which are the parameters of $g_{i}$:
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$ and
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$
(corresponding to the penalties in Table II), are shared across layers of
NuSPAN-$2$. The projection operator refers to the constraints imposed
according to Table II through the Pytorch [71] clamp function. The operators
$\oplus$ and $\odot$ represent sum and product, all considered element-wise.

Input: $\boldsymbol{y}$, $\boldsymbol{h}$,
$L={\left\|\boldsymbol{h}\right\|}_{2}^{2},k_{max}$,
$\\{\boldsymbol{\omega}_{i}:\boldsymbol{0}<\omega_{i}<\boldsymbol{1},\sum_{i=1}^{m}\boldsymbol{\omega}_{i}=\boldsymbol{1}\\}$,
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$,
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$

Initialize: $\boldsymbol{x}^{(0)}=\boldsymbol{0}$

while  _$k <k_{max}$_ do

$\boldsymbol{z}^{(k+1)}=\boldsymbol{x}^{(k)}+(1/2L)\left(\boldsymbol{h}^{\prime}*(\boldsymbol{y}-\boldsymbol{h}*{\boldsymbol{x}^{(k)}})\right)$

$\boldsymbol{x}^{(k+1)}=\sum\limits_{i=1}^{m}\boldsymbol{\omega}_{i}\odot\mathcal{P}_{g_{i}}\left(\boldsymbol{z}^{(k+1)}\right)$

$k=k+1$

Output: $\hat{\boldsymbol{x}}\leftarrow\boldsymbol{x}^{(k_{max})}$

Algorithm 3 Nonuniform Sparse Proximal-Averaged Thresholding Algorithm –
Type-$2$ (NuPATA-$2$)

Input: $\boldsymbol{y}$, $L={\left\|\mathbf{H}\right\|}_{2}^{2},k_{max}$,
$\\{\boldsymbol{\omega}_{i}:\boldsymbol{0}<\omega_{i}<\boldsymbol{1},\sum_{i=1}^{m}\boldsymbol{\omega}_{i}=\boldsymbol{1}\\}$,
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$,
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$

Initialize:
$\mathbf{W},\mathbf{S},\boldsymbol{x}^{(0)}=\sum\limits_{i=1}^{m}\boldsymbol{\omega}_{i}\odot\mathcal{P}_{g_{i}}(\mathbf{W}{\boldsymbol{y}})$

while  _$k\leq k_{max}$_ do

$\mathbf{c}^{(k+1)}=\mathbf{W}{\boldsymbol{y}}+\mathbf{S}{\boldsymbol{x}^{(k)}}$

$\boldsymbol{x}^{(k+1)}=\sum\limits_{i=1}^{m}\boldsymbol{\omega}_{i}\odot\mathcal{P}_{g_{i}}(\mathbf{c}^{(k+1)})$

$k=k+1$

$\hat{\boldsymbol{x}}_{\mathcal{S}(\hat{\boldsymbol{x}})}=\mathbf{H}_{\mathcal{S}(\hat{\boldsymbol{x}})}^{\dagger}\boldsymbol{y}$

Output: $\hat{\boldsymbol{x}}\leftarrow\boldsymbol{x}^{(k_{max})}$

Algorithm 4 Nonuniform Sparse Proximal Average Network – Type-$2$ (NuSPAN-$2$)

##  IV Proximal Average Network for Reflectivity Inversion — Nonuniform Sparse
Model (Type-$2$)

The generalized variant of the composite-regularization problem given in Eq.
(6) is given as follows:

$\displaystyle\arg\min_{\boldsymbol{x},\,\\{\boldsymbol{\omega}_{i}\\}}$
$\displaystyle\,{\frac{1}{2}{\left\|\boldsymbol{h}*\boldsymbol{x}-\boldsymbol{y}\right\|}_{2}^{2}}+\sum\limits_{i=1}^{m}\langle\boldsymbol{\omega}_{i},q_{i}(\boldsymbol{x})\rangle$
$\displaystyle\text{s.t.}\sum\limits_{i=1}^{m}{\boldsymbol{\omega}_{i}=\boldsymbol{1}},\boldsymbol{0}<\boldsymbol{\omega}_{i}<\boldsymbol{1},\;\forall{i},$
(12)

where
$q_{i}(\boldsymbol{x})=[g_{i}(x_{1}),g_{i}(x_{2}),\ldots,g_{i}(x_{n})]^{\textsc{T}}$,
$g_{i}$ is the scalar penalty provided in Table II,
$\langle\cdot,\cdot\rangle$ denotes inner-product, $\boldsymbol{0}$ is the
null vector, $\boldsymbol{1}$ is the vector of all ones, and
$\boldsymbol{\omega}_{i}=[\omega_{i1},\omega_{i2},\ldots,\omega_{in}]^{\textsc{T}}$.

###  IV-A Optimization algorithm – Type-$2$

We can follow the strategy mentioned in Section III-A for optimizing the cost
given in Eq. (12). The update for $\boldsymbol{x}$ at $(k+1)^{\text{th}}$
iteration is given by

$\displaystyle\boldsymbol{x}^{(k+1)}$
$\displaystyle=\arg\,\min_{\boldsymbol{x}}H\left(\boldsymbol{x},\boldsymbol{x}^{(k)}\right),$
$\displaystyle=\arg\,\min_{\boldsymbol{x}}\frac{1}{2\eta}\left\|\boldsymbol{x}-\left(\boldsymbol{x}^{(k)}-\eta\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{(k)})\right)\right\|_{2}^{2}$
$\displaystyle\hskip
85.35826pt+\sum\limits_{i=1}^{m}\langle\boldsymbol{\omega}_{i},q_{i}(\boldsymbol{x})\rangle.$
(13)

Since the above problem is complex and non-convex, we solve the approximate
variant of the problem corresponding to Eq. (13) based on the proximal average
strategy [27].

$\displaystyle\boldsymbol{x}^{(k+1)}$
$\displaystyle=\arg\min_{\boldsymbol{x}}\sum\limits_{i=1}^{m}\frac{\omega_{i}}{2\eta}{\left\|\boldsymbol{x}-\left(\boldsymbol{x}^{(k)}-\eta\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{(k)})\right)\right\|}_{2}^{2}$
$\displaystyle\hskip
56.9055pt+\sum\limits_{i=1}^{m}\langle\boldsymbol{\omega}_{i},q_{i}(\boldsymbol{x})\rangle,$
$\displaystyle=\sum\limits_{i=1}^{m}\boldsymbol{\omega}_{i}\odot\mathcal{P}_{g_{i}}(\mathbf{u}^{(k)}).$
(14)

Algorithm (3) illustrates the NuPATA-$2$ formulation for solving the
optimization problem (12) using the proximal average strategy. As in the
previous section, we unfold the update in (14) into a layer of a neural
network.

###  IV-B Nonuniform Sparse Proximal Average Network – Type-$2$ (NuSPAN-$2$)

Representing the NuPATA-$2$ update in (14) as a layer in a neural network, we
obtain the nonuniform sparse proximal average network (NuSPAN-$2$). The
structure for each layer in NuSPAN-$2$ is given by

$\boldsymbol{x}^{(k+1)}=\sum_{i=1}^{m}\boldsymbol{\omega}_{i}\odot\mathcal{P}_{g_{i}}(\mathbf{W}\boldsymbol{y}+\mathbf{S}\boldsymbol{x}^{(k)}),$
(15)

where $\mathbf{W}$ and $\mathbf{S}$ definitions are provided in Section III-B.
Similar to NuSPAN-$1$, we optimize the learnable parameters of NuSPAN-$2$
subject to the constraints given in Eq. (12), by minimizing the smooth
$\ell_{1}$ cost defined between the true reflectivity $\boldsymbol{x}$ and
prediction $\hat{\boldsymbol{x}}(\boldsymbol{\theta})$, where
$\boldsymbol{\theta}$ denotes the learnable parameters such as the matrices
$\mathbf{W}$ and $\mathbf{S}$, the parameters of $g_{i}(\cdot)\,\forall i$:
$\\{\boldsymbol{\lambda$,$\mu$,$\nu}\\}\in\mathbb{R}^{n}_{>0}$ and
$\\{\boldsymbol{\gamma}\in\mathbb{R}^{n}_{>1},\boldsymbol{a}\in\mathbb{R}^{n}_{>2}\\}$,
and the weights
$\\{\boldsymbol{\omega}_{i}:\boldsymbol{0}<\omega_{i}<\boldsymbol{1},\sum_{i=1}^{m}\boldsymbol{\omega}_{i}=\boldsymbol{1}\\}$.
Algorithm 4 gives the NuSPAN-$2$ algorithm for $k_{max}$ layers, and Figure 6
gives the structure of layers in NuSPAN-$2$.

##  V Experimental Results

We demonstrate the efficacy of the proposed networks, namely, NuSPAN-$1$ and
NuSPAN-$2$ on both synthetic and simulated datasets as well as on real data in
comparison with the benchmark techniques such as BPI [10, 9], FISTA [11, 12],
and SBL-EM [68, 16]. The performance of the proposed approaches is quantified
based on objective metrics computed between ground-truth sparse vector
$\boldsymbol{x}$ and the predicted sparse vector $\hat{\boldsymbol{x}}$ listed
in the following section.

###  V-A Objective Metrics

We evaluate the performance in terms of amplitude and support recovery
metrics, which are crucial to estimate the amplitudes and locations of
reflection coefficients.

  * •

Correlation Coefficient (CC) [72]:

$\displaystyle\text{CC}=\frac{\langle{\boldsymbol{x}},\hat{{\boldsymbol{x}}}\rangle-\langle{\boldsymbol{x}},{{\boldsymbol{1}}}\rangle\langle\hat{{\boldsymbol{x}}},{{\boldsymbol{1}}}\rangle}{\sqrt{\left(\left\|{{\boldsymbol{x}}}\right\|_{2}^{2}-{\langle{\boldsymbol{x}},{{\boldsymbol{1}}}\rangle}^{2}\right)\left(\left\|\hat{{\boldsymbol{x}}}\right\|_{2}^{2}-{\langle\hat{{\boldsymbol{x}}},{{\boldsymbol{1}}}\rangle}^{2}\right)}},$

  * •

Relative Reconstruction Error (RRE) and Signal-to-Reconstruction Error Ratio
(SRER) (in dB)

$\displaystyle\text{RRE}=\frac{{\left\|\hat{\boldsymbol{x}}-\boldsymbol{x}\right\|}_{2}^{2}}{{\left\|\boldsymbol{x}\right\|}_{2}^{2}},\text{SRER}=10\mathop{\log\limits_{10}\left(\frac{{\left\|{\boldsymbol{x}}\right\|_{2}^{2}}}{{\left\|{\hat{\boldsymbol{x}}-{\boldsymbol{x}}}\right\|_{2}^{2}}}\right)}.$

  * •

Probability of Error in Support (PES) is given by

$\displaystyle\text{PES}=\frac{1}{t}\mathop{\sum\limits_{i=1}^{t}\left(\frac{{{\text{max}}(\left|{\mathcal{S}({\hat{\boldsymbol{x}_{i}}})}\right|,\left|{\mathcal{S}({\boldsymbol{x}_{i}})}\right|)-\left|{\mathcal{S}({\hat{\boldsymbol{x}_{i}}})\mathop{\cap}{\mathcal{S}}({\boldsymbol{x}_{i}})}\right|}}{{{\text{max}}(\left|{\mathcal{S}({\hat{\boldsymbol{x}_{i}}})}\right|,\left|{\mathcal{S}({\boldsymbol{x}_{i}})}\right|)}}\right)},$

where $|\cdot|$ denotes the cardinality, and $\mathcal{S}(\cdot)$ denotes the
support of the argument.

Figure 7:  Sample results for a synthetic $1$-D seismic trace. (a) True
seismic trace; (b) noisy seismic trace; (c) true reflectivity; (d)-(h)
Recovered reflectivity (blue circles) compared with true reflectivity (red
crosses). The proposed techniques, NuSPAN-$1$ (g) and NuSPAN-$2$ (h),
distinguish between closely-spaced spikes around $150$ ms; the benchmark
techniques BPI, FISTA, and SBL-EM (d)-(f) predict a single reflector instead.

###  V-B Training Phase

The synthetic training data that appropriately represents observed seismic
data is generated as recommended by [53]. We generate synthetic training data
of size $5\times 10^{5}$ seismic traces, each consisting of $300$ samples,
obtained by convolving $1$-D reflectivity profiles with a Ricker wavelet. Each
reflectivity profile contains $200$ samples with amplitudes ranging from
$-1.0$ to $1.0$, padded with zeros before convolution with the wavelet to have
the same length as the seismic trace [4]. The sparsity factor, which is the
ratio of the number of non-zero elements to the total number of elements, is
set to $0.05$. We also report results for sparsity factors $0.10$, $0.15$, and
$0.20$ in the supplementary document. The reflector/spike locations are chosen
uniformly at random, without any minimum spacing constraint between any two
spikes. Reflectivity values are then assigned to these locations, picked
randomly from the amplitude range of $-1.0$ to $1.0$. The amplitude increment,
sampling interval, wavelet frequency, and the initial hyperparameters vary
depending on the dataset. The optimum number of layers also varies with the
dataset, and we consider $10$, $15$, or $20$ layer-models in our experiments.
Our models are trained with a batch size of $200$, use the ADAM optimizer
[73], and with the learning rate set to $1\times 10^{-3}$, and input
measurement SNR set to $10$ dB — to ensure robustness against noisy testing
data [53]. We consider models trained on $1$-D data to operate trace-by-trace
on all the datasets: synthetic $1$-D data, synthetic $2$-D wedge models [74],
simulated Marmousi$2$ data [75], and real data.

###  V-C Testing Phase – Synthetic $1$-D & $2$-D Data

####  V-C1 Synthetic $1$-D Traces

We validate the performance of the proposed models over $1000$ realizations of
synthetic $1$-D traces, with a $30$ Hz Ricker wavelet, $1$ ms sampling
interval, $0.2$ amplitude increments, and the minimum spacing between
reflection coefficients (spikes) $1$ ms. Comparisons across several methods
are reported in Figure 7 and Table III. Figure 7 shows the results for a
sample synthetic $1$-D trace out of the $1000$ test realizations for which the
results are reported in Table III. NuSPAN-$1$ and NuSPAN-$2$ resolve closely-
spaced reflection coefficients, whereas the benchmark techniques (BPI, FISTA,
and SBL-EM) predict a single reflector. From Table III, observe that
NuSPAN-$1$ and NuSPAN-$2$ (Nonuniform Sparse Proximal Average Network, Type
$1$ and $2$) recover amplitudes with higher accuracy than the benchmark
methods, and NuSPAN-$1$ exhibits superior support recovery. The low
computation time of the proposed models is crucial for processing the large
volume of data in reflection seismic processing.

TABLE III: Performance evaluation in terms of objective metrics averaged over
$1000$ test realizations of synthetic $1$-D traces. The proposed NuSPAN-$1$
and NuSPAN-$2$ show superior amplitude recovery in terms of CC, RRE, and SRER
and significantly lower computation time, while NuSPAN-$1$ outperforms in
support recovery in terms of PES. The best performance is highlighted in
boldface. The second best is underlined.

Method | CC | RRE | SRER | PES | Time ($\mathrm{s}$)
---|---|---|---|---|---
BPI | ${0.5499}$ | ${0.7290}$ | ${1.8687}$ | ${0.9704}$ | ${503.8519}$
FISTA | ${0.5473}$ | ${0.7203}$ | ${1.8391}$ | $\underline{0.8112}$ | ${33.2539}$
SBL-EM | ${0.5501}$ | ${0.7682}$ | ${1.8276}$ | ${0.9704}$ | ${1076.1833}$
NuSPAN-$1$ | $\underline{0.5979}$ | $\underline{0.6354}$ | $\underline{2.2038}$ | $\mathbf{0.7104}$ | $\mathbf{0.1778}$
NuSPAN-$2$ | $\mathbf{0.6050}$ | $\mathbf{0.6274}$ | $\mathbf{2.2508}$ | ${0.9563}$ | $\underline{0.1870}$

####  V-C2 Synthetic $2$-D Wedge Models

Synthetic $2$-D wedge models are considered to evaluate resolving capability
on thin beds [74]. A wedge model typically consists of two interfaces, one
horizontal and another inclined, with the polarity (N: Negative, P: Positive)
of reflection coefficients same or opposite, giving four possible types of
wedge models. Here, we report results from an odd wedge model (NP) with
negative polarity on the upper horizontal interface (N) and positive on the
lower inclined interface (P), and an even (NN) wedge model, while those for
the remaining two variants are given in the supplementary. In our experimental
setup, each model consists of $26$ traces, with separation between reflectors
of the two interfaces increasing from $0$ ms to $50$ ms in $2$ ms increments.
The amplitudes of the reflectors are fixed as $\pm~{}0.5$ based on the
polarity.

Results for the NP odd wedge model are given in Table IV and Figure 8.
NuSPAN-$1$ outperforms other benchmark techniques in terms of the objective
metrics, namely, CC and RRE, and shows superior support recovery, measured in
terms of PES. Figure 8 highlights the fact that BPI, FISTA, and SBL-EM fail to
resolve the locations of closely-spaced reflectors, which is evident from the
divergence observed below the tuning thickness of $13$ ms (wedge thickness
between $5$-$6$ m) [76].

Table V and Figure 9 give the results for the NN even wedge model. For the
even wedge model, the baselines perform marginally better than NuSPAN-$1$ in
terms of amplitude recovery, measured through CC, RRE and SRER, but NuSPAN-$1$
outperforms in terms of support recovery (Table V). Figure 9 shows that the
baselines and the proposed networks recover the reflectivity profile well for
the even wedge model.

Although NuSPAN-$1$ outperforms the baselines in all the objective metrics
except SRER, we observe a drop in the performance of NuSPAN-$1$ and NuSPAN-$2$
in the synthetic $2$-D case over the $1$-D case in Section V-C1, especially in
terms of SRER. We hypothesize that this could be attributed to the mismatch
between the training and testing conditions. Further, we observe a disparity
between the performance of the proposed networks in the case of the odd (NP)
vs. even (NN) wedge model. This discrepancy could be attributed to the
suboptimal amplitude recovery of reflection coefficients in the region where
destructive interference is observed between the two interfaces of the wedge
models (as highilghted in Figure 9). Nonetheless, the support recovery
evaluated in terms of PES remains comparable in the case of both wedge models
(compare Tables IV and V). However, these aspects need further investigation.

TABLE IV: Metrics for a synthetic $2$-D odd (NP) wedge model. NuSPAN-$1$
outperforms in amplitude recovery in terms of CC and RRE and support recovery
in terms of PES. The best performance is highlighted in boldface. The second
best is underlined.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | ${0.7865}$ | ${0.2919}$ | $\mathbf{12.3197}$ | ${0.9933}$
FISTA | ${0.7786}$ | ${0.3008}$ | $\underline{12.0683}$ | $\underline{0.7090}$
SBL-EM | $\underline{0.8064}$ | $\underline{0.2693}$ | ${10.1708}$ | ${0.9933}$
NuSPAN-$1$ | $\mathbf{0.8736}$ | $\mathbf{0.2553}$ | ${8.4383}$ | $\mathbf{0.5994}$
NuSPAN-$2$ | ${0.7858}$ | ${0.3441}$ | ${5.3164}$ | ${0.9896}$

Figure 8:  Results for a synthetic $2$-D odd (NP) wedge model. True (a)
seismic traces and (c) reflectivity. (b) Noisy seismic traces with input SNR
$10$ dB. (d)-(h) Recovered reflectivity profiles show that NuSPAN-$1$ (g) and
NuSPAN-$2$ (h) resolve reflectors $<5$ m thickness, whereas the baselines
(d)-(f) fail to do so, evident from the diverging interfaces highlighted by
the rectangle in black.

###  V-D Testing Phase – Marmousi$2$ Model

The Marmousi$2$ model [75] is widely used in reflection seismology to
calibrate algorithms in structurally complex settings. The model (width
$\times$ depth: $17$ km $\times~{}3.5$ km), an expanded version of the
original Marmousi model (width $\times$ depth: $9.2$ km $\times~{}3$ km) [77],
has a $2$ ms sampling interval, with traces at an interval of $12.5$ m. We
obtained the reflectivity profile from the P-wave velocity and density models
(1) and convolved them with a $30$ Hz Ricker wavelet to generate the
measurement.

Figure 10 shows the result for a region of the model with a gas-charged sand
channel [75]. NuSPAN-$1$ and NuSPAN-$2$ preserve the lateral continuity
better, evident from an observation of the insets in Figure 10. BPI and FISTA
introduce false interfaces due to the interference of multiple events at the
ends of the channel. The testing times mentioned in Table III for synthetic
data, when computed for the Marmousi$2$ model, further highlight the advantage
of the proposed approach for seismic processing (Table VI).

TABLE V: Results for a synthetic $2$-D odd (NN) wedge model. NuSPAN-$1$ offers
superior support recovery in terms of PES. The benchmark techniques (BPI,
FISTA, and SBL-EM) show higher amplitude recovery accuracy quantified in terms
of CC, RRE, and SRER.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | ${0.8504}$ | $\underline{0.2831}$ | $\mathbf{12.4768}$ | ${0.9933}$
FISTA | $\mathbf{0.8791}$ | $\mathbf{0.2524}$ | $\underline{11.7202}$ | $\underline{0.6681}$
SBL-EM | $\underline{0.8519}$ | ${0.3280}$ | ${10.7475}$ | ${0.9933}$
NuSPAN-$1$ | ${0.8208}$ | ${0.3963}$ | ${9.0663}$ | $\mathbf{0.5885}$
NuSPAN-$2$ | ${0.7591}$ | ${0.4605}$ | ${5.5752}$ | ${0.9896}$

Figure 9:  Results for a synthetic $2$-D even (NN) wedge model. True (a)
seismic traces and (c) reflectivity; (b) Noisy seismic traces with input SNR
$10$ dB; (d)-(h) Recovered reflectivity. The underperformance of NuSPAN-$1$
and NuSPAN-$2$ reported in Table V could be due to the poor amplitude
recovery, as a result of the destructive interference between the two
interfaces, in the region highlighted in black. TABLE VI: Metrics for a
portion of the Marmousi$2$ model. NuSPAN-$2$ outperforms in terms of CC, while
NuSPAN-$1$ shows superior support recovery in terms of PES. The reported time
is for the complete model, showing significantly lower computation time for
NuSPAN-$1$ and NuSPAN-$2$. The best performance is highlighted in boldface.
The second best is underlined.

Method | CC | RRE | SRER | PES | Time ($\mathrm{h}$)
---|---|---|---|---|---
BPI | ${0.9473}$ | ${0.0875}$ | ${14.8024}$ | ${0.9724}$ | ${16.2671}$
FISTA | ${0.9407}$ | ${0.1017}$ | ${13.8579}$ | $\underline{0.7146}$ | ${5.5636}$
SBL-EM | $\underline{0.9549}$ | $\mathbf{0.0684}$ | $\mathbf{18.1890}$ | ${0.9724}$ | ${101.2659}$
NuSPAN-$1$ | ${0.9376}$ | ${0.1032}$ | ${14.2592}$ | $\mathbf{0.3693}$ | $\mathbf{0.0737}$
NuSPAN-$2$ | $\mathbf{0.9591}$ | $\underline{0.0715}$ | $\underline{15.0973}$ | ${0.9662}$ | $\underline{0.1024}$

Figure 10:  P-wave velocity profile for the (a) Marmousi$2$ model and (b) the
portion corresponding to Table VI. (c) True reflectivity (ground truth), and
(d)-(h) predicted reflectivity. The inset plots are zoomed-in portions of the
selected area. NuSPAN-$1$ and NuSPAN-$2$ preserve lateral continuity better.

###  V-E Testing Phase – Real Data

We validate on real data from the field, a $3$-D volume from the Penobscot
$3$-D survey off the coast of Nova Scotia, Canada [78]. We pick a portion of
the $3$-D volume, with $201$ Inlines (from inline $1150$-$1350$) and $121$
Xlines (between $1040$-$1160$), chosen such that the region includes two wells
(wells L-$30$ and B-$41$) [79]. The sampling interval for the dataset is $4$
ms, and the region chosen includes $800$ samples between $0$-$3196$ ms along
the time/depth axis. A $25$ Hz Ricker wavelet fits the data well, also
observed by [79].

The recovered reflectivity profiles are shown in Figure 11, along with the
reflectivity profiles calculated from the sonic logs of well L-$30$ overlaid
in black. Figure 11 shows that inverted reflectivity for BPI and FISTA is
smooth and missing details, with relatively poor amplitude recovery (following
convention, red indicates positive reflectivity, and blue, negative). SBL-EM
results provide a more detailed image for characterization by recovering the
sparse reflectivity profiles, but the method fails to remove the noise. The
NuSPAN-$1$ and NuSPAN-$2$ recovered reflectivity profiles show better
amplitude recovery while removing the noise from the observations. Further,
the NuSPAN variants also preserve lateral continuity, especially for closely
spaced interfaces, as seen around $1.1$ s on the time axis in Figure 11.

Figure 11:  (a) Observed seismic data, and (b)-(f) predicted reflectivity
profiles for the inset marked in (a) for Xline $1155$ of the Penobscot $3$-D
survey [78]. The overlaid waveforms in black show the recorded (a) well
seismic and (b)-(f) reflectivity profiles, respectively, at well L-30 [79].
NuSPAN-$1$ and NuSPAN-$2$ show superior amplitude recovery while also removing
noise from the observations.

##  VI Conclusions

We considered the problem of seismic reflectivity inversion based on
systematically learned non-convex sparse-prior. We proposed a nonuniform
sparse model based on composite regularization and solved it based on a deep-
unrolled architecture for prior learning given training data. Using the
proposed framework, we solved the problem of seismic reflectivity inversion,
where the challenge lies in improving the resolution to characterize the
subsurface. The proposed techniques outperform benchmark reconstructions given
by state-of-the-art techniques on synthetic, simulated, and real datasets in
terms of objective metrics.

## Acknowledgments

This work is supported by the Ministry of Earth Sciences, Government of India;
Centre of Excellence in Advanced Mechanics of Materials, Indian Institute of
Science (IISc), Bangalore; and Science and Engineering Research Board (SERB),
India.

## References

  * [1] P. M. Shearer, _Introduction to Seismology_. Cambridge Univ. Press, 2009.
  * [2] D. W. Oldenburg, T. Scheuer, and S. Levy, “Recovery of the acoustic impedance from reflection seismograms,” _Geophysics_ , vol. 48, no. 10, pp. 1318–1337, 1983.
  * [3] Ö. Yilmaz, _Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data_. Society of Exploration Geophysicists, 2001.
  * [4] B. Russell, “Machine learning and geophysical inversion — A numerical study,” _The Leading Edge_ , vol. 38, no. 7, pp. 512–519, 2019.
  * [5] H. L. Taylor, S. C. Banks, and J. F. McCoy, “Deconvolution with the $\ell_{1}$ norm,” _Geophysics_ , vol. 44, no. 1, pp. 39–52, 1979.
  * [6] A. J. Berkhout, “Least-squares inverse filtering and wavelet deconvolution,” _Geophysics_ , vol. 42, no. 7, pp. 1369–1383, 1977.
  * [7] H. W. J. Debeye and P. Van Riel, “$\ell_{p}$-norm deconvolution,” _Geophysical Prospecting_ , vol. 38, no. 4, pp. 381–403, 1990.
  * [8] A. Tarantola, _Inverse Problem Theory and Methods for Model Parameter Estimation_. Society for Industrial and Applied Mathematics, 2005.
  * [9] R. Zhang and J. Castagna, “Seismic sparse-layer reflectivity inversion using basis pursuit decomposition,” _Geophysics_ , vol. 76, no. 6, pp. R147–R158, 2011.
  * [10] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” _SIAM Review_ , vol. 43, no. 1, pp. 129–159, 2001.
  * [11] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” _SIAM J. Imaging Sciences_ , vol. 2, no. 1, pp. 183–202, 2009.
  * [12] D. O. Pérez, D. R. Velis, and M. D. Sacchi, “Inversion of prestack seismic data using FISTA,” _Mecánica Computacional_ , vol. 31, no. 20, pp. 3255–3263, 2012.
  * [13] D. O. Pérez, D. R. Velis, and M. D. Sacchi, “High-resolution prestack seismic inversion using a hybrid FISTA least-squares strategy,” _Geophysics_ , vol. 78, no. 5, pp. R185–R195, 2013.
  * [14] C. Li, X. Liu, K. Yu, X. Wang, and F. Zhang, “Debiasing of seismic reflectivity inversion using basis pursuit de-noising algorithm,” _J. Applied Geophysics_ , vol. 177, p. 104028, 2020.
  * [15] F. Li, R. Xie, W.-Z. Song, and H. Chen, “Optimal seismic reflectivity inversion: Data-driven ${\ell_{p}}$-loss-${\ell_{q}}$-regularization sparse regression,” _IEEE Geoscience and Remote Sensing Lett._ , vol. 16, no. 5, pp. 806–810, 2019.
  * [16] C. Yuan and M. Su, “Seismic spectral sparse reflectivity inversion based on SBL-EM: experimental analysis and application,” _J. Geophysics and Engineering_ , vol. 16, no. 6, pp. 1124–1138, 2019.
  * [17] E. J. Candès, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $\ell_{1}$ minimization,” _J. Fourier Analysis and Applications_ , vol. 14, no. 5-6, pp. 877–905, 2008.
  * [18] T. Zhang, “Analysis of multi-stage convex relaxation for sparse regularization,” _J. Machine Learning Research_ , vol. 11, no. 3, 2010.
  * [19] I. Selesnick, “Sparse regularization via convex analysis,” _IEEE Trans. Signal Process._ , vol. 65, no. 17, pp. 4481–4494, 2017.
  * [20] S. J. Wright, R. D. Nowak, and M. A. Figueiredo, “Sparse reconstruction by separable approximation,” _IEEE Trans. Signal Process._ , vol. 57, no. 7, pp. 2479–2493, 2009.
  * [21] C.-H. Zhang, “Nearly unbiased variable selection under minimax concave penalty,” _The Annals of Statistics_ , vol. 38, no. 2, pp. 894–942, 2010\.
  * [22] J. Woodworth and R. Chartrand, “Compressed sensing recovery via nonconvex shrinkage penalties,” _Inverse Problems_ , vol. 32, no. 7, p. 075004, 2016\.
  * [23] J. Fan and R. Li, “Variable selection via nonconcave penalized likelihood and its oracle properties,” _J. American Statistical Association_ , vol. 96, no. 456, pp. 1348–1360, 2001.
  * [24] W. Zhong and J. Kwok, “Gradient descent with proximal average for nonconvex and composite regularization,” _Proc. AAAI Conf. Artificial Intelligence_ , vol. 28, no. 1, 2014.
  * [25] U. S. Kamilov, “A parallel proximal algorithm for anisotropic total variation minimization,” _IEEE Trans. Image Process._ , vol. 26, no. 2, pp. 539–548, 2016.
  * [26] H. H. Bauschke, R. Goebel, Y. Lucet, and X. Wang, “The proximal average: basic theory,” _SIAM J. Optimization_ , vol. 19, no. 2, pp. 766–785, 2008.
  * [27] Y. Yu, “Better approximation and faster algorithm using the proximal average,” _Adv. Neur. Info. Process. Syst._ , vol. 26, pp. 458–466, 2013\.
  * [28] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” _Proc. 27th Intl. Conf. on Machine Learning_ , pp. 399–406, 2010.
  * [29] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” _IEEE Signal Process. Mag._ , vol. 38, no. 2, pp. 18–44, 2021.
  * [30] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” _Communications on Pure and Applied Mathematics_ , vol. 57, no. 11, pp. 1413–1457, 2004.
  * [31] D. Mahapatra, S. Mukherjee, and C. S. Seelamantula, “Deep sparse coding using optimized linear expansion of thresholds,” 2017, arXiv:1705.07290 [Online].
  * [32] S. Mukherjee, D. Mahapatra, and C. S. Seelamantula, “DNNs for sparse coding and dictionary learning,” _Proc. NIPS Bayesian Deep Learning Workshop_ , 2017\.
  * [33] J. Zhang and B. Ghanem, “ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing,” _Proc. IEEE Conf. on Computer Vision and Pattern Recognition_ , pp. 1828–1837, 2018.
  * [34] M. Borgerding, P. Schniter, and S. Rangan, “AMP-inspired deep networks for sparse linear inverse problems,” _IEEE Trans. Signal Process._ , vol. 65, no. 16, pp. 4293–4308, 2017.
  * [35] H. Sreter and R. Giryes, “Learned convolutional sparse coding,” _Proc. IEEE Intl. Conf. on Acoustics, Speech and Signal Process._ , pp. 2191–2195, 2018\.
  * [36] R. Liu, S. Cheng, L. Ma, X. Fan, and Z. Luo, “Deep proximal unrolling: Algorithmic framework, convergence analysis and applications,” _IEEE Trans. Image Process._ , vol. 28, no. 10, pp. 5013–5026, 2019.
  * [37] P. K. Pokala, A. G. Mahurkar, and C. S. Seelamantula, “FirmNet: A sparsity amplified deep network for solving linear inverse problems,” _Proc. IEEE Intl. Conf. on Acoustics, Speech and Signal Process._ , pp. 2982–2986, 2019\.
  * [38] Y. Li, M. Tofighi, J. Geng, V. Monga, and Y. C. Eldar, “Efficient and interpretable deep blind image deblurring via algorithm unrolling,” _IEEE Trans. on Computational Imaging_ , vol. 6, pp. 666–681, 2020.
  * [39] N. Shlezinger, J. Whang, Y. C. Eldar, and A. G. Dimakis, “Model-based deep learning,” 2020, arXiv:2012.08405 [Online].
  * [40] P. K. Pokala, P. K. Uttam, and C. S. Seelamantula, “ConFirmNet: Convolutional FirmNet and application to image denoising and inpainting,” _Proc. IEEE Intl. Conf. on Acoustics, Speech and Signal Process._ , pp. 8663–8667, 2020\.
  * [41] L. Wang, C. Sun, M. Zhang, Y. Fu, and H. Huang, “DNU: deep non-local unrolling for computational spectral imaging,” _Proc. IEEE /CVF Conf. Comp. Vis. and Patt. Recog._ , pp. 1658–1668, 2020.
  * [42] B. Tolooshams, S. Dey, and D. Ba, “Deep residual autoencoders for expectation maximization-inspired dictionary learning,” _IEEE Trans. on Neural Networks and Learning Systems_ , 2020.
  * [43] D. Jawali, P. K. Pokala, and C. S. Seelamantula, “CORNET: Composite-regularized neural network for convolutional sparse coding,” _Proc. IEEE Intl. Conf. on Image Processing_ , pp. 818–822, 2020.
  * [44] B. Tolooshams, S. Mulleti, D. Ba, and Y. C. Eldar, “Unfolding neural networks for compressive multichannel blind deconvolution,” 2021, arXiv:2010.11391 [Online].
  * [45] S. Yuan and S. Wang, “Spectral sparse Bayesian learning reflectivity inversion,” _Geophysical Prospecting_ , vol. 61, no. 4, pp. 735–746, 2013\.
  * [46] K. J. Bergen, P. A. Johnson, M. V. de Hoop, and G. C. Beroza, “Machine learning for data-driven discovery in solid earth geoscience,” _Science_ , vol. 363, no. 6433, 2019.
  * [47] A. Adler, M. Araya-Polo, and T. Poggio, “Deep learning for seismic inverse problems: Toward the acceleration of geophysical analysis workflows,” _IEEE Signal Process. Mag._ , vol. 38, no. 2, pp. 89–119, 2021.
  * [48] W. Lewis and D. Vigh, “Deep learning prior models from seismic images for full-waveform inversion,” _Proc. SEG Technical Program Expanded Abstracts_ , pp. 1512–1517, 2017.
  * [49] A. Richardson, “Seismic full-waveform inversion using deep learning tools and techniques,” 2018, arXiv:1801.07232 [Online]. .
  * [50] O. Ovcharenko, V. Kazei, M. Kalita, D. Peter, and T. Alkhalifah, “Deep learning for low-frequency extrapolation from multioffset seismic data,” _Geophysics_ , vol. 84, no. 6, pp. R989–R1001, 2019.
  * [51] H. Sun and L. Demanet, “Extrapolated full-waveform inversion with deep learning,” _Geophysics_ , vol. 85, no. 3, pp. R275–R288, 2020.
  * [52] M. Araya-Polo, J. Jennings, A. Adler, and T. Dahlke, “Deep-learning tomography,” _The Leading Edge_ , vol. 37, no. 1, pp. 58–66, 2018.
  * [53] Y. Kim and N. Nakata, “Geophysical inversion versus machine learning in inverse problems,” _The Leading Edge_ , vol. 37, no. 12, pp. 894–901, 2018\.
  * [54] W. Wang, F. Yang, and J. Ma, “Velocity model building with a modified fully convolutional network,” _Proc. SEG Technical Program Expanded Abstracts_ , pp. 2086–2090, 2018.
  * [55] Y. Wu, Y. Lin, and Z. Zhou, “Inversionnet: Accurate and efficient seismic waveform inversion with convolutional neural networks,” _Proc. SEG Technical Program Expanded Abstracts 2018_ , pp. 2096–2100, 2018.
  * [56] F. Yang and J. Ma, “Deep-learning inversion: A next-generation seismic velocity model building method,” _Geophysics_ , vol. 84, no. 4, pp. R583–R599, 2019.
  * [57] A. Adler, M. Araya-Polo, and T. Poggio, “Deep recurrent architectures for seismic tomography,” _Proc. 81st EAGE Conference and Exhibition_ , vol. 2019, no. 1, pp. 1–5, 2019.
  * [58] V. Das, A. Pollack, U. Wollner, and T. Mukerji, “Convolutional neural network for seismic impedance inversion,” _Geophysics_ , vol. 84, no. 6, pp. R869–R880, 2019.
  * [59] M. J. Park and M. D. Sacchi, “Automatic velocity analysis using convolutional neural network and transfer learning,” _Geophysics_ , vol. 85, no. 1, pp. V33–V43, 2020.
  * [60] B. Wu, D. Meng, L. Wang, N. Liu, and Y. Wang, “Seismic impedance inversion using fully convolutional residual network and transfer learning,” _IEEE Geoscience and Remote Sensing Lett._ , vol. 17, no. 12, pp. 2140–2144, 2020.
  * [61] W. Wang and J. Ma, “Velocity model building in a crosswell acquisition geometry with image-trained artificial neural networks,” _Geophysics_ , vol. 85, no. 2, pp. U31–U46, 2020.
  * [62] R. Biswas, M. K. Sen, V. Das, and T. Mukerji, “Prestack and poststack inversion using a physics-guided convolutional neural network,” _Interpretation_ , vol. 7, no. 3, pp. SE161–SE174, 2019.
  * [63] M. Alfarraj and G. AlRegib, “Semi-supervised learning for acoustic impedance inversion,” _Proc. SEG Technical Program Expanded Abstracts_ , pp. 2298–2302, 2019.
  * [64] M. Araya-Polo, S. Farris, and M. Florez, “Deep learning-driven velocity model building workflow,” _The Leading Edge_ , vol. 38, no. 11, pp. 872a1–872a9, 2019.
  * [65] Y. Wang, Q. Ge, W. Lu, and X. Yan, “Seismic impedance inversion based on cycle-consistent generative adversarial network,” _Proc. SEG Technical Program Expanded Abstracts_ , pp. 2498–2502, 2019.
  * [66] L. Mosser, O. Dubrule, and M. J. Blunt, “Stochastic seismic waveform inversion using generative adversarial networks as a geological prior,” _Mathematical Geosciences_ , vol. 52, no. 1, pp. 53–79, 2020.
  * [67] Z. Zhang and Y. Lin, “Data-driven seismic waveform inversion: A study on the robustness and generalization,” _IEEE Trans. on Geoscience and Remote sensing_ , vol. 58, no. 10, pp. 6900–6913, 2020.
  * [68] D. P. Wipf and B. D. Rao, “Sparse Bayesian learning for basis selection,” _IEEE Trans. Signal Process._ , vol. 52, no. 8, pp. 2153–2164, 2004.
  * [69] S. Mache, P. K. Pokala, K. Rajendran, and C. S. Seelamantula, “DuRIN: A deep-unfolded sparse seismic reflectivity inversion network,” 2021, arXiv:2104.04704 [Online].
  * [70] M. A. T. Figueiredo, J. M. Bioucas-Dias, and R. D. Nowak, “Majorization–minimization algorithms for wavelet-based image restoration,” _IEEE Trans. Image Process._ , vol. 16, no. 12, pp. 2980–2991, 2007.
  * [71] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” _Adv. Neur. Info. Process. Syst._ , vol. 32, pp. 8024–8035, 2019\.
  * [72] D. Freedman, R. Pisani, and R. Purves, _Statistics (international student edition)_. WW Norton & Company, New York, 2007.
  * [73] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014, arXiv:1412.6980 [Online].
  * [74] W. Hamlyn, “Thin beds, tuning, and AVO,” _The Leading Edge_ , vol. 33, no. 12, pp. 1394–1396, 2014.
  * [75] G. S. Martin, R. Wiley, and K. J. Marfurt, “Marmousi2: An elastic upgrade for Marmousi,” _The Leading Edge_ , vol. 25, no. 2, pp. 156–166, 2006.
  * [76] H. Chung and D. C. Lawton, “Frequency characteristics of seismic reflections from thin beds,” _Canadian J. Exploration Geophysics_ , vol. 31, no. 1, pp. 32–37, 1995.
  * [77] R. Versteeg, “The Marmousi experience: velocity model determination on a synthetic complex data set,” _The Leading Edge_ , vol. 13, no. 9, pp. 927–936, 1994.
  * [78] O. S. R., dGB Earth Sciences, “Penobscot 3D - Survey,” 2017, data retrieved from https://terranubis.com/datainfo/Penobscot.
  * [79] E. Bianco, “Geophysical tutorial: Well-tie calculus,” _The Leading Edge_ , vol. 33, no. 6, pp. 674–677, 2014.

## Supplementary Material

In this document, we provide additional experimental validation, with results
on synthetic and simulated data for the proposed NuSPAN-$1$ and NuSPAN-$2$, in
comparison with the benchmark techniques, namely, basis-pursuit inversion
(BPI) [10, 9], fast iterative shrinkage-thresholding algorithm (FISTA) [11,
12] and expectation-maximization-based sparse Bayesian learning (SBL-EM) [68,
16].

###  -A Experimental Results and Discussion

In this section, we provide additional results to demonstrate the efficacy of
the proposed networks, NuSPAN-$1$ and NuSPAN-$2$, in comparison with the
benchmark techniques, BPI, FISTA, and SBL-EM, on synthetic $1$-D seismic
traces and $2$-D wedge models [74], and the simulated $2$-D Marmousi$2$ model
[75]. We quantify the performance based on objective metrics defined in
Section V-A of the main document, namely, Correlation Coefficient (CC),
Relative Reconstruction Error (RRE), Signal-to-Reconstruction Error Ratio
(SRER), and Probability of Error in Support (PES).

####  -A1 Testing Phase – Synthetic $1$-D Traces

Table VII shows the computational time during testing for NuSPAN-$1$ and
NuSPAN-$2$ in comparison with the benchmark techniques, which is computed over
$100$, $200$, $500$, and $1000$ test realizations of synthetic $1$-D traces.
The proposed networks, namely, NuSPAN-$1$ and NuSPAN-$2$, require lower
computational time in comparison to other techniques. FISTA, the next best to
our techniques, requires two orders of magnitude more computation time than
ours. Low computation times are significant in the context of reflection
seismic processing, where the amount of data to be handled is large. We note
that NuSPAN-$1$ and NuSPAN-$2$ are trained on a large number of synthetic
seismic traces before testing and have longer training times, but a much
shorter testing time.

Tables VIII, IX, and X compare the performance of the proposed NuSPAN-$1$ and
NuSPAN-$2$ with that of the benchmark techniques for sparsity factor $0.10$,
$0.15$, and $0.20$. The results for NuSPAN-$1$ and NuSPAN-$2$ for all the
tested sparsity factors are obtained by considering the model trained with
sparsity factor of $0.20$, whereas the parameters of the benchmark techniques
are tuned for each sparsity factor. This demonstrates the robustness of the
proposed NuSPAN to mismatch in sparsity factor, which is a significant
advantage in the context of seismic reflectivity inversion, as the accurate
estimation of the sparsity of the seismic reflections is challenging [15, 16].

TABLE VII: Computational time (in seconds) during testing computed over $100$,
$200$ $500$, and $1000$ realizations of synthetic $1$-D traces. This
comparison shows that the proposed techniques require significantly lower
compute time as compared with the benchmark methods during inference. The best
performance is highlighted in boldface. The second best performance scores are
underlined. BPI: Basis-Pursuit Inversion; FISTA: Fast Iterative Shrinkage-
Thresholding Algorithm; SBL-EM: Expectation-Maximization-based Sparse Bayesian
Learning; NuSPAN-$1$ and NuSPAN-$2$: Nonuniform Sparse Proximal Average
Network – Type $1$ and Type $2$.

Method | Time ($\mathrm{s}$) for # of test realizations
---|---
$\mathbf{100}$ | $\mathbf{200}$ | $\mathbf{500}$ | $\mathbf{1000}$
BPI | ${45.0384}$ | ${93.0857}$ | ${238.8020}$ | ${503.8519}$
FISTA | ${3.4676}$ | ${6.8286}$ | ${17.5604}$ | ${33.2539}$
SBL-EM | ${155.0183}$ | ${313.2058}$ | ${781.2321}$ | ${1076.1833}$
NuSPAN-$1$ | $\mathbf{0.0420}$ | $\mathbf{0.0627}$ | $\mathbf{0.1505}$ | $\mathbf{0.1778}$
NuSPAN-$2$ | $\underline{0.0493}$ | $\underline{0.0707}$ | $\underline{0.1536}$ | $\underline{0.1870}$

TABLE VIII: Metrics averaged over $1000$ test realizations of synthetic $1$-D
traces, for sparsity factor $\mathbf{0.10}$. NuSPAN-$2$ shows superior
amplitude recovery accuracy mentioned in terms of CC, RRE, and SRER.
NuSPAN-$1$ and NuSPAN-$2$ are trained with sparsity factor $0.20$ and tested
for sparsity factor $0.10$, whereas the benchmark techniques (BPI, FISTA, and
SBL-EM) are tuned and tested for the same sparsity factor of $0.10$.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | ${0.3764}$ | ${0.8590}$ | ${0.6749}$ | ${0.9424}$
FISTA | $\underline{0.3988}$ | $\underline{0.8383}$ | $\underline{0.7770}$ | $\mathbf{0.8880}$
SBL-EM | ${0.3827}$ | ${0.8520}$ | ${0.7036}$ | ${0.9424}$
NuSPAN-$1$ | ${0.3874}$ | ${0.8472}$ | ${0.7258}$ | $\underline{0.9130}$
NuSPAN-$2$ | $\mathbf{0.4042}$ | $\mathbf{0.8340}$ | $\mathbf{0.7950}$ | ${0.9131}$

TABLE IX: Metrics averaged over $1000$ test realizations of synthetic $1$-D
traces for sparsity factor $\mathbf{0.15}$. The proposed NuSPAN-$1$ and
NuSPAN-$2$ show superior amplitude and support recovery, although they are
trained with sparsity factor $0.20$ but tested on a different sparsity factor.
The best performance is highlighted in boldface. The second best performance
scores are underlined.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | ${0.3602}$ | ${0.8693}$ | ${0.6126}$ | ${0.9155}$
FISTA | ${0.3727}$ | ${0.8585}$ | ${0.6679}$ | ${0.8759}$
SBL-EM | ${0.3646}$ | ${0.8658}$ | ${0.6299}$ | ${0.9155}$
NuSPAN-$1$ | $\underline{0.3804}$ | $\underline{0.8528}$ | $\underline{0.6965}$ | $\mathbf{0.8725}$
NuSPAN-$2$ | $\mathbf{0.3904}$ | $\mathbf{0.8448}$ | $\mathbf{0.7377}$ | $\underline{0.8726}$

TABLE X: Metrics averaged over $1000$ test realizations of synthetic $1$-D
traces for sparsity factor $\mathbf{0.20}$. NuSPAN-$1$ and NuSPAN-$2$
outperform the benchmark techniques in both amplitude and support recovery.
Here, all the methods are tuned/trained with the same sparsity factor as they
are tested for.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | ${0.3550}$ | ${0.8735}$ | ${0.5910}$ | ${0.8895}$
FISTA | ${0.3619}$ | ${0.8687}$ | ${0.6144}$ | ${0.8781}$
SBL-EM | ${0.3240}$ | ${0.8999}$ | ${0.4602}$ | ${0.8899}$
NuSPAN-$1$ | $\underline{0.3750}$ | $\underline{0.8574}$ | $\underline{0.6724}$ | $\mathbf{0.8339}$
NuSPAN-$2$ | $\mathbf{0.3812}$ | $\mathbf{0.8523}$ | $\mathbf{0.6984}$ | $\underline{0.8340}$

####  -A2 Testing Phase – Synthetic $2$-D Wedge Models

In Section V-C2 of the main document, we have shown the results for one odd
(NP) and one even (NN) wedge model, where N and P denote the polarity,
negative and positive, respectively, of the two interfaces of the wedge
models. Here, we present results for two wedge model variants, one odd (PN –
Table XI and Figure 12) and one even (PP – Table XII and Figure 13).

TABLE XI: Results for a synthetic $2$-D odd (PN) wedge model. NuSPAN-$1$ shows
superior amplitude recovery in terms of CC, and outperforms the benchmark
techniques in support recovery quantified in terms of PES.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | ${0.8002}$ | ${0.2711}$ | $\mathbf{13.0565}$ | ${0.9933}$
FISTA | ${0.7984}$ | $\underline{0.2684}$ | $\underline{12.8473}$ | $\underline{0.7058}$
SBL-EM | $\underline{0.8096}$ | $\mathbf{0.2622}$ | ${10.3299}$ | ${0.9933}$
NuSPAN-$1$ | $\mathbf{0.8650}$ | ${0.2695}$ | ${8.1181}$ | $\mathbf{0.6019}$
NuSPAN-$2$ | ${0.7877}$ | ${0.3440}$ | ${5.3725}$ | ${0.9896}$

Figure 12:  Results for a synthetic $2$-D odd (PN) wedge model. True (a)
seismic traces and (c) reflectivity; (b) Noisy seismic traces with input SNR
$10$ dB. (d)-(f) Recovered reflectivity signatures for BPI, FISTA, and SBL-EM
show that these methods fail to resolve reflectors with wedge thickness
$<~{}5$ m, evident from the diverging interfaces highlighted by the rectangle
in black. (g)-(h) The recovered reflectivity profiles for NuSPAN-$1$ and
NuSPAN-$2$ show better resolution of reflector locations at wedge thickness
$<~{}5$ m. TABLE XII: Results for a synthetic $2$-D even (PP) wedge model.
NuSPAN-$1$ outperforms the benchmark techniques in support recovery in terms
of PES, whereas BPI, FISTA, and SBL-EM show better amplitude recovery accuracy
mentioned in terms of CC, RRE, and SRER.

Method | CC | RRE | SRER | PES
---|---|---|---|---
BPI | $\underline{0.8586}$ | $\underline{0.2608}$ | $\mathbf{12.9680}$ | ${0.9933}$
FISTA | $\mathbf{0.8791}$ | $\mathbf{0.2274}$ | $\underline{12.3515}$ | $\underline{0.6673}$
SBL-EM | ${0.8491}$ | ${0.3205}$ | ${10.559}$ | ${0.9933}$
NuSPAN-$1$ | ${0.8037}$ | ${0.3689}$ | ${8.6339}$ | $\mathbf{0.5680}$
NuSPAN-$2$ | ${0.7516}$ | ${0.4540}$ | ${5.45024}$ | ${0.9896}$

Figure 13:  Results for a synthetic $2$-D even (PP) wedge model. True (a)
seismic traces and (c) reflectivity; (b) Noisy seismic traces with input SNR
$10$ dB; (d)-(h) Recovered reflectivity. All methods resolve the even wedge
model well. NuSPAN-$1$ does not introduce spurious supports as the other
methods, seen as low-amplitude spikes in regions away from the two interfaces
in the results of the other methods.

####  -A3 Testing Phase – Simulated Marmousi$2$ model

Initial evaluations on the Marmousi$2$ model [75] showed very low PES for BPI
and SBL-EM (Table XIII), possibly due to the spurious support estimates
introduced by these methods complementing the low-amplitude spikes ($<~{}1\%$
of the absolute of the maximum amplitude) in the Marmousi$2$ model. When these
spikes were muted, the CC, RRE, and SRER for all the methods improved, but the
PES, now evaluated only on significant interfaces, was much higher for BPI and
SBL-EM, and lower for NuSPAN-$1$ (Table VI in the main document).

In Table XIII, we also provide the training and testing times of the benchmark
methods, namely, BPI, FISTA and SBL-EM, and the proposed NuSPAN-$1$ and
NuSPAN-$2$, on the complete Marmousi$2$ model. Comparing the three learning-
based methods, i.e., SBL-EM, NuSPAN-$1$, and NuSPAN-$2$, shows that the
proposed approaches have lower combined training and testing time than SBL-EM.
Although the combined time for FISTA is low, NuSPAN-$2$ achieves higher
accuracy in terms of both amplitude and support recovery after training on a
larger synthetic dataset.

TABLE XIII: Results for a portion of the Marmousi$2$ model corresponding to
Table VI and Figure 10 in Section V-D of the main document, without muting
low-amplitude spikes. NuSPAN-$2$ provides the best amplitude recovery. The
spurious supports introduced by BPI and SBL-EM coincide with some of the low-
amplitude spikes, leading to very low PES values. The reported times are for
the complete model (training times on one RTX 2080 Ti GPU). Training times are
applicable only to NuSPAN-$1$ and NuSPAN-$2$. The best performance is
highlighted in boldface. The second best performance scores are underlined.

Method | CC | RRE | SRER | PES | Time ($\mathrm{h}$)
---|---|---|---|---|---
Training | Testing | Training + Testing
BPI | ${0.9472}$ | ${0.0878}$ | ${14.7626}$ | $\mathbf{0.0181}$ | - | $16.2671$ | ${16.2671}$
FISTA | ${0.9406}$ | ${0.1021}$ | ${13.8227}$ | ${0.912}$ | - | $5.5636$ | $\mathbf{5.5636}$
SBL-EM | $\underline{0.9548}$ | $\mathbf{0.0686}$ | $\mathbf{18.1117}$ | $\mathbf{0.0181}$ | - | $101.2659$ | $101.2659$
NuSPAN-$1$ | ${0.9375}$ | ${0.1035}$ | ${14.2261}$ | ${0.9650}$ | $\underline{9.8069}$ | $\mathbf{0.0737}$ | $9.8806$
NuSPAN-$2$ | $\mathbf{0.9590}$ | $\underline{0.0718}$ | $\underline{15.0579}$ | $\underline{0.2021}$ | $\mathbf{8.6007}$ | $\underline{0.1024}$ | $\underline{8.7031}$