# A hybrid approach to seismic deblending: when physics meets self-supervision

Nick Luiken, Matteo Ravasi, Claire E. Birnie
King Abdullah University of Science and Technology,
$\left\\{\texttt{nicolaas.luiken, matteo.ravasi,
<EMAIL_ADDRESS>
###### Abstract

To limit the time, cost, and environmental impact associated with the
acquisition of seismic data, in recent decades considerable effort has been
put into so-called simultaneous shooting acquisitions, where seismic sources
are fired at short time intervals between each other. As a consequence, waves
originating from consecutive shots are entangled within the seismic
recordings, yielding so-called blended data. For processing and imaging
purposes, the data generated by each individual shot must be retrieved. This
process, called deblending, is achieved by solving an inverse problem which is
heavily underdetermined. Conventional approaches rely on transformations that
render the blending noise into burst-like noise, whilst preserving the signal
of interest. Compressed sensing type regularization is then applied, where
sparsity in some domain is assumed for the signal of interest. The domain of
choice depends on the geometry of the acquisition and the properties of
seismic data within the chosen domain. In this work, we introduce a new
concept that consists of embedding a self-supervised denoising network into
the Plug-and-Play (PnP) framework. A novel network is introduced whose design
extends the blind-spot network architecture of [28] for partially coherent
noise (i.e., correlated in time). The network is then trained directly on the
noisy input data at each step of the PnP algorithm. By leveraging both the
underlying physics of the problem and the great denoising capabilities of our
blind-spot network, the proposed algorithm is shown to outperform an industry-
standard method whilst being comparable in terms of computational cost.
Moreover, being independent on the acquisition geometry, our method can be
easily applied to both marine and land data without any significant
modification.

##  1 Introduction

Reflection seismology [40] is a geophysical technique that uses reflected
seismic waves to characterize the Earth’s subsurface. This method comprises of
a controlled source of seismic energy and an array of receivers that record
the pressure (or displacement) induced by the reflected waves. The physical
instrumentation varies depending on the acquisition environment: for example,
an airgun is used as seismic source in marine settings whilst a vibroseis is
the preferred choice in land settings. Similarly, hydrophones (and possibly
accelerometers) are towed behind a boat inside a so-called streamer cable or
geophones (combined with hydrophones) can be laid directly on the seafloor for
marine acquisition; whilst geophones are used for land acquisition. After the
introduction of 3D seismic [11], today’s conventional seismic acquisition
campaigns may last several weeks up to a few months [7, 26, 10]. In an attempt
to improve acquisition efficiency, and therefore limit the time, cost, and
associated environmental impact, [6, 9, 3, 36] introduced a new paradigm in
seismic acquisition usually referred to as simultaneous shooting. Simply put,
consecutive sources are fired at short time intervals, thereby minimizing the
overall acquisition time. This comes at the cost of recording entangled
seismic data, also called blended data, where the waves originating from one
source tend to overlap with those originating from previous and subsequent
sources. To render such data suitable for subsequent steps of seismic
processing and imaging, the interference between consecutive shots must be
removed such that the contribution of each individual source (also referred to
as a shot gather) is retrieved. This process is called deblending.
Mathematically speaking, blending is an operator that takes the otherwise
conventionally acquired data and combines them together using a pre-defined
sequence of time delays. In theory, deblending can be achieved by solving the
associated inverse problem; however, as this problem is heavily
underdetermined, choosing an appropriate regularization is fundamental to
achieve a successful inversion. Historically, the design of suitable
regularizers is motivated by the effect of the adjoint of the blending
operator on the blended data. In fact, the resulting data can be seen as a
superposition of coherent signal (i.e, reflections from the shot whose firing
time has been properly accounted for) and trace-wise, burst-like noise (i.e,
reflections from all other interfering shots whose firing times have not been
properly accounted for).
The recent success of deep learning in various scientific disciplines has
attracted the interest of the geophysical community, resulting in many
opportunities and new challenges [50]. A common hurdle when applying deep
learning methods to geophysical inverse problems is the lack of trustworthy
training data. More specifically, whilst training data should consist of
clean, representative ground truth examples that resemble the solution to the
inverse problem at hand, such a solution generally does not exist a priori.
Two approaches commonly adopted to circumvent this problem are to either
generate synthetic data or to use state-of-the-art algorithms to produce
input-output pairs to train a network on; in both cases, transfer learning
[41, 34] or domain adaptation [2, 12] techniques are then required to
generalize the network capabilities to unseen field data. A major drawback of
this approach is that synthetic data may not resemble field data accurately
enough to be considered a representative dataset: this is well-known in the
geophysical community and has been a major criticism for decades when new
methods are tested only on synthetic data. This also represents a serious
roadblock to the application of deep learning methods in geophysics.
Additionally, in most geophysical applications the underlying physics, or at
least part of the physics, is well understood. Pure, end-to-end machine
learning methods tend to ignore these well-studied physical principles,
thereby discarding important a priori knowledge of the problem they are tasked
to solve.

#### Our contribution

We introduce a novel algorithm for seismic deblending, which blends the
physics of the underlying physical process with a state-of-the-art self-
supervised denoiser into a single, well-crafted inverse process. This is
specifically achieved within the framework of Plug-and-Play (PnP) priors. Our
network architecture is inspired by the blind-spot network of [28] and
modified to handle trace-wise coherent noise. The network is trained on-the-
fly at each PnP iteration in a self-supervised manner, completely bypassing
the need for ground truth data. Our numerical experiments illustrate that the
proposed algorithm can outperform a state-of-the-art conventional method.
Finally, we show that our algorithm is independent on the underlying structure
of the seismic data and can be used easily for different acquisition set-ups -
a clear advantage over conventional methods.

##  2 Background

#### The seismic data layout

Seismic data are commonly acquired by firing a source at a given time and
recording the reflections arising from the interaction between the emitted
seismic wave and changes in subsurface properties. Conceptually, seismic data
can be arranged as a three dimensional tensor (or a cube), having the
dimensions of the number of sources $n_{s}$, number of receivers $n_{r}$, and
number of time samples $n_{t}$: $d_{c}(x_{s},x_{r},t)$. Slicing this cube in
different directions give raise to so-called seismic gathers: more
specifically, when slicing across the source axis, we obtain the data recorded
by all receivers for a single shot, usually called a Common Shot Gather (CSG);
conversely, by slicing across the receiver axis we obtain the data generated
by all shots for a single receiver. When the receivers move alongside the
source (i.e., marine case) the resulting gather is called the Common Channel
Gather (CCG). For static receivers (i.e., ocean-bottom or land acquisition),
the seismic gather is know as the Common Receiver Gather (CRG). Both scenarios
will later be considered.

#### Blended acquisition

In practice, to be able to collect data where no overlap exists between
consecutive shots, each shot has to be fired with an appropriate time delay,
such that all reflections from one shot have been recorded by the receivers
before the next shot is fired. This dictates the overall acquisition time and
greatly limits any possible acquisition speed-up. Alternatively, in blended
acquisition, shots are fired at shorter intervals. This means that each
individual CSG contains recordings from both the nominal as well as the
previous and subsequent shots. In this work, we consider the so-called
continuous blending setting. This approach is state-of-the-art in marine
seismic acquisition due to the fact it is easy to implement in the field. It
is in fact simply achieved by firing the airgun towed by the acquisition
vessel at short time intervals, and continuously recording the waves returning
to the receiver array as depicted in figure 1. The recorded data $d_{b}$ can
be simply described as the superposition of all of the unblended, or clean,
data shifted in time by the given time delay $t_{i}=i\cdot T+\Delta t_{i}$.
Here, $T$ is the nominal firing interval and $\Delta t_{i}$ is a random dither
applied to the nominal firing time of each shot. The blended data can thus be
described as a function of the clean data

$d_{b}=Bd_{c}:=\left[B_{1},\ldots,B_{n_{s}}\right]d_{c}=B_{1}d_{c,1}+\ldots+B_{n_{s}}d_{c,n_{s}}$
(1)

where the blending operator is a horizontal stack of time-shift operators
$B_{i}$, and the clean data is a vector where all vectorized shot gathers,
$d_{c,i}=vec(d_{c}(x_{s,i},x_{r},t))$, are stacked together. Note that each
$B_{i}$ time-shift operator has the property that $B_{i}^{H}B_{i}=I$, and a
composition of time-shift operators is again a time-shift operator.

Figure 1: Schematic illustration of a seismic simultaneous shooting
acquisition. a) Cartoon of a seismic acquisition campaign in continuous
blending mode. A single vessel towing a source (red star) and an array of
receivers (blue triangles) moves from right to left and fires energy into the
ground at dithered periodic time samples. For each shot, reflections
originated from shallow subsurface layers are immediately recorded by the
receivers, whilst those produced by deeper reflectors are recorded later in
time alongside the shallow reflections from the next firing shot. This
phenomenon leads to the blending of independent shot gathers. b) A short time
window of the continuously blended seismic data. Dashed vertical color lines
represent the nominal firing times (i.e., $i\cdot T$), whilst the solid color
lines represent the actual firing times with dithering. Color rectangles refer
to every individual shot gather that we wish to separate from the other
overlapping gathers. c) Pseudo-deblended data for a single receiver (white
dashed line in panel b).

#### Pseudo-deblending

To better understand how to design effective regularization strategies for the
deblending problem, we first have to consider the action of $B^{H}$ on the
blended data. For the $i^{th}$ shot gather, the result of $B_{i}^{H}Bd_{c}$
can be written as

$\displaystyle
B_{i}^{H}\\!\left(B_{1}d_{c,1}+\ldots+B_{n_{s}}d_{c,n_{s}}\right)$
$\displaystyle=d_{c,i}+\left(B_{i}^{H}\\!B_{1}d_{c,1}+\ldots+B_{i}^{H}\\!B_{n_{s}}d_{c,n_{s}}\right)$
(2)

Therefore the action of $B_{i}^{H}$ on the blended data produces the original
$i^{th}$ shot gather alongside randomly shifted versions of all the other shot
gathers. Whilst these randomly shifted shot gathers are coherent and look like
standard seismic signal in the CSG domain, they appear as trace-wise coherent
noise in the CRG (or CCG) domain as shown in the Experiments section. Because
the application of the adjoint of the blending operator retrieves the true
signal, albeit with some additional noise, its action is usually called
pseudo-deblending. As a consequence of this, it is now clear that to retrieve
the various $d_{c,i}$, an effective regularization must filter the trace-wise
noise in CRGs (or CCGs) whilst preserving the coherent signal. Conventional
approaches identify a domain in which the signal can be easily discriminated
from the noise, and more specifically the signal in such domain is sparse
whilst the noise is not. Examples of such a kind are the hyperbolic Radon
transform for CRGs [25], the patched Fourier transform for CCGs [1], or the
Curvelet transform [31]. It should be noted here that the performance and
computational cost of such conventional approaches varies greatly. As
explained later, a self-supervised denoiser can also be empowered to filter
this noise, and, in fact, its denoising capabilities will be shown to be
superior to those of commonly used, hand-crafted priors.

#### Deblending by denoising

Deblending by denoising aims to deblend the data by directly denoising the
pseudo-deblended CRGs. Deblending by denoising can be justified by noting that
$BB^{H}$ is a diagonal matrix. Therefore, making the substitution
$d_{c}=B^{H}\\!z$ leads to

$Bd_{c}=BB^{H}\\!z=Dz,$

where $D$ is a diagonal matrix. The structure of $D$ depends on the blending
strategy. For continuous blending, the diagonal elements $d_{ii}$ count the
number of overlapping shots in each time-space sample of the blended data,
which usually ranges from two to three. We can therefore approximate $D=I$,
based on the argument that $D$ is a simple amplitude correction. The solution
is then given by

$z=d_{b}\Longrightarrow d_{c}=B^{H}d_{b}.$

By doing so, the clean data, $d_{c}$, is approximated by the pseudo-deblended
data, $d_{b}$, which however still contains the blending noise. Therefore, we
can solve the following denoising problem to fully deblend the data:

$\min_{d_{c}}\|d_{c}-B^{H}d_{b}\|+\mathcal{R}(d_{c}),$ (3)

where $\mathcal{R}(\cdot)$ is any chosen regularization. Because the
distribution of the blending noise is far from being Gaussian, a typical
choice for the data fidelity norm is $\|\cdot\|_{1}$. However, the
$\ell_{1}$-norm in the data fidelity renders the solution non-trivial. Note
that while in principle we could choose to denoise $z$, in this domain we do
not have an appropriate regularization term, making it impossible to denoise.

#### Deblending by inversion

Alternatively, one could retrieve the clean data by solving the (heavily)
underdetermined inverse problem,

$\min_{d_{c}}\frac{1}{2}\|Bd_{c}-d_{b}\|_{2}^{2}+\mathcal{R}(d_{c}).$ (4)

As the incoherent noise is fully explained by the blending operator, the
$\ell_{2}$-norm for the data fidelity is an appropriate choice. From the above
equation, it is not obvious how the noise appears in the data since we are not
dealing directly with the pseudo-deblended data. The noise is however
introduced into the solution through the gradient of the data term,
$B^{H}(\\!Bd_{c}-\\!d_{b})$. In sparse inversion, the objective is minimized
iteratively by thresholding the gradient at every iteration, which can be
interpreted as a noise removal step. The literature has shown that deblending
by inversion is superior in comparison to deblending by denoising in terms of
the overall quality of reconstruction.

#### Self-supervised denoising: Incoherent noise

Self-supervised denoising is a recently introduced concept in the deep
learning literature. This family of methods are designed in such a way that
noisy images can be used as both the input and label to train a neural network
to act as a denoiser, thereby bypassing the need for clean data as labels.
Noise2Noise represents the first such method not relying on ground truth
labels [29]. The network is forced to infer the signal from pairs of noisy
data. For applications where such pairs are unavailable, an alternative was
proposed in the concurrent works of [27] and [5], who introduced Noise2Void
and Noise2Self, respectively. In both cases, the same image is used as input
and label: under the assumption that the noise is incoherent whilst the signal
is coherent, the network can naturally learn to infer only the signal from its
neighbouring pixels. This mechanism bears similarity to classical statistical
denoisers such as non-local means [15] or BM3D [18], however it relies on a
more flexible nonlinear mapping using learned convolutional filters. More
specifically, to denoise a particular pixel, [27] replace the pixel of the
input image with a randomly selected neighbouring pixel. As this pre-
processing step introduces randomicity in the central pixel of the receptive
field of the network, the network should not learn anything from it and
naturally learns to infer the signal from its neighbours (since the noise is
assumed to be incoherent). Rather than directly replacing the pixel of
interest, [5] pre-process the input image with a blind-spot (or donut)
convolutional filter, so that the network once again cannot rely on the
central pixel to predict itself. A key limitation of both approaches lies in
the fact that the self-supervised loss can be evaluated only at the pixels
that have been corrupted, making the training of these denoisers relatively
slow. An alternative approach to blind-spot networks was introduced by [28].
Instead of corrupting the middle pixel, their network is explicitly designed
to have a receptive field with a hole in the middle. This is achieved by
combining padding and cropping with a standard convolution layer (i.e., to
create a causal filter) and by rotating the input image four times prior to
feeding it through the network. After the rotated inputs have been fed through
the network, they are rotated back, concatenated, and combined by a series of
$1\times 1$ convolutions prior to evaluating the loss at every pixel of the
output image. A schematic description of this network is depicted in figure
2a.

#### Self-supervised denoising: Coherent noise

Both Noise2Void and Noise2Self operate under the assumption that the noise is
independent and identically distributed. [14] shows that the denoising quality
of Noise2Void is degraded when the noise is structured. This shortcoming of
Noise2Void is solved by masking pixels along the direction of the noise: the
authors dub their method Structured Noise2Void. For the seismic blending
problem, the noise that we are interested to suppress is also structured: more
specifically, the blending noise shows correlation along the time axis. We,
therefore, extend here the efficient implementation of [28] to suppress
structured noise in seismic data, by simply using the original and flipped
(over the time axis) version of the image as input. This produces a network
whose receptive field is masked over an entire time trace, see figure 2b. In
the following, we will call this network StructBS.

Figure 2: (a) The blind-spot network of [28], whose receptive field excludes
the center pixel. (b) Our newly proposed blind-spot network, whose receptive
field excludes an entire direction instead of just the middle pixel. Impulse
responses are created by feeding the respective networks with unitary weights
and zero biases with a image containing a unitary spike in the middle.

##  3 Related work

#### Simultaneous shooting

Simultaneous shooting was first pioneered by [6] and has gained popularity in
recent years [9]. Although the first attempts at deblending were mostly by
means of denoising [33], recent research has reveled the superiority of
deblending by inversion [1]. Since then, research has been devoted to finding
appropriate regularization terms. Some approaches involve median-filtering
[23, 22, 24], rank-reduction methods [17, 54], sparse regularization [30, 31,
55, 56, 58], and deep learning [42, 57, 48]. All the deep learning approaches
to date use CNNs and require pre-training. [4] uses the RED framework
introduced in [37], which is similar to the PnP framework. The difference is
that RED explicitly incorporates the denoiser into the objective function. The
authors propose the use of two conventional regularization techniques as a
denoiser, the patched Fourier transform [1] and the singular-spectral-analysis
filter [17], instead of applying them as a sparse penalty.

#### Self-supervised seismic denoising

Seismic data are a prime example of a noisy data type where no clean, ground
truth labels are available. As such, the application of self-supervised
denoisers has recently been proposed for the suppression of different types of
noise present in seismic data. Following the Noise2Void methodology, [13] use
blind-spot networks for the suppression of random noise in post-stack seismic
data. Expanding on this, [32] adapted the methodology of StucturedNoise2Void
[14] for the suppression of trace-wise noise in seismic shot gathers,
originating from poorly coupled receivers and/or dead sensors. The method that
is most closely related to ours is the one presented in [48] \- both with
respect to application and methodology. The authors propose to use a self-
supervised denoising network to deblend the data by denoising. To produce
satisfactory results they require a number of additional pre- and post-
processing steps. In our work, we incorporate a deep learning based denoiser
in deblending by inversion, thereby leveraging both the underlying physics and
the power of neural networks. Moreover, no pre- and post-processing is
required.

#### The Plug-and-Play framework

The Plug-and-Play framework was pioneered by [47]. The authors considered a
number of popular denoisers, including BM3D [19], K-SVD [20], PLOW [16] and
q-GGMRF [44]. In subsequent works, the denoisers have been replaced by pre-
trained neural networks, most notably CNN and DnCNN. An extensive list of
references is provided in [51], and include [37, 53, 35, 45, 21, 46, 30, 43,
52]. This research focuses mostly on properly training the neural networks
such that they generalize well to different noise levels. The novelty of our
method is that there is no pre-training with clean data required.

##  4 Method

Equipped with a self-supervised denoiser, a straightforward approach to
deblending is to directly denoise the pseudo-deblended data. However,
deblending by denoising is known to be sub-optimal in comparison to deblending
by inversion. On the other hand, because a denoiser cannot be naturally added
as a constraint to the objective function in equation 4, it is not immediately
clear how to incorporate the denoiser into the inversion process.

#### PnP derivation

The PnP framework is derived from the ADMM algorithm, which we will shortly
derive here. The ADMM algorithm is generally used to solve inverse problems of
the form

$\min_{x}\mathcal{D}(\mathcal{M}(x),d)+\mathcal{R}(x),$

where $\mathcal{M}$ is the forward model, $d$ is the measured data,
$\mathcal{D}$ is a data fidelity term that is generally smooth, and
$\mathcal{R}$ is a convex, possibly non-smooth regularization term. To account
for the non-smoothness of $\mathcal{R}$, an auxiliary variable $y=x$ is
introduced, yielding the optimization problem

$\min_{x,y}\mathcal{D}(\mathcal{M}(x),d)+\mathcal{R}(y)\textnormal{ subject to
}x=y.$

The ADMM solves this problem by forming the so-called augmented Lagrangian,

$\max_{u}\min_{x,y}\mathcal{D}(\mathcal{M}(x),d)+\mathcal{R}(y)+\frac{\rho}{2}\|x-y\|^{2}_{2}+u^{T}(x-y),$

where $u$ is the Lagrange multiplier and $\rho$ is a scalar. This problem is
solved by alternatively minimizing over $x$ and $y$, and maximizing over $u$.
This yields the following scheme:

$\displaystyle x_{k+1}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{x}\left\\{\mathcal{D}(\mathcal{M}(x),\>d)+\frac{\rho}{2}\|x-y_{k}+u_{k}\|_{2}^{2}\right\\}$
$\displaystyle y_{k+1}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{y}\left\\{\mathcal{R}(y)+\frac{\rho}{2}\|x_{k+1}-y+u_{k}\|_{2}^{2}\right\\}$
$\displaystyle u_{k+1}$ $\displaystyle=$ $\displaystyle
u_{k}+x_{k+1}-y_{k+1}.$

The key observation of [47] is that the $y$-update can be interpreted as a
denoising inverse problem. As such, the authors propose to drop the user-
defined regularization $\mathcal{R}(\cdot)$ and instead plug in a denoiser of
choice in the $y$-update of the ADMM iterations. Equipped with a powerful
self-supervised denoiser that can filter the blending noise, we embed this
into the PnP framework. Our proposed algorithm reads as follows:

$\displaystyle x_{k+1}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{x}\left\\{\frac{1}{2}\|Bx-
d_{b}\|_{2}^{2}+\frac{\rho}{2}\|x-y_{k}+u_{k}\|_{2}^{2}\right\\}$
$\displaystyle y_{k+1}$ $\displaystyle=$
$\displaystyle\textnormal{StructBS}_{\theta}(x_{k+1}+u_{k})$ $\displaystyle
u_{k+1}$ $\displaystyle=$ $\displaystyle u_{k}+x_{k+1}-y_{k+1}.$

where $x$ is used here for simplicity in place of $d_{c}$, and the $x$-update
is solved by means of LSQR.

##  5 Experiments

In the following, our algorithm is tested on the openly available Mobil AVO
viking graben line 12 marine dataset
111https://wiki.seg.org/wiki/Mobil_AVO_viking_graben_line_12. As the data has
been originally acquired in a conventional fashion, we create the blending
operator and blend the data ourselves. In addition to containing all the
challenging features of a field dataset, this also provides us with a ground
truth $d_{c}$ onto which to assess the quality of our reconstruction. In this
example, the original dataset is composed of $n_{s}=64$ sources, $n_{r}=120$
receivers, and $n_{t}=1024$ samples (i.e., the total recording time per shot
equals 4 seconds). For the continuous blending operator, we choose a fixed
firing interval of $T=2$ seconds, with added random delays selected uniformly
in the interval $\Delta t_{i}\sim[-1,1]$ seconds. This overlap is quite
challenging as generally half of the signal overlaps with either that of the
previous or that of the next shot. Moreover, some pseudo-deblended shot
gathers exhibit contributions from three consecutive shots. Finally, the
relative mean-square error,
$RMSE=\|d_{c}-d_{c,true}\|_{2}/\|d_{c,true}\|_{2}$, is chosen as a metric of
comparison in all of our numerical examples. All experiments are performed on
a Intel(R) Xeon(R) CPU @ 2.10GHz equipped with a single NVIDIA GEForce RTX
3090 GPU.

#### Deblending by inversion with FISTA

To begin with, our newly proposed methodology is compared with the state-of-
the-art deblending algorithm of [1] that solves the deblending problem as a
sparsity promoting inversion

$z_{\star}=\operatorname*{arg\,min}_{z}\|BFz-
d_{b}\|_{2}^{2}+\lambda\|z\|_{1},\quad d_{b}=Fz_{\star},$ (5)

where $F$ is a linear operator that performs a patched two-dimensional Fourier
transform, and $z_{\star}$ is the solution in the Fourier domain that is
ultimately transformed back to the original time-space domain of the seismic
data. In our experiment, the size and number of patches as well as the
regularization parameter $\lambda$ are selected by hand to yield optimal
results. Moreover, the FISTA [8] solver is used with an adaptive decreasing
sequence, $\lambda_{k}$ (as this has been shown in the literature to
outperform a fixed $\lambda$ for this specific problem). We choose the
sequence $\lambda_{k}=\left(\frac{6}{5}e^{-0.05k}+6\right)\lambda_{0}$, which
was once again fine-tuned to give the best performance. The final error is
roughly $9.8\%$; hereon in this represents the benchmark against which we will
assess the effectiveness of our self-supervised PnP algorithm.

#### Deblending by inversion with PnP

Next, our PnP algorithm is applied to the same dataset. We choose 30 outer
iterations, 3 inner iterations, $\rho=1$ and 30 denoiser epochs. The choice of
these hyperparameters will be justified in the ablation study. We also use the
U-Net architecture in [28] and the $L_{1}$ norm for the self-supervised
training loss because it is more appropriate for burst-like noise. Since the
denoiser is trained on all the CCGs, the size of our training data is 120 and
we use a batch size of 8. This leads to a solution that has an overall error
of roughly $6.7\%$, which is approximately $3\%$ lower than the conventional
method. Figure 3a displays the relative error as function of the FISTA
iterations. One of the patches into which the pseudo-deblended data has been
decomposed alongside its Fourier spectrum are shown in figure 3d and 3b,
respectively. The presence of burst-like noise translates into a non-sparse
Fourier spectrum. On the other hand, the noise is removed in the deblending
results in 3e as a result of a much sparser Fourier spectrum obtained during
the inversion process as shown figure 3c.

Figure 3: a) Error for the patched Fourier inversion. b) Extracted patch in
the Fourier domain for the pseudo-deblended data. c) Extracted patch in the
Fourier domain for the deblended data. d) Extracted patch of the pseudo-
deblended data. e) Extracted patch of the deblended data.

As a visual comparison, figure 4 displays the results for a given CCG (top)
and CSG (bottom) for both the conventional and proposed approaches. It is
noteworthy that our algorithm shows a clear improvement in terms of denoising
capabilities, as visible in the displayed CCG. Especially after $t=2s$, where
the signal is weak and blending noise dominates, the conventional approach
tends to be more prone to signal leakage compared to our PnP algorithm.
Finally, the computational cost of the conventional algorithm can be
quantified in terms of the number of forward and adjoint operations for both
the blending ($B$) and patched Fourier ($F$) operators: in our example, this
amounts to 200 forward and adjoint passes. On the other hand, our method
requires 90 forward and adjoint passes for the blending operator and a total
of 900 training epochs for the network. Considering that all computations
(apart from the network related ones) are performed on the CPU, the two
algorithms are comparable in terms of overall computational time (2h and
34mins for the conventional algorithm and 1h and 51mins for the PnP
algorithm).

Figure 4: Deblending results for one CCG (top) and CSG (bottom). Although both
algorithms can successfully remove most of the blending noise, our algorithm
is less prone to signal leakage and provides better amplitude fidelity - a key
factor in seismic data processing.

#### PnP solution progression

To show how the solution progresses as a function of outer iterations, we
select a receiver gather and show $x_{k}$ in figure 5.

Figure 5: Progression of $x_{k}$ for one receiver gather.

We clearly see the noise level drop, which explains the need to keep training
the network as shown in the paper. Lastly, we show a denoised receiver gather
after every epoch of training, for outer iterations 1, 10, 20 and 30 in
figures 12 13, 14 and 15 in the appendix.

#### PnP on ocean-bottom cable scenario

In this section we show that our algorithm generalizes well to other
acquisition geometries. More specifically, we mimic here an ocean-bottom cable
acquisition scenario where both CSGs and CRGs contain hyperbolic events.
Despite the denoising process is here applied on data that presents a
different structure when it comes to the coherent signal, our PnP algorithm is
still very successful and achieves an overall error of $6.7\%$. We choose 3
inner iterations, $\rho=1$, and 20 denoiser epochs and let the algorithm run
for a total of 40 outer iterations. The reconstructions are shown in figure 6.

Figure 6: Deblending results for one CRG (top) and CSG (bottom) in ocean-
bottom configuration.

###  5.1 Ablation study

This section provides an extensive analysis of some of the key components of
the proposed PnP methodology and their impact on the overall solution of the
deblending inverse problem.

#### PnP iterations

To begin with, we assess the importance of the PnP iterations compared to
simply training the self-supervised denoiser on pseudo-deblended data and
applying it directly to the entire dataset. Although not shown here, the
result of this one-shot denoising produces a solution with an overall error of
roughly $19\%$. This is much worse than both the conventional and PnP method
and therefore considered not suitable.

#### The $x$-update

As previously explained, the $x$-update requires the solution of the linear
system

$\min_{x}\frac{1}{2}\left\|\begin{bmatrix}B\\\
\sqrt{\rho}I\end{bmatrix}x-\begin{bmatrix}d\\\
\sqrt{\rho}(y_{k}-u_{k})\end{bmatrix}\right\|_{2}^{2},$ (6)

which can be efficiently accomplished via LSQR. The convergence rate of LSQR
depends on the spectrum of the blending operator $B$, specifically on its
condition number [39, section 6.11.3]. As the singular values of $B$ are the
number of overlapping shots at each time step in the blended data, in our
experiment the condition number is 3. Nevertheless, the number of inner
iterations represents a hyperparameter that should be assessed. In order to do
so, $\rho$ and the number of epochs used to train the denoiser are fixed
whilst the number of inner iterations for our PnP algorithm is varied between
1, 3, and 5. In all cases, the relative error as a function of outer
iterations is computed and shown in figure 7a.

Figure 7: a) The error for fixed $\rho$ and denoiser epochs as a function of
the inner iterations. b) Error for fixed number of inner iterations and
denoiser epochs and variable $\rho$.

As expected, increasing the number of inner iterations does not improve the
overall solution. Somewhat surprisingly, the number of inner iterations
exhibits a regularizing behaviour, in the sense that there seems to be an
optimal number (here 3) above and below which the solution is poorer. When
using fewer iterations, the overall error decreases slower in the first few
outer iterations and quickly plateaus at around $20\%$. On the other hand,
when using more iterations, the initial convergence is as fast as that of the
optimal value, however the overall solution is of poorer quality.
Another important hyperparameter is the augmented Lagrangian scalar $\rho$.
For our algorithm, $\rho$ takes the role of a regularization parameter that
controls the discrepancy between $x_{k+1}$ and $y_{k}+u_{k}$. We test three
values, $\rho=0.1,1\textnormal{ and }10$, and choose the number of denoiser
epochs that gives the lowest error. Again, we compare the relative error as a
function of outer iterations (figure 7b). In the very first $x$-update,
$y_{k}$ and $u_{k}$ are 0. Therefore, the linear system in equation 6 amounts
to solving a Tikhonov regularized problem with regularization parameter
$\rho$. For $\rho=10$, the solution shrinks to values close to zero because
$\sigma_{\text{max}}(B)\ll\rho$. Therefore, to obtain a meaningful result, we
set $\rho=0$ in the first outer iteration and then switch to $\rho=10$.
Nevertheless, looking at the error curve, it is clear that $\rho=10$ is too
large: this is not unexpected because $\sigma_{\text{max}}(B)\ll\rho$ and the
regularization terms dominates the inversion at each outer step. Clearly,
$\rho=1$ is the best choice, which is interesting because it is also the value
that lies in the same order of magnitude of the singular values of $B$. This
implies that there is a perfect balance between data misfit and
regularization, which seems to be beneficial to the PnP algorithm. As stated
before, the parameter $\rho$ also controls the discrepancy between $x$ and
$y$. Theoretically, if the PnP algorithm were to converge, $u\to u_{\star}$ as
$k\to\infty$, then $x_{k}=y_{k}$ as $k\to\infty$. Therefore, the difference
between $x_{k}$ and $y_{k}$ is a good measure of the convergence of the
algorithm. Figure 8 shows the progression of $x_{k}$ and $y_{k}$ over the
number of outer iterations.

Figure 8: RMSE of $x_{k}$ and $y_{k}$ for different $\rho$.

In the scenario where $\rho$ is small, $x_{k}$ and $y_{k}$ are not close,
indicating that the algorithm does not converge. Interestingly, $x_{k}$ has a
lower error than its denoised counterpart $y_{k}$. On the other hand, the
choice $\rho=10$ is clearly too large, which is to be expected from the fact
that $\sigma_{\text{max}}=3$. In this case, whilst after a few iterations
$x_{k}$ and $y_{k}$ become similar, the overall reconstruction error remains
very high. For $\rho=1$, the iterates $x_{k}$ and $y_{k}$ do seem to converge
to a satisfactory solution. Although the blending operator depends on the
firing times during acquisition, its sensitivity to $\rho$ is dictated by its
spectrum. Moreover, the noise level of the blending noise is likely to be
similar for different blending scenarios: the noise will always be of the same
order of magnitude as that of the signal. As the spectrum of different
blending operators will also be similar, it seems safe to conclude that
$\rho=1$ is a choice that is likely to work for our algorithm in general.
Additionally, although not applied in our numerical examples, the difference
between $x_{k}$ and $y_{k}$ would represent a good stopping criterion for our
algorithm.

#### The $y$-update

In our implementation we propose to start with a randomly initialized network
and train it for a fixed number of epochs at every outer iteration. A warm
start strategy is employed such that the weights of the network at a given
outer iteration are initialized to those obtained at the end of the training
of the previous outer iteration. The efficacy of this approach is shown in
figure 9a, where we compare on-the-fly training with and without warm starts,
where the latter re-initializes the network at every $y$-update. From the
error curves, we can safely conclude that warm starting the network is clearly
beneficial. Since there is no theoretical justification for this particular
strategy, we consider a few other alternative strategies. The first strategy
is to use a pre-trained network. Here pre-training is achieved by denoising
the pseudo-deblended data in a self-supervised manner; this approach could
greatly reduce the computational cost of the overall algorithm since we do not
need to train the network at every iteration. A comparison of the relative
error with that of the proposed, on-the-fly training shown in figure 9b
reveals that after a few outer iterations, the network is unable to further
remove the remaining noise in the data.

Figure 9: a) Error for network training with and without warm starts. b) Error
when running the PnP algorithm with a network pre-trained on the pseudo-
deblended data.

Another option is to stop training the network after a few outer iterations.
Ideally, the network will have learnt how to remove the noise encountered
during the first iterations, and extra training will not improve the denoiser
capabilities. We run experiments where we stop the training after a fixed
number of outer iterations to see whether there is an added benefit to
continuing training the network. Results are shown in figure 10a.

Figure 10: a) Error for on-the-fly training where training is stopped after a
certain number of outer iterations. b) Error for different number of training
epochs. c) Error when using the network at the end of the PnP algorithm for
the entire process versus our proposed on-the-fly training strategy.

In all of the scenarios we clearly see that stopping the training after a
certain number of outer iterations leads to a stagnation in the error, or even
worse to an increase in the error at later iterations. This behaviour is known
as semiconvergence in the inverse problems community. Both results are perhaps
not surprising, as the input to the network at every iteration contains a
different noise level compared to that of earlier iterations: the noise in
$x_{k}$ constantly reduces during the overall inversion. Therefore, the
network is required to learn a slightly different task at each time. In figure
10b, we assess the impact of the number training epochs for the denoiser. We
clearly observe that at some point the curves for 30 and 20 epochs start to
coincide, meaning that there is no additional gain in the extra 10 epochs of
training. The curve for training with 10 epochs seems stagnant after the first
three iterations, but eventually it picks up momentum and goes down again.
Note that, in terms of overall epochs, the cost of performing 30 outer
iterations with 10 epochs each is the same as using 10 outer iterations with
30 epochs each. However, every outer iteration carries an additional cost of
three inner iterations for the $x$-update, which is not negligible as it
requires evaluating the forward and adjoint of the blending operator. In
general, it seems beneficial to perform more epochs in the early outer
iterations, although there is a limit after which the error starts to
stagnate. Moreover, training with 60 epochs leads to overfitting.
To further investigate the generalization capabilities of the network for
blending problems, the network weights are saved after the last outer
iteration of the PnP algorithm. The PnP algorithm is then re-run using the
saved network without performing any on-the-fly training. Figure 10c shows
that this strategy fails, illustrating that the network may have forgotten how
to deal with the higher noise levels encountered in the early iterations. This
results highlights the importance of using a self-supervised denoiser that can
be easily and cheaply trained on-the-fly. The use of pre-trained denoising
networks such as DnCNN may instead require training multiple networks with
different noise levels, unless a bias-free, non-blind network is used [51].
Finally, in order to assess the influence of network initialization on the
final deblending results, we run our algorithm for 10 different random
initializations of the network weights and biases. We show the resulting error
curves as function of outer iterations in figure 11a. Moreover, figure 11b
displays a box plot of the final relative errors compared to that of the
conventional patched Fourier approach. We can observe that apart from one
seed, all the others tend to produce a final deblending result of superior
quality to the benchmark algorithm.

Figure 11: a) The average relative error plus and minus one standard deviation
as function of outer iterations. Min and max error curves are also displayed
based on the error after 30 iterations. b) Box plot of the relative error
after 30 outer iterations for 10 different seeds. A horizontal red line
indicates the error by the conventional patched Fourier approach.

##  6 Limitations and Conclusions

#### Limitations

Our algorithm requires the setting of a number of hyperparameters, namely the
number of inner and outer iterations, the parameter $\rho$, and the number of
epochs for the self-supervised denoiser. The number of epochs seems to have a
major impact on the quality of the deblending process and the overall
convergence properties of our algorithm. Perhaps an adaptive strategy may
reduce the computational cost without impacting the quality of the
reconstruction. Since we warm start each training process, one could
hypothesize that some of the shallow features learned at one iteration may be
still valid for the subsequent denoising task. Moreover, additional
hyperparameters that have not been explored in this work are associated with
the network itself, e.g. the number of layers, the activation function, batch
size, etc. This is also a direction for further research. Another possible
drawback of our approach is the lack of convergence guarantees. PnP methods
have been shown to converge under certain assumptions, one of them being that
the physical operator in the data misfit term, in our case the blending matrix
$B$, has to be strongly convex [38]. This requirement is never met for the
deblending problem, as $B$ is always underdetermined. Empirically, we observe
that $x_{k}$ and $y_{k}$ tend to converge to similar values for some carefully
selected hyperparameters, indicating that at least in our experiments the
algorithm converges successfully. Moreover, the distance between $x_{k}$ and
$y_{k}$ may be used as a suitable stopping criterion for the outer iterations,
which is another hyperparameter that is currently set upfront. Similarly, to
obtain convergence guarantees for the PnP method, the denoiser has to be
Lipschitz continuous; when a neural network is used, this means that spectral
normalization is required during training. In [49], it was shown that PnP
algorithms can be convergent when combined with carefully pre-trained
denoisers that satisfy such condition.

#### Societal impact

Blended acquisition greatly reduces the time required to acquire seismic data,
thereby limiting the impact of seismic acquisitions on the environment. Apart
from shooting at shorter intervals, there is no difference compared to
conventional acquisition. Moreover, since a larger number of shots can be
fired in the same overall acquisition time, recent research has suggested that
the energy emitted by each source could be lowered. This may provide
acquisition solutions that are more environmentally friendly for marine life.

#### Conclusions

We have introduced a novel hybrid algorithm for seismic deblending, combining
the physics of the blending operator with a self-supervised denoiser that is
naturally embedded into the Plug-and-Play framework. We have adapted the
network architecture in [28] to enforce an extended blind spot along an entire
axis (time, in our case) instead of single pixels. Because the denoiser is
self-supervised, our approaches bypasses the need for ground truth labels that
are usually unavailable for seismic applications. Experiments on a field
dataset have shown that the proposed method can outperform a state-of-the-art,
sparsity-based algorithm. Moreover, our algorithm is independent on the type
of acquisition, which is usually an issue for conventional algorithms.
Although our algorithm requires the setting of a number of hyperparameters, we
have argued that the number of inner iterations and $\rho$ can most likely be
set to a fixed number and this easily generalizes to different seismic
acquisitions. However, the network architecture and the number of epochs may
require tuning for different acquisition setups. We hope to address these
issues by having an adaptive strategy for setting the number epochs in future
work.

## References

  * [1] R. Abma, D. Howe, M. Foster, I. Ahmed, M. Tanis, Q. Zhang, A. Arogunmati, and G. Alexander.  Independent simultaneous source acquisition and processing.  Geophysics, 80(6):WD37–WD44, 2015.
  * [2] T. Alkhalifah, H. Wang, and O. Ovcharenko.  Mlreal: Bridging the gap between training on synthetic data and real data applications in machine learning.  arXiv, 2021.
  * [3] C. Bagaini.  Acquisition and processing of simultaneous vibroseis data.  Geophysical Prospecting, 58:81–99, 2010.
  * [4] B. Bahia, R. Rongzhi Lin, and M. Sacchi.  Regularization by denoising for simultaneous source separation.  Geophysics, 86(6):1942–2156, 2021.
  * [5] J. Batson and L. Royer.  Noise2self: Blind denoising by self-supervision.  International Conference on Machine Learning, pages 524–533, 2019\.
  * [6] C.J. Beasley, R.E. Chambers, and Z. Jiang.  A new look at simultaneous sources.  In 68th Meeting, SEG Expanded Abstracts,, volume 2022, pages 133–135. European Association of Geoscientists & Engineers, 1998.
  * [7] G. Beaudoin and A.A. Ross.  Field design and operation of a novel deepwater, wide-azimuth node seismic survey.  The Leading Edge, 26:385–544, 2007.
  * [8] A. Beck and M. Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM Journal on Imaging Sciences, 2:183–202, 2009.
  * [9] A. J. G. Berkhout.  Changing the mindset in seismic data acquisition.  The Leading Edge, 27:924–938, 2008.
  * [10] M. Berraki, S. Buizard, J. Ramirez, R.M. Elde, D. Eckert S.S. Roy, and J.-F. Synnevag.  Grane prm - from acquisition to interpretation in record time.  In 79th EAGE Conference and Exhibition 2022, volume 2017, pages 1–5. European Association of Geoscientists & Engineers, 2017.
  * [11] B. Biondi.  3D Seismic Imaging Seismology (2nd ed.).  SEG Books, 2006.
  * [12] C. Birnie and T. Alkhalifah.  Leveraging domain adaptation for efficient seismic denoising.  In Energy in Data Conference, Austin, Texas, 20–23 February 2022, pages 11–15. Energy in Data, 2022.
  * [13] C. Birnie, M. Ravasi, L. Sixiu, and T. Alkhalifah.  The potential of self-supervised networks for random noise suppression in seismic data.  Artificial Intelligence in Geosciences, 2:47–59, 2021.
  * [14] C. Broaddus, A. Krull, M. Weigert, U. Schmidt, and G. Myers.  Removing structured noise with self-supervised blind-spot networks.  2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI), pages 159–163, 2020.
  * [15] A. Buades, B. Coll, and J.-M. Morel.  A non-local algorithm for image denoising.  IEEE Conference on Computer Vision and Pattern Recognition, 2:60–65, 2005.
  * [16] P. Chatterjee and P. Milanfar.  Patch-based near-optimal image denoising.  Image Processing, IEEE Transactions on, 21(4):161–175, 2012.
  * [17] J. Cheng and M.D. Sacchi.  Separation and reconstruction of simultaneous source data via iterative rank reduction.  Geophysics, 80(4):V57–V66, 2015.
  * [18] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian.  Image denoising by sparse 3-d transform-domain collaborative filtering.  IEEE Transactions on Image Processing, 16:2080–2095, 2007.
  * [19] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian.  Image denoising by sparse 3-d transform-domain collaborative filtering.  IEEE Transactions on image processing, 16(8):2080–2095, 2007.
  * [20] N. Elad and M. Aharon.  Image denoising via sparse and redundant representations over learned dictionaries.  Image Processing, IEEE Transactions on, 15(12):3726–3745, 2006\.
  * [21] S. Gu, R. Timofte, and L. van Gool.  Integrating local and non-local denoiser priors for image restoration.  In International Conference on Pattern Recognition, pages 2923–2928, 2018.
  * [22] W. Huang, R. Wang, X. Gong, and Y. Chen.  Iterative deblending of simultaneous-source seismic data with structuring median constraint.  IEEE Geoscience and Remote Sensing Letters, 15(1):58–62, 2017.
  * [23] S. Huo, Y. Luo, and P.G. Kelamis.  Simultaneous sources separation via multidirectional vector-median filtering.  Geophysics, 77(4):123–131, 2012.
  * [24] S. Huo, Y. Luo, and P.G. Kelamis.  Deblending using a space-varying median filter.  Exploration Geophysics, 46(4):332–341, 2015.
  * [25] A. Ibrahim and M.D. Sacchi.  Simultaneous source separation using a robust radon transform.  Geophysics, 79:V1–V11, 2014.
  * [26] P.R.S. Johann, E.A. Thedy, F.A. Gomes, and M.C. Schinelli.  4d seismic in brazil: Experiences in reservoir monitoring.  Offshore Technology Conference, 2006.
  * [27] A. Krull, T.-O. Buchholz, and F. Jug.  Noise2void - learning denoising from single noisy images.  Proc. IEEE Conf. Comput. Vision Pattern Recognit., 2019.
  * [28] S. Laine, T. Karras, J. Lehtinen, and T. Aila.  High-quality self-supervised deep image denoising.  Advances in Neural Information Processing Systems, pages 6970–6980, 2019.
  * [29] J. Lehtinen, J. Munkberg, J. Hasselgren, S. Laine, T. Karras, M. Aittala, and T. Aila.  Noise2noise: Learning image restoration without clean data.  International Conference on Machine Learning (ICML), 2019.
  * [30] C. Li, C.C. Mosher, and Y. Ji.  An amplitude-preserving deblending approach for simultaneous sources.  Geophysics, 84(3):V185–V196, 2019.
  * [31] C. Li, C.C. Mosher, and Y. Ji.  Randomized marine acquisition with compressive sampling matrices.  Geophysical Prospecting, 60(3):648–662, 2019.
  * [32] S. Liu, C. Birnie, and T. Alkhalifah.  Coherent noise suppression via a self-supervised deep learning scheme.  In 83rd EAGE Conference and Exhibition 2022, volume 2022, pages 1–5. European Association of Geoscientists & Engineers, 2022.
  * [33] A. Mahdad, P. Doulgeris, and G. Blacquiere.  Separation of blended data by iterative estimation and subtraction of blending interference noise.  Geophysics, 76:Q9–Q17.
  * [34] S. Mandelli, V. Lipari, P. Bestagini, and S. Tubaro.  Interpolation and denoising of seismic data using convolutional neural networks.  arXiv, 2019.
  * [35] T. Meinhardt, M. Moller, C. Hazirbas, and D. Cremers.  Learning proximal operators: Using denoising networks for regularizing inverse imaging problems.  In Proceedings of the IEEE International Conference on Computer Vision, pages 1781–1790, 2017.
  * [36] I. Moore, B. Dragoset, T. Ommundsen, D. Wilson, C. Ward, and D. Eke.  Simultaneous source separation using dithered sources.  SEG Technical Program Expanded Abstracts, 2008.
  * [37] Y. Romano, M. Elad, and P. Milanfar.  The little engine that could: Regularization by denoising (RED).  SIAM Journal on Imaging Sciences, 10(4):1804–1844, 2017.
  * [38] Ernest Ryu, Jialin Liu, Sicheng Wang, Xiaohan Chen, Zhangyang Wang, and Wotao Yin.  Plug-and-play methods provably converge with properly trained denoisers.  In International Conference on Machine Learning, pages 5546–5557. PMLR, 2019.
  * [39] Yousef Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [40] R. E. Sheriff and L.P. Geldart.  Exploration Seismology (2nd ed.).  Cambridge University Press, 1995.
  * [41] A. Siahkoohi, M. Louboutin, and F. J. Herrmann.  The importance of transfer learning in seismic modeling and imaging.  Geophysics, 84(6):A47–A52, 2019.
  * [42] J. Sun, S. Slang, T. Elboth, T.L. Greiner, S. McDonald, and L.-J. Gelius.  A convolutional neural network approach to deblending seismic data.  Geophysics, 85(4):WA13–WA26, 2020.
  * [43] Y. Sun, J. Liu, and U. Kamilov.  Block coordinate regularization by denoising.  In Advances in Neural Information Processing Systems, pages 380–390, 2019.
  * [44] J.-B. Thibault, K.D. Sauer, C.A. Bouman, and J. Hsieh.  A three-dimensional statistical approach to improved image quality for multislice helical ct.
  * [45] T. Tirer and R. Giryes.  Image restoration by iterative denoising and backward projections.  IEEE Transactions on Image Processing, 28(3).
  * [46] T. Tirer and R. Giryes.  Super-resolution via image-adapted de- noising CNNs: Incorporating external and internal learning.  IEEE Signal Processing Letters, 26(7).
  * [47] S.V Venkatakrishnan, C.A. Bouman, and B. Wohlberg.  Plug-and-play priors for model based reconstruction.  2013 IEEE Global Conference on Signal and Information Processing, 2013.
  * [48] S. Wang, W. Hu, P. Yuan, X. Wu, Q. Zhang, P. Nadukandi, G.O. Botero, and J. Chen.  Seismic deblending by self-supervised deep learning with a blind-trace network.  SEG/AAPG/SEPM First International Meeting for Applied Geoscience & Energy, 2021.
  * [49] Xiajian Xu, Yu Sun, Jiaming Liu, Brendt Wohlberg, and Ulugbek S. Kamilov.  Provable convergence of plug-and-play priors with mme denoisers.  IEEE Signal Processing Letters, pages 1280–1284, 2020.
  * [50] Siwei Yu and Jianwei Ma.  Deep learning for geophysics: Current and future trends.  Reviews of Geophysics, 59(3):e2021RG000742, 2021.
  * [51] K. Zhang, Y. Li, W. Zuo, L. Zhang, L. van Gool, and R. Timofte.  Plug-and-play image restoration with deep denoiser prior.  In IEEE TPAMI, 2021.
  * [52] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang.  Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising.  IEEE Transactions on Image Processing, pages 3142–3155, 2017.
  * [53] K. Zhang, W. Zuo, S. Gu, and L. Zhang.  Learning deep cnn denoiser prior for image restoration.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 3929–3938, 2017.
  * [54] H. Zhou, W. Mao, D. Zhang, Q. Ge, and H. Wang.  Deblending of simultaneous source with rank-reduction and thresholding constraints.  79th EAGE Conference and Exhibition 2017, 2017(1):1–5, 2017.
  * [55] Y. Zhou, W. Chen, and J. Gao.  Separation of seismic blended data by sparse inversion over dictionary learning.  Journal of Applied Geophysics, 106:146–153, 2014.
  * [56] Y. Zhou, J. Gao, W. Chen, and P. Frossard.  Seismic simultaneous source separation via patchwise sparse representation.  IEEE Transactions on Geoscience and Remote Sensing, 54(9):5271–5284, 2016.
  * [57] S. Zu, J. Cao, S. Qu, and Y. Chen.  Iterative deblending for simultaneous source data using the deep neural network.  Geophysics, 85(2):V131–V141, 2020.
  * [58] S. Zu, H. Zhou, R. Wu, W. Mao, and Y. Chen.  Hybrid-sparsity constrained dictionary learning for iterative deblending of extremely noisy simultaneous-source data.  IEEE Transactions on Geoscience and Remote Sensing, 57(4):2249–2262, 2018.

##  Appendix A Denoiser progression over iterations

Figure 12: Progression of the denoiser for outer iteration 1. Figure 13:
Progression of the denoiser for outer iteration 10. Figure 14: Progression of
the denoiser for outer iteration 20. Figure 15: Progression of the denoiser
for outer iteration 30.


# Minimal design of the elephant trunk as an active filament

Bartosz Kaczmarski, Sophie Leanza, Renee Zhao  Department of Mechanical
Engineering, Stanford University, Stanford, CA-94305, United States  Derek E.
Moulton  Mathematical Institute, University of Oxford, Oxford, UK  Ellen Kuhl
Department of Mechanical Engineering, Stanford University, Stanford, CA-94305,
United States  Alain Goriely  Mathematical Institute, University of Oxford,
Oxford, UK

###### Abstract

One of the key problems in active materials is the control of shape through
actuation. A fascinating example of such control is the elephant trunk, a
long, muscular, and extremely dexterous organ with multiple vital functions.
The elephant trunk is an object of fascination for biologists, physicists, and
children alike. Its versatility relies on the intricate interplay of multiple
unique physical mechanisms and biological design principles. Here we explore
these principles using the theory of active filaments and build,
theoretically, computationally, and experimentally, a minimal model that
explains and accomplishes some of the spectacular features of the elephant
trunk.

Introduction.—In nature and engineering, one of the main tasks of active
filamentary structures is to explore their surroundings. The simplest problem
for the activation of slender structures is curvature generation, a well-
understood phenomenon studied by Timoshenko a century ago in the context of
bi-metallic strips Timoshenko (1925). In Timoshenko’s initial setting,
curvature changes only generate planar shapes. A more involved problem is the
generation of torsion which, through curvature coupling, allows for arbitrary
three-dimensional shapes as found in plants Moulton  _et al._ (2020a), but
also in the octopus arms and elephant trunks Dagenais  _et al._ (2021). Here
we take inspiration from the elephant trunk as a design paradigm to couple
curvature, twist, and torsion, and create a minimal activation model using
three uniform actuators.

The elephant trunk is an elongated, muscular, and highly flexible proboscis
made of specialized muscles, fascia, and skin; it represents a unique
adaptation of the upper lip and the nose, see Fig. 1. It is made of 17 major
muscle groups, eight on each side and one for the nasal cavity, with more than
90,000 muscle fascicles controlled by 50-60,000 facial nucleus neurons, which
together give it its exceptional strength and versatility of movement Kaufmann
_et al._ (2022); Purkart  _et al._ (2022); Longren  _et al._ (2023).

Figure 1: Muscular architecture within the elephant trunk. (a) Skinned
elephant trunk with longitudinal muscles on the dorsal side and oblique
muscles wrapping helically along the length on the ventral side Boas and
Paulli (1908). (b) Cross ection of a trunk by Cuvier (1850) with longitudinal
(blue), oblique (red, green) and radial muscles around the two nasal cavities
Cuvier and Laurillard (1850). (c) Minimal model of a trunk with three muscular
bundles.

Unsurprisingly, the biomechanics of the elephant trunk have attracted the
attention of many scientists Wu  _et al._ (2018); Schulz  _et al._ (2022) with
the ultimate goal to understand and replicate some of its functions Guan  _et
al._ (2020). Yet, some of the basic principles of the trunk have been
overlooked. Here, in contrast to most previous studies that use a great number
of segmented actuators to replicate function Trivedi  _et al._ (2008); Zhang
_et al._ (2023), we aim to unravel important properties of the trunk by
considering an ideal trunk with a minimum number of actuators that still
reproduces some of its key functions. Towards this goal, we show that simple
physical models of active materials can help improve our understanding of
these spectacular structures. Inspired by the muscular architecture of the
trunk, we propose a minimal model of three actuators and investigate its
control capabilities. In particular, we study the design space of our
idealized trunk and identify the fundamental trade-off between some desirable
properties. Then, we validate our model predictions experimentally using
liquid crystal elastomer actuators.

Model set-up.—Despite its massive size, the elephant trunk can be modeled as a
slender biological filament that deforms by muscular activation of
longitudinal, radial, and helical muscles (Fig. 1). To model the deformation
induced by local fiber contraction, we adopt the active filament theory
Kaczmarski  _et al._ (2023a, b, 2022), which describes a general three-
dimensional tubular structure equipped with an arbitrary fiber activation
tensor $\mathbf{G}$ dictated by the local orientation of the muscle fibers
throughout the continuum. Assuming that initially the trunk is straight, any
deformation can be written as $\bm{\chi}=\bm{r}+\sum_{i=1}^{3}\varepsilon
e_{i}\mathbf{d}_{i}$ Moulton  _et al._ (2020b), where $\bm{r}(Z)$ is the
centerline of the deformed tubular structure w.r.t. the arc length $Z$ of the
initial structure, such that $\bm{r}^{\prime}=\zeta\mathbf{d}_{3}$, where
$\zeta$ is the centerline stretch, $\varepsilon e_{i}$ are the cross-sectional
reactive strains with $\varepsilon$ the slenderness ratio (a characteristic
radius divided by the length $L$), and
$\\{\mathbf{d}_{1},\mathbf{d}_{2},\mathbf{d}_{3}\\}$ form an orthonormal basis
that defines the orientation of each cross-section so that
$\mathbf{d}_{i}^{\prime}=\zeta\mathbf{u}\times\mathbf{d}_{i}$, where
$\mathbf{u}=\sum{\mathsf{u}}_{i}\mathbf{d}_{i}$ is the Darboux vector
specifying the filament’s three curvatures. To model activation, we decompose
the deformation gradient Goriely (2017),
$\mathbf{F}=\text{Grad}{\bm{\chi}}=\mathbf{A}\cdot\mathbf{G}$, into an elastic
contribution $\mathbf{A}$ and an activation $\mathbf{G}$. Using a systematic
energy minimization procedure, we then obtain the intrinsic curvatures
$\hat{{\mathsf{u}}}_{1}$, $\hat{{\mathsf{u}}}_{2}$, $\hat{{\mathsf{u}}}_{3}$,
and the extension $\hat{\zeta}$ of the trunk in the absence of external loads
by integrating the contribution of each fiber contraction over the cross-
section Kaczmarski  _et al._ (2022). In the case of a ring of helical fibers
with helical angle $\alpha$ and radii $R_{1}$ and $R_{0}$ (see Fig. 2), the
curvatures are:

$\displaystyle\hat{{\mathsf{u}}}_{1}(Z)$
$\displaystyle=-\frac{4}{3R_{0}^{4}}A\delta_{1}\sin\left(\varphi-\frac{Z}{R_{0}}\tan\alpha\right),$
$\displaystyle\hat{{\mathsf{u}}}_{2}(Z)$
$\displaystyle=-\frac{4}{3R_{0}^{4}}A\delta_{2}\cos\left(\varphi-\frac{Z}{R_{0}}\tan\alpha\right),$
$\displaystyle\hat{{\mathsf{u}}}_{3}(Z)$
$\displaystyle=\frac{2}{R_{0}^{4}}\delta_{3}a_{0},\quad\hat{\zeta}=1+\frac{1}{2R_{0}^{2}}\delta_{0}a_{0},$

where $\delta_{i}$ are functions of the filament geometry, fiber architecture,
and mechanical properties (see Supplementary Material). For a given periodic
muscle activation function $\gamma=\gamma(\theta)$, $A$, $\varphi$, $a_{0}$
are

$\displaystyle a_{1}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{2\pi}\\!\\!\\!\\!\gamma(\theta)\cos\theta\,\text{d}\theta,\quad
b_{1}=\frac{1}{\pi}\int_{0}^{2\pi}\\!\\!\\!\\!\gamma(\theta)\sin\theta\,\text{d}\theta,$
$\displaystyle a_{0}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{2\pi}\\!\\!\\!\\!\gamma(\theta)\,\text{d}\theta,\quad
A=\sqrt{a_{1}^{2}+b_{1}^{2}},\quad\tan(\varphi)=-\frac{b_{1}}{a_{1}}.$

We restrict our attention to discrete muscle fiber bundles with inputs
$\gamma_{1},\ldots,\gamma_{N}$ and angular extents defined by the parameter
$\sigma$ as shown in Fig. 2ce.

Figure 2: Active filament model. (a) Filament with embedded helically-arranged
fibers. (b) and (d) Realistic continuous muscle activation $\gamma(\theta)$.
(c) and (e) Idealized discrete piecewise constant activation with $N=3$
actuators. Figure 3: Design space of the elephant trunk: slenderness,
tapering, and fiber helicity. Deformations I, II, III, IV and their respective
curvature-activation curves for: (a) varied slenderness
($S_{\text{r}}=5,10,15,20$, $\gamma=-0.6$) with a single longitudinal bundle
($\sigma=90^{\circ}$); (b) varied tapering angle
($\phi=1^{\circ},3^{\circ},5^{\circ},7^{\circ}$, $\gamma=-2.5$) for the same
single bundle; and (c) varied helicity
($\Omega\approx{10.0^{\circ},29.2^{\circ},48.3^{\circ},67.5^{\circ}}$,
$\gamma=-3$, $\sigma=45^{\circ}$ in both bundles) for two symmetrically
arranged, equally activated helical fiber bundles. Dashed regions highlight
physiological ranges $S_{\text{r}}\in[6.7,11.7]$, and
$\phi\in[5^{\circ},6^{\circ}]$ for the elephant trunk, and for $\Omega$ far
from the longitudinal case. Young’s modulus is constant, Poisson’s ratio is
$\nu=1/2$, and $R_{0}-R_{1}=L/50$ at $Z=0$.

Design principles.—We apply the active filament theory to study the main
design principles that govern the configuration space of the elephant trunk:
slenderness, tapering, and helicity. First, we consider a single longitudinal
muscular bundle. Upon activation, the filament curls to minimize its elastic
energy with dimensionless curvature $L\kappa$, extending Timoshenko’s result
to large deformations. In Fig. 3a, we show that slenderness has a significant
effect on the deformation. Indeed, by varying the ratio $S_{\text{r}}$ of
filament length $L$ to the outer radius $R_{0}$, the same muscular activation
will yield a larger curvature for filaments with a larger slenderness. This is
a natural consequence of the scaling of the second moment of area of the
filament’s cross-section as $\sim 1/S_{\text{r}}^{4}$, which dominates the
less pronounced effect of the decreased fiber activation area for more slender
filaments. Fig. 3a further shows that, for sufficiently small fiber magnitudes
($|\gamma|<1$), the dimensionless curvature $L\kappa$ is an approximately
linear function of the activation in the entirety of the evaluated slenderness
range $S_{\text{r}}\in[5,20]$.

Second, the geometry of the elephant trunk is highly tapered, with a larger
cross-sectional radius at the proximal end and increasingly smaller radii
towards the distal end, effectively inducing a gradually varying slenderness.
For uniform activations in a longitudinal muscle, the local curvature of the
trunk increases with the decreasing trunk radius. The larger the decrease in
radius the larger the increase in curvature, as demonstrated in Fig. 3b by the
increase in maximal curvature $L\kappa_{\text{max}}$ with the increasing
tapering angle $\phi$. In the physiological range of tapering angles
$\phi\in[5^{\circ},6^{\circ}]$ (structures III and IV), we observe the typical
curling motion that the elephants use to wrap their trunks around objects.

Third, the elephant trunk contains oblique muscle fibers that we model through
two opposite helical fibers. In addition to creating curvature, a single
helical fiber bundle can generate torsion. Combining the effect of two
identical helical fibers but of opposite handedness cancels the torsional
contribution and creates pure bending of the filament Goriely and Tabor (2013)
as shown Fig. 3c. An increase in the fiber revolution $\Omega$ results in
larger maximal curvature and increased variability in the bending curvature
(see Supplementary Material for the definition of the angle $\Omega$). The
influence of $\Omega$ on the resulting deformation is similar to the effect of
varying the tapering angle $\phi$: the rate of increase of curvature along $Z$
is amplified by increasing either $\Omega$ or $\phi$, which might be
responsible for generating the curling deformations commonly observed in
elephant trunks motions.

Figure 4: Reachability cloud for a minimal elephant trunk design with three
actuators. (a) The point cloud consists of 2 million points, RGB-coded by
activation magnitude within the range
$\gamma_{1},\gamma_{2},\gamma_{3}\in[-5,0]$ (red, green, and blue,
respectively). (b) A design with only 2 helical actuators (shown in red)
restricts the reachability cloud to a 2D surface region compared to 3D cloud
of the minimal design (in blue). (c) Similarly, a design with 3 longitudinal
actuators has a small reachability cloud with a large interior unreachable
region.

Minimal elephant trunk-inspired design.—Our study in Fig. 3 suggests that a
single muscular bundle or two equal opposite helical bundles can generate a
variety of deformation modes. By combining these modes, we can create a
minimal design of a tubular structure made up of one longitudinal and two
helical fibers. The longitudinal fiber generates pure curvature as inspired by
the massive upper muscle group in the elephant trunk, see Fig. 1. The two
symmetrically arranged helical fibers generate either twist, torsion, or
curvature depending on their relative activation magnitudes; see Fig. 4 and
the design parameters in the Supplementary Material.

To quantify the performance of this design, we created a comprehensive three-
dimensional reachability point cloud in Fig. 4a to illustrate the space of end
points accessible by the trunk by randomly sampling muscular activations.
Notably, the reachability cloud consists of 2 million points and takes less
than 1 minute to compute on a standard desktop computer, at a rate of roughly
100,000 configurations per second thanks to the semi-analytical nature of the
model and unlike traditional finite-element code.

The resulting cloud spans an extensive three-dimensional space within the
entire $360^{\circ}$ arc around the longitudinal axis of the initial
configuration. Our minimal trunk design can reach points to the front, back,
right, and left of the initial configuration, as well as anywhere in between
these regions. For comparison, the trunk of an Asian elephant can generate
compressive strains up to 33% Wilson  _et al._ , which corresponds to a
stretch of 0.67 and an activation of $\gamma\approx-3.33$ in our minimal
design. Importantly, a tight concave hull of the reachability cloud with
normalized volume $\tilde{V}_{\text{cloud}}=V_{\text{cloud}}/L^{3}\approx
0.85$ amounts to about $74\%$ of the convex hull with
$\tilde{V}_{\text{conv}}\approx 1.14$, which implies that only $26\%$ of the
conservative volume bounded by the cloud is inaccessible by the design.
Comparison with other designs containing either 2 helical actuators or 3
longitudinal ones shows in Fig. 4bc that our minimal design is vastly superior
in terms of reachability. Specifically, removing the blue longitudinal
actuator reduces the reachability cloud to a 2D surface (Fig. 4b), and
replacing the helical actuators with longitudinal ones decreases the cloud
volume by $14\%$ to $\tilde{V}_{\text{cloud}}\approx 0.73$ due to a large
inaccessible void created by the removal of helical fibers, amounting to
almost $42\%$ of the new convex hull (Fig. 4c). Indeed, the torsional modes
activated by the helical fibers empower the design to not only reach a larger
space, but also perform a variety of complex motions, including avoiding
obstacles or grasping and rotating objects through a curling motion, a common
motion performed by elephant trunks when manipulating branches of trees.

Figure 5: Experimental validation of predicted deformation patterns.
Activation of the single longitudinal, a single helical, both helical, and a
pair of longitudinal and helical fibers, from left to right. The predicted
deformation patterns, top, agree with the experimental deformation modes of a
soft slender structure activated by liquid-crystal elastomer fibers, bottom.

Experimental validation.—To conclude our study, we validate our predicted
deformation patterns with experimental deformation patterns of a soft slender
structure actuated by liquid crystal elastomer fibers, which contract upon
heating. Briefly, we created an inactive acrylate polymer-based slender
cylinder via molding. Three active liquid crystal elastomer fibers were
fabricated through 3D printing for mesogen alignment. Before curing, we
embedded a coiled copper wire in each of the fibers to enable active
contraction through Joule heating. We then cured the fibers with ultraviolet
light to fix the mesogen alignment and lock the wire in place. Finally, we
bonded the resulting fibers to the elastic cylinder, giving us the complete
experimental model.

Fig. 5 shows a direct comparison of four representative modes from the
activation of (i) the single longitudinal fiber, (ii) a single helical fiber,
(iii) both helical fibers, and (iv) a pair of longitudinal and helical fibers,
from left to right. For a more accurate comparison with the experiment, the
model predictions in Fig. 5 include the effects of gravity (see Supplementary
Material). The four configurations in Fig. 5 illustrate the richness of
possible deformation modes: configuration (i) exhibits pure bending achieved
by activating the longitudinal actuator only; configuration (ii) exhibits high
torsion through the activation of one of the helical actuators; configuration
(iii) exhibits pure bending through the two equal and opposite helical
actuators, with the bending direction changing halfway following the fiber
arrangement; finally, configuration (iv) represents a combination of bending
and torsion, generating an out-of-plane motion induced by the interplay of the
longitudinal and helical actuators. The simulated configurations agree with
the experimental deformation patterns, which supports the utility of our
parametric studies and reachability considerations.

Conclusion.—Unlike traditional robotic arms with finite degrees of freedom,
soft structures like the elephant trunk or the octopus arm have infinitely
many, hence creating “a body of pure possibility” according to the philosopher
Peter Godfrey Godfrey (2016). This universe of possibilities is made possible
through fundamental physical mechanisms coupling geometry to mechanics,
amplified by the slenderness of the structure and the softness of matter.
While curvature generation in a one-dimensional structure is easy to
understand and create through longitudinal contraction, the creation of twist
and torsion requires dedicated helical fibers. To create both left- and right-
handed structures, the design requires two helical bundles with opposite
handedness. Hence, the simplest design is a three-actuator architecture,
inspired by the elephant trunk. Here we demonstrate, analytically,
computationally, and experimentally, that this minimal design is sufficient to
reach a large portion of space. Further, the large range of deformation types
empowered by this minimal design suggests that it could play an important role
in motion planning tasks.

###### Acknowledgements.

Funding: This work has been funded by EPSRC grant EP R020205 and the NSF grant
CMMI 2318188.
Authors contribution: BK, EK, DM, AG designed the study and developed the
theory, BK performed the computations, SL, RZ designed the actuators and
performed the experiments. All authors analyzed the data, discussed the
results, and wrote the paper.

## References

  * Timoshenko (1925) S. P. Timoshenko, J. Opt. Soc. Am. 11, 233 (1925).
  * Moulton  _et al._ (2020a) D. E. Moulton, H. Oliveri, and A. Goriely, Proceedings of the National Academy of Sciences 117, 32226 (2020a).
  * Kaufmann  _et al._ (2022) L. V. Kaufmann, U. Schneeweiß, E. Maier, T. Hildebrandt, and M. Brecht, Science Advances 8, eabq2789 (2022).
  * Purkart  _et al._ (2022) L. Purkart, J. M. Tuff, M. Shah, L. V. Kaufmann, C. Altringer, E. Maier, U. Schneeweiß, E. Tunckol, L. Eigen, S. Holtze, _et al._ , Current Biology 32, 904 (2022).
  * Longren  _et al._ (2023) L. L. Longren, L. Eigen, A. Shubitidze, O. Lieschnegg, D. Baum, J. A. Nyakatura, T. Hildebrandt, and M. Brecht, Current Biology (2023).
  * Dagenais  _et al._ (2021) P. Dagenais, S. Hensman, V. Haechler, and M. C. Milinkovitch, Current Biology 31, 4727 (2021).
  * Boas and Paulli (1908) J. E. V. Boas and S. Paulli, (1908).
  * Cuvier and Laurillard (1850) G. Cuvier and C. L. Laurillard,  _Anatomie comparée, recueil de planches de myologie dessinées par Cuvier… ou Laurillard, publié sous les auspices de ministre de l’instruction par Laurillard et Mercier:…_ (Dusacy, 1850).
  * Wu  _et al._ (2018) J. Wu, Y. Zhao, Y. Zhang, D. Shumate, S. Braccini Slade, S. V. Franklin, and D. L. Hu, Journal of the Royal Society Interface 15, 20180377 (2018).
  * Schulz  _et al._ (2022) A. K. Schulz, M. Boyle, C. Boyle, S. Sordilla, C. Rincon, S. Hooper, C. Aubuchon, J. S. Reidenberg, C. Higgins, and D. L. Hu, Proceedings of the National Academy of Sciences 119, e2122563119 (2022).
  * Guan  _et al._ (2020) Q. Guan, J. Sun, Y. Liu, N. M. Wereley, and J. Leng, Soft robotics 7, 597 (2020).
  * Trivedi  _et al._ (2008) D. Trivedi, C. D. Rahn, W. M. Kier, and I. D. Walker, Applied bionics and biomechanics 5, 99 (2008).
  * Zhang  _et al._ (2023) J. Zhang, Y. Li, Z. Kan, Q. Yuan, H. Rajabi, Z. Wu, H. Peng, and J. Wu, Soft Robotics 10, 636 (2023).
  * Kaczmarski  _et al._ (2023a) B. Kaczmarski, D. E. Moulton, A. Goriely, and E. Kuhl, Computer Methods in Applied Mechanics and Engineering 408, 115939 (2023a).
  * Kaczmarski  _et al._ (2023b) B. Kaczmarski, A. Goriely, E. Kuhl, and D. E. Moulton, IEEE Robotics and Automation Letters 8, 936 (2023b).
  * Kaczmarski  _et al._ (2022) B. Kaczmarski, D. E. Moulton, E. Kuhl, and A. Goriely, Journal of the Mechanics and Physics of Solids 164, 104918 (2022).
  * Moulton  _et al._ (2020b) D. E. Moulton, T. Lessinnes, and A. Goriely, Journal of the Mechanics and Physics of Solids , 104022 (2020b).
  * Goriely (2017) A. Goriely,  _The Mathematics and Mechanics of Biological Growth_ (Springer Verlag, New York, 2017).
  * Goriely and Tabor (2013) A. Goriely and M. Tabor, Proc. Roy. Soc. Lond. A 469 (2013).
  * (20) J. F. Wilson, U. Mahajan, S. A. Wainwright, and L. J. Croner,  113, 79.
  * Godfrey (2016) S. Godfrey, William Collins, Glasgow (2016).