---
title: Time-jittered marine acquisition --- low-rank v/s sparsity
author: |
    R. Kumar``^1``, H. Wason``^{1}`` and F. J. Herrmann``^1``\
    ``^{1}`` Seismic Laboratory for Imaging and Modeling (SLIM), University of British Columbia, Canada\
bibliography:
	- icml.bib   
---


## Abstract:

Time-jittered simultaneous marine acquisition has been recognized as an economic way of improving the spatial sampling, and speedup acquisition, where a single (or multiple) source vessel fires at -- jittered source locations and instances in time. It has been shown in the past that this problem can be setup as a -- compressed sensing problem, where conventional seismic data is reconstructed from blended data via a sparsity-promoting optimization formulation. While the recovery quality of deblended data is very good, the recovery process is computationally very expensive. In this paper, we present a computationally efficient rank-minimization algorithm to deblend the seismic data. The proposed algorithm is suitable for large-scale seismic data, since it avoids SVD computations and uses a low-rank factorized formulation instead. Results are illustrated with simulations of time-jittered marine acquisition, which translates to jittered source locations for a given speed of the source vessel, for a single source vessel with two airgun arrays.


## Introduction 
Simultaneous (or blended) marine acquisition mitigates the challenges posed by conventional marine acquisition in terms of sampling and make large area acquisition \emph{relatively} inexpensive when compared to conventional acquisition [@beasley98ana;@deKok2002eage;@Berkhout2008cms;@Beasley2008anl;@hampson2008aus]. While the final objective of deblending is to get the interference-free shot gathers, the recovery of late reflections is challenging as they can be overlaid by interfering seismic responses from other shots. [@Mansour11TRssma, @wason2013EAGEobs] have shown that this challenge can be effectively addressed through a combination of tailored single- (or multiple-) source, blended acquisition design and curvelet-based sparsity-promoting recovery. The idea was to design a pragmatic time-jittered marine acquisition scheme where acquisition related costs are no longer determined by the Nyquist sampling criteria, but by the transform-domain sparsity of the data. While the proposed sparsity-promoting approach to deblend the data is very effective, it poses computational challenges since the curvelet-based sparsity-promoting methods are not computationally tractable (both in terms of speed and memory usage) for large-scale 5D seismic data volumes. We refer to deblending as interpolating the sub-Nyquist jittered shot positions to a fine regular grid while unraveling the overlapping shots. To avoid this, we propose rank-minimization based methods to deblend the time-jittered simultaneous data volumes. More recently, rank-minimization based techniques have been used successfully for
source separation [@maraschini2012EAGEssi; @cheng2013SEGsss; @wason2014SEGsss]. The general idea is to exploit the low-rank structure of seismic data when it is organized in a matrix, in some transform domain. Low-rank structure refers to a small number of nonzero singular values, or quickly decaying singular values. 
Following [@aravkin2014SISCfmd], in time-jittered simultaneous marine acquisition acquisition, we observe that monochromatic frequency slices of the fully sampled data matrix (conventional acquisition) have low-rank structure in the midpoint-offset (m-h) domain as compared to source-receiver (s-r) domain, whereas, blending increases the rank of the resulting frequency slice in the m-h domain. Hence, the blended data in time-jittered simultaneous marine acquisition is separated (and in some cases interpolated) into its constituent source components using a fast rank-minimization approach that combines the (SVD-free) matrix factorization approach recently developed by [@Lee2010] with the Pareto curve approach proposed by [@BergFriedlander:2008]. The experimental results demonstrate the successful implementation of the proposed methodology. In addition, we also make comparisons with the sparsity-promoting based deblending techniques.

##Low-rank extension
 Let ``X_0`` be a matrix in ``\mathbb{C}^{n_m\times n_h \times n_f}`` and ``\mathcal{A}`` be a linear measurement operator that maps from
  ``\mathbb{C}^{n_m\times n_h \times n_f} \rightarrow \mathbb{C}^{n_r\times n_s \times n_f}``, where ``n_r`` represents number of receivers, ``n_s`` represents number of sources,
  ``n_m`` is the number of midpoints, ``n_h`` is the number of offsets and ``n_f`` is the number of frequencies.  
  [@RechtFazelParrilo2010] showed that under
  certain general conditions on the operator ``\mathcal{A}``, the
  solution to the rank-minimization problem can be found by solving
  the following nuclear-norm minimization problem:
```math #BPDN
\min_{X} ||X||_*  \quad \text{s.t. } \|\mathcal{A} (X) - b\|_2 \leq \epsilon, \tag{BPDN$_\epsilon$}
```
  where ``b`` is a set of measurements, ``\left\| X\right\|_* =
  \|\sigma\|_1``, and ``\sigma`` is the vector of singular values. The sampling-transformation operator ``\mathcal{A}`` is composed of the product of 
  a restricted-measurement operator ``RM`` and the transformation operator ``\mathcal{S}`` such that ``{\mathcal{A}} = RM\mathcal{S}^H``,
  where ``^H`` denotes the Hermitian transpose. In the case of seismic data interpolation, ``RM`` is separable into a ``n \times N`` restriction matrix (``R``) and a ``N \times N``
  measurment matrix ``M``, where ``N = n_t \times n_r \times n_s`` and ``n = n_t \times n_r \times n_{sub}``. Therefore, [#BPDN] is solved in the frequency domain where each frequency slice is reconstructed independently. 
  However, in the time-jittered simultaneous marine acquisition, a single- (or multiple-) source vessel maps the entire survey area while firing sequential shots at randomly time-dithered instances[@Mansour11TRssma] resulting in
 a single long \emph{supershot} of length ``n \leq N``. As a result, ``RM`` turns into a non-separable operator, thus, we can not 
independently solve over different monochromatic frequency slices. Therefore, we redefine our nuclear-norm minimization formulation as
```math
\min_{X} \sum_{i}^{n_f} ||X_i||_*  \quad \text{s.t. } \|\mathcal{A} (X) - b\|_2 \leq \epsilon,
```
where ``X_i \in \mathbb{C}^{n_m\times n_h} `` represents a monochromatic slice of unblended seismic data, ``b`` is the blended data, index ``{i}`` runs over frequencies,
``{\mathcal{A}}= {RMF^H}\mathcal{S}^H``, ``F`` is the Fourier transform-operator along the time-axis and ``M`` is the blending operator. In order to efficiently solve [#BPDN], we use an extension of
  the SPG``\ell_1`` solver [@BergFriedlander:2008] developed for
  the [#BPDN] problem in [@AravkinBurkeFriedlander:2012]. The
  SPG``\ell_1`` algorithm finds the solution to the [#BPDN] problem by
  solving a sequence of LASSO subproblems:

```math #LASSO
\min_{X} \|\mathcal{A} (X) - b\|_2 \quad \text{s.t. } \sum_{i}^{n_f}||X_i||_* \leq \tau, \tag{LASSO$_\tau$}
``` 
  where ``\tau`` is updated by traversing the Pareto curve. Solving each
  LASSO subproblem requires a projection onto the nuclear norm ball
  ``{\sum_{i}^{n_f}} \left\| X_i\right\|_* \leq \tau`` in every iteration by performing a
  singular value decomposition on each monochromatic slice and then thresholding the singular
  values. In the case of large-scale seismic problems, it becomes
  prohibitively expensive to carry out such a large number of
  SVDs. Instead, we adopt a recent factorization-based approach to
  nuclear-norm minimization [@Rennie2005;@Lee2010;@RechtRe:2011].
  The factorization approach parametrizes the matrix ``X``
  ``\in \mathbb{C}^{n_h\times n_m \times n_f}`` as the product of two low-rank factors
  ``L \in \mathbb{C}^{n_m\times k \times n_f}`` and ``R \in \mathbb{C}^{n_h\times k \times n_f}`` such that,
``X = {LR^H}``.
   The optimization scheme can then be carried out using the factors
   ``L, R`` instead of ``X``, thereby
   significantly reducing the size of the decision variable from ``2n_mn_hn_f``
   to ``2k(n_m +n_h)n_f`` when ``k \ll n_m,n_h``. Following [@Rennie2005], the
   nuclear norm obeys the relationship:

```math #nucInequality
{\sum_{i}^{n_f}}\|X\|_* \leq  {\sum_{i}^{n_f}}\frac{1}{2}\left\|L_i; R_i\right\|_F^2 =: {\sum_{i}^{n_f}}\Phi(L_i,R_i),
```
where ``\|\cdot\|_F^2`` is the Frobenius norm of the matrix (sum of the squared entries). Consequently, the LASSO subproblem can be replaced by 

```math #subproblem
\min_{L,R}  \|\mathcal{A}(X) - b\|_2 \quad \text{s.t.} \quad {\sum_{i}^{n_f}}\Phi(L_i,R_i) \leq \tau\;, 
```
  where the projection onto ``\sum_{i}^{n_f}\Phi(L_i,R_i) \leq \tau`` is
  easily achieved by multiplying each factor (``L_i, R_i``) by the scalar ``\sqrt{2\tau/\sum_{i}^{n_f}\Phi(L_i,R_i)}``. 
  
 
##Results
In order to show the efficacy of the proposed formulation, we use two realistic field data sets from the Gulf of Suez,
where one data set is sampled at the source (and receiver) sampling of 25.0m and the other data set is sampled at the source
(and receiver) sampling of 12.5m, with ``n_s = 128`` shots,
``n_r = 128`` receivers and ``n_t = 1024`` time samples. The acquisition involves a single source vessel with two airgun arrays 
firing at a 50.0m jittered grid with the receivers (OBC) recording continuously. Detailed explanation of the acquisition design can be found in [@wason2013EAGEobs]. 
We solve the factorization based rank-minimization formulation to recover the 
conventionally sampled seismic line (from the time-jittered data). We fixed the rank ``(k)`` of two low-rank factors to be 30. Note here that, in the proposed framework the rank value does not have a 
major effect on the reconstruction as long as the chosen rank is not too big or not too small [@aravkin2014SISCfmd].
 Figure [#mode] visualizes the different steps of the sampling-transformation operator ``\mathcal{A}``.
Figures [#results](e, f) show that the proposed rank-minimization based techniques 
effectively deblend the jittered data along with interpolation to the finer grid of 12.5m. This observation is also supported by the residual plots which show very small coherent energy (Figures [#results](g, h)). We further compare the proposed methodology to the sparsity-promoting
based techniques (Figures [#results](a, b)). As shown in Table [#SNRtable], although, the reconstruction quality is
 similar for both the methods, the computational speed up is by a factor of 8 and the storage of deblended data is reduced by a factor of 27 using rank-minimization based methods.

#### Table: {#SNRtable}

|   | SNR | | Time(hr.) | |Storage (gb)|
|--- |--- |--- |--- |--- |--- |--- |--- |
|            | Curvelet| low-rank|Curvelet| low-rank|Curvelet| low-rank|
| 50 m to 25 m  | 14.6 | 14.5 | 14 | 1.9 | 2 | 0.09 |
| 50 m to 12.5 m   | 11.3 | 12.6 | 16 | 2.1 | 2 | 0.09 |

:Qualitative (SNR) and Quantitative (computation time) analysis of different deblending algorithms.  

#### Figure: {#mode}
![](Fig/sop.png){#isame width=100%}
Visualization of sampling-transformation operator steps during the forward-operation of ``\mathcal{A}``. The output of sampling-transform operator is blended jittered marine data (showing only 30 seconds of the jittered data volume).

#### Figure: {#results}
![](Fig/SR_true){#isame width=40%}
![](Fig/MH_true){#isame width=40%}
![](Fig/SR_adj){#isame width=40%}
![](Fig/MH_adj){#isame width=40%}
![](Fig/SVD_true){#isame width=40%}
![](Fig/SVD_adj){#isame width=40%}
Monochromatic frequency slice at 10 Hz in the s-r domain and m-h domain for data sampled at 25m. (a,b) conventional data, and (c,d) blended data set. Decay of singular values for a frequency slice
(at 10 Hz) in the source-receiver (s-r) and midpoint-offset (m-h) domain of (e) conventional, and (f) blended seismic data set. Notice that, decay of singular values decay of fully sampled data matrix (conventional acquisition) have low-rank structure in the m-h domain as compared to s-r domain, whereas, blending increases the rank of the resulting frequency slice in the m-h domain.

#### Figure: {#results}
![](Fig/curv_rec){#isame width=40%}
![](Fig/curv_shot){#isame width=40%}
![](Fig/curv_diff_rec){#isame width=40%}
![](Fig/curv_diff_shot){#isame width=40%}
![](Fig/rec_rec){#isame width=40%}
![](Fig/rec_shot){#isame width=40%}
![](Fig/diff_rec){#isame width=40%}
![](Fig/diff_shot){#isame width=40%}
(a,b) Sparsity-promotion based recovery (SNR = 11.3 dB). (e,f) Rank-minimization based recovery (SNR = 12.6 dB). (c,d, g,h) Corresponding residual for the data sampled at 12.5m.

## Conclusions 

We propose a factorization based rank-minimization formulation for simultaneous source separation (deblending) and interpolation of seismic data in the case of time-jittered simultaneous marine acquisition. We showed that the proposed methodology is computationally feasible, both in terms of speed and memory usage, for the large-scale seismic data volumes. The separation quality of the proposed rank-minimization based technique is comparable to that of the sparsity-promoting based techniques. The experimental results demonstrate the potential benefit of this methodology. In the current formulation, we assume that the source vessel fires on the grid (discrete jittering) as opposed to firing off the grid (more realistic scenario). Future works involve the incorporation of irregularity along the acquisition grid and extension of this framework to 3-D seismic data acquisition. 

## Acknowledgements

This work was financially supported in part by the Natural Sciences
and Engineering Research Council of Canada Discovery Grant (RGPIN
261641-06) and the Collaborative Research and Development Grant DNOISE
II (CDRP J 375142-08). This research was carried out as part of the
SINBAD II project with support from the following organizations: BG
Group, BGP, Chevron, ConocoPhillips, CGG, ION GXT, Petrobras, PGS,
Statoil, Sub Salt Solutions, Total SA, WesternGeco, Woodside.


```math_def
\def\argmin{\mathop{\rm arg\,min}}
\def\vec{\mbox{```mathrm{vec}``}}
\def\ivec{\mbox{```mathrm{vec}^{-1}``}}

%\newcommand{\Id}{\mbox{``\tensor{I}\,``}}

\def\pector#1{\mathrm{\mathbf{#1}}} 
\def\cector#1{#1} 
\def\censor#1{#1} 
\def\vector#1{\mathbf{#1}}
\def\fvector#1{{\widehat{\vector{#1}}}}
\def\evector#1{{\widetilde{\vector{#1}}}}
\def\pvector#1{{\breve{\vector{#1}}}}
\def\pector#1{\mathrm{#1}}
\def\ctensor#1{\bm{\mathcal{#1}}}
%\def\tensorm#1{{#1}}
\def\tensorm#1{\bm{#1}}
%\def\tensor#1{\bm{#1}}7
\def\tensor#1{\vector{#1}}
%\def\cector#1{{\vector{\underline{#1}}}}
\def\hensor#1{\tensor{#1}}
\def\uensor#1{\underline{\bm{#1}}}
%\def\hensor#1{{\boldsymbol{\mathsf{#1}}}}
\def\hector#1{\vector{#1}}
%\def\hector#1{{\boldsymbol{\mathsf{#1}}}}
% \def\hfector#1{\hat{\boldsymbol{\mathsf{#1}}}}
% \def\ctensor#1{{\tensor{\underline{\underline{\mathsf{#1}}}}}}
\DeclareMathOperator\minim{\hbox{minimize}}
\def\ftensor#1{{\widehat{\tensor{#1}}}}
\def\calsor#1{{\boldsymbol{\mathcal{#1}}}}
\def\optensor#1{{\boldsymbol{\mathcal{#1}}}}
\def\hvector#1{\hat{\boldsymbol{\mathbf{#1}}}}
\def\uvector#1{{\underline{\vector{#1}}}}
\def\utensor#1{{\underline{\tensor{#1}}}}

%

```