---
title: Matrix completion on unstructured grids --- 2-D seismic data regularization and interpolation
author: |
    R. Kumar``^1``, O. Lopez``^{2}`` , E. Esser``^{1}`` and F. J. Herrmann``^1``\
    ``^{1}`` Seismic Laboratory for Imaging and Modeling (SLIM)\
	``^2`` Department of Mathematics\
	University of British Columbia, Canada
bibliography:
	- icml.bib   
---


## Abstract:

Seismic data interpolation via rank-minimization techniques has been recently introduced in the seismic community. All the existing rank-minimization techniques assume the recording locations to be on a regular grid, e.g. sampled periodically, but seismic data are typically irregularly sampled along spatial axes.  Other than the irregularity of the sampled grid, we often have missing data. In this paper, we study the effect of grid irregularity to conduct matrix completion on a regular grid for unstructured data. We propose an improvement of existing rank-minimization techniques to do regularization. We also demonstrate that we can perform seismic data regularization and interpolation simultaneously. We illustrate the advantages of the modification using a real seismic line from the Gulf of Suez to obtain high quality results for regularization and interpolation, a key application in exploration geophysics.


## Introduction 
Due to budgetary and/or physical constraints, Wide-azimuth (WAZ) data sets are coarsely sampled at irregular locations i.e., the samples are unstructured because they do not fall on a periodic grid.  Seismic data regularization and/or interpolation are one of the keys pre-processing steps to provides the regularly sampled data, with a dense sampling rate (without aliasing) in all the spatial directions.
Various methodologies have been proposed to perform regularization and interpolation. The simplest regularization process is called bin-centering, which moves traces from their recording locations to locations on a regular grid. Based upon the method used to perform bin-centering, amplitude of the traces may be altered and some traces may be discarded that can lead to serious errors, for example, as reported by @hennenfent2010GEOPnct. More recently, rank-minimization based techniques have been introduced to interpolate the seismic data [@oropeza:V25; @kumar2013EAGEsind]. The key idea of rank-minimization is to exploit the low-rank structure of seismic data in some "transform-domain" when organized in a matrix. The low-rank structure corresponds to a small number of nonzero singular values or quickly decaying singular values. @kumar2013EAGEsind showed that the monochromatic frequency slices of the fully sampled data matrix on a regular grid have low-rank structure in the midpoint-offset (m-h) domain, while sub-sampled seismic-data matrices do not. Missing traces increase the rank or make the singular values decay less quickly in the m-h domain, an essential feature for rank-minimization techniques to be effective. 
The existing rank-minimization techniques assume the input data is on a regularly sampled grid. As a result, these methods are less efficient when applied to an unstructured grid because discarding (binning) the actual recording locations of the input traces introduces errors. 
The objective of this paper is to establish the benefits of incorporating the unstructured sampling operator inside rank-minimization framework via performing regularization and interpolation of seismic data that are regularly sampled along receivers, time, and irregularly sampled along sources. We demonstrate the efficacy of the proposed extension on a real seismic line from the Gulf of Suez.

##Regularization
During acquisition, seismic data is acquired in an irregular fashion i.e., sources and receivers may lie on an unstructured grid. The goal of data regularization is to preserve the low-rank structure of seismic data in transform-domain while creating a regular sampled data at the specified bins from the irregular input data.
To illustrate the effects of regularization on the low-rank structure of seismic data in the m-h domain, we plot the decay of singular values in the m-h domain. We extract a fully sampled monochromatic frequency slice from a regularly sampled seismic line at high frequency (35 Hz) and transformed it into the m-h domain as shown in Figures [#svd]a. The 355 sources and receivers position for this frequency slice are at the centre of 12.5m, contiguous bin. From this fully sampled frequency slice; we generate 355 new sources at uniformly random locations within each bin. The data at unstructured grid is either binned using nearest-neighbor interpolation (Figure [#svd]b) or regularized using unstructured sampling operator based rank-minimization (Figure [#svd]c). Notice that, the decay of singular values is faster for the original data at structured grid in m-h domain. Regularization using bin-centering slow down the decay of singular values in m-h domain since binning \emph{breaks} the continuity along the wavefields, while unstructured sampling operator based rank-minimization \emph{preserves} the continuity along the wavefields, therefore, \emph{does not} changes the decay of singular values in m-h domain (Figure [#svd]d). Hence, incorporation of the unstructured sampling operator in rank-minimization perpetuates the low-rank structure of seismic data in the transform-domain.

##Regularization and Interpolation
In reality, seismic data is under sampled either along sources or receivers on irregular grids. Thus, regularization and missing trace interpolation problem can be perceived as a matrix-completion problem.
Let ``X_0`` in ``\mathbb{C}^{n\times m}`` be a regularly sampled data matrix and let ``\mathcal{A}`` be a linear measurement operator that maps from ``\mathbb{C}^{n\times m} \rightarrow \mathbb{C}^p`` with ``p \ll {n\times m}``. \cite{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 linear measurement operator is defined as ``\mathcal{A} := {{R}}{{\mathcal{S}}}^H ``, where ``{R}`` is the sampling operator, ``{{\mathcal{S}}}`` is the transformation operator from the source-receiver domain to the midpoint-offset domain and ``^H`` denotes the Hermitian transpose. The [#BPDN] formulation require regularly sampled data along all spatial axes, which is challenging in practice since a set of measurements ``b`` is always under sampled on an irregular grid. For that reason, we replace ``{R}`` with an unstructured sampling operator ``{{{U}}}`` which output the irregularly sampled data in the physical-domain by applying the fast Fourier transform followed by inverse non-equispaced fast Fourier transform [@Potts2001; @kunis06nfft].
To solve the nuclear-norm minimization problem, we combined the Pareto curve approach for optimizing [#BPDN] formulations with the SVD-free matrix factorization methods, following @aravkin2013SISCLR. 

##Results
The data set we use is a real seismic line from the Gulf of Suez with ``N_n = 355`` sources, ``N_m = 355`` receivers on a periodic interval of 12.5m and ``N_t=1024``  time samples with a temporal sampling of 0.004s. The 355 sources positions of this reference seismic line are at the center of 12.5m, contiguous bin. Most of the energy of the seismic line is concentrated inside 10-60 Hz frequency band. From this regularly sampled source positions (Figures [#data]a), we generate 355 new sources at uniformly random locations within each bin (Figures [#data]b). We use an unstructured sampling operator based on bi-cubic splines [@Keys]. This operator models taking samples from a densely periodically sampled data at source locations that are not on a grid. The nominal spatial sampling remains 12.5m. In all the subsequent experiments, we use 150 iterations of SPG``\ell_1`` for all frequency slices.
To map the data to a structured grid and to avoid binning errors, we conduct the regularization and compare it to the bin-centering techniques. Figure [#reg]a shows the resulting 12.5 m binned dataset using nearest-neighbor interpolation. Figure [#reg]c shows the detail along late times. The wiggle traces are the binned data and the grayscale image in the background is the ground truth. We can see that binning institute large errors (Figure [#reg]b), yielding a signal-to-noise ratio (SNR) of 7.3 dB since wavefronts need to be smooth and that needs to be dealt with properly in the sampling and binning does not do that properly. Figures [#reg]d, [#reg]f show the regularization results using unstructured [#BPDN] with a SNR of 18.6 dB. We can see that difference plot (Figure [#reg]e) shows very low-residual between the regularized data (Figure [#reg]d) and the ground truth (Figure [#data]a).
As mentioned before, seismic data are often inadequately or irregularly sampled along spatial axes. To imitate the under sampled data on an irregular grid, the fully sampled source locations at unstructured grid (Figure [#data]b) is used to generate 177 new sources at uniformly random locations (Figure [#data]c) where the minimum distance between two consecutive sources is 12.5m and the maximum distance is 130m. To show the benefits of incorporating the unstructured sampling operator in the rank-minimization, we also conduct a combined interpolation and regularization using rank-minimization techniques [@kumar2013EAGEsind], where we accomplish binning and interpolation simultaneously. Figure [#interp]a shows the regularized and interpolated data on unstructured grids using [#BPDN] with a SNR of 19.3 dB. Figure [#interp]c shows the details along late times. The difference plot (Figure [#interp]b) shows very low-residual, which explains that, the incorporation of grid irregularity in rank-minimization helps to achieve better regularization and interpolation results.


#### Figure: {#svd}
![](Fig1/D120mh){#isame width=40%}
![](Fig1/D250irreg){#isame width=40%}
![](Fig1/D250irreg_spgl1){#isame width=40%}
![](Fig1/svd_fullysampled_bin_hss){#isame width=40%}
Impact of regularization on singular values decay in the transform-domain. Monochromatic frequency slice at 35 Hz in the m-h domain, (a) ground truth, (b) binning using nearest-neighbor interpolation, (c) unstructured sampling operator based rank-minimization. (d) Singular values decay in the m-h domain. Notice that, the decay of singular values is \emph{faster} for the original data at a regularly sampled grid in m-h domain, and that regularization using bin-centering \emph{slow down} the decay of singular values in m-h domain, while unstructured sampling operator based rank-minimization \emph{does not} changes the decay of singular values in m-h domain. Therefore, incorporation of the unstructured sampling operator in rank-minimization preserve the low-rank structure of seismic data in the transform-domain.

#### Figure: {#data}
![](Fig1/wav){#isame width=40%}
![](Fig1/wav_irreg){#isame width=40%}
![](Fig1/wav_irreg_samp){#isame width=40%}
Ground truth (common receiver gather).
  (a) structured grid.
  (b) unstructured grid.
  (c) ``50\%`` subsampled data at unstructured grid.

#### Figure: {#reg}
![](Fig1/bwav){#isame width=40%}
![](Fig1/diffbwav){#isame width=40%}
![](Fig1/zbwav){#isame width=40%}
![](Fig1/bshot_reg){#isame width=40%}
![](Fig1/difbshot_reg){#isame width=40%}
![](Fig1/zbshot_reg){#isame width=40%}
Regularization. (a) Recovered data, binned using nearest-neighbor interpolation at nominal spatial sampling of 12.5 m. (d) Recovery using proposed rank-minimization on unstructured grid at 12.5 m. (c) and (f) zooms on late times. The grayscale image in the background is the ground truth, the wiggle traces are the binned data of (a) and (d), respectively. (b) and (e) are difference plots at the same scale as (a) and (d).
  

#### Figure: {#interp}
![](Fig1/wav_interp){#isame width=40%}
![](Fig1/diffnfft){#isame width=40%}
![](Fig1/zwav_interp){#isame width=40%}
Regularization and Interpolation. (a,b) Recovery and difference using proposed rank-minimization on unstructured grids with a SNR of 19.3 dB. (c) zooms on late times. The grayscale image in the background is the ground truth, the wiggle traces are the regularized-interpolated data of (a) respectively. We can see that incorporation of structured grid irregularity in rank-minimization helps to achieve better regularization and interpolation results.
  
## Conclusions 

We demonstrate the benefits of incorporating the grid irregularity for an unstructured data in the rank-minimization to achieve data regularization and interpolation. We establish that irregularity along spatial axes benefits the unstructured sampling operator based rank-minimization to perform better regularization of fully sampled irregular data resulting in minimal loss of coherent energy, since unstructured sampling operator preserves the low-rank structure of data in the transform-domain. We also illustrate that the proposed alteration of rank-minimization can be used to perform regularization and interpolation, simultaneously, from irregular and/or aliased data.  

## 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}}}}

%

```