---
title: A simulation-free seismic survey design by maximizing the spectral gap
author: |
    Yijun Zhang^1^, Mathias Louboutin^1^ , Ali Siahkoohi^1^, Ziyi Yin^1^, Rajiv Kumar^2^, and Felix J. Herrmann^1^ \
    ^1^ Georgia Institute of Technology \
    ^2^ Schlumberger \
bibliography:
    - abstract.bib
---


## Abstract:

Due to the tremendous cost of seismic data acquisition, methods have been developed to reduce the amount of data acquired by designing optimal missing trace reconstruction algorithms. These technologies are designed to record as little data as possible in the field, while providing accurate wavefield reconstruction in the areas of the survey that are not recorded. This is achieved by designing randomized subsampling masks that allow for accurate wavefield reconstruction via matrix completion methods. Motivated by these recent results, we propose a simulation-free seismic survey design that aims at improving the quality of a given randomized subsampling using a simulated annealing algorithm that iteratively increases the spectral gap of the subsampling mask, a property recently linked to the quality of the reconstruction. We demonstrate that our proposed method improves the data reconstruction quality for a fixed subsampling rate on a realistic synthetic dataset.

\vspace*{-0.5cm}
## Introduction
\vspace*{-0.25cm}

Due to relatively recent breakthroughs in compressive sensing [@candes2006robust], seismic data is increasingly gathered randomly along spatial coordinates to decrease acquisition productivity [@mosher2014increasing; @kumar2015efficient; @chiu20193d]. Wavefield reconstruction is used to recover fully sampled data from randomly subsampled observed seismic data [@hennenfent2008simply; @kumar2015efficient; @zhang2020wavefield]. To remove the imprint of large gaps in uniform random sampling, @gilles2008sampling proposed jittered subsampling, which by controlling the maximum gap size of subsampled data creates favorable conditions for seismic wavefield recovery based on sparsity promotion in a transformed domain made of localized atoms including curvelets [@herrmann2008curvelet]. While uniform random [@candes2009exact; @candes2010power] and random jittered subsampling schemes [@herrmann2008non; @hennenfent2008simply] are relatively straightforward to generate, these sampling strategies are almost certainly suboptimal and have shown to be improvable by solving certain optimization problem [@mosher2014increasing; @manohar2018data; @li2017compressive]. For instance, @mosher2014increasing\, @li2017compressive\, and later @titova2019mutual improved the reconstruction quality by devising a global optimization scheme that uses the mutual coherence.

In addition to wavefield reconstruction with optimized sampling schemes, @mosher2014increasing also proposed a simulation-based acquisition design to support the use of compressive sensing in seismic data acquisition. For time-lapse seismic, @guo2020data also used a data-driven approach where the acquisition is optimized by using prior information on the seismic data [@manohar2018data]. While these methods have lead to promising results, they either require significant computational resources to determine the optimal source-receiver layout using combined wavefield simulations and recoveries or require detailed information on the to-be-collected seismic data. Neither is feasible for the design of optimized sampling strategies in 3D.

To overcome this difficulty, this abstract provides a global optimization strategy for determining improved source-receiver layouts suitable for wavefield reconstructions based on matrix completion [@recht2010guaranteed; @kumar2015efficient; @kumar2017enabling] without the need to carry out expensive wavefield simulations. Similar to @li2017compressive and @mosher2014increasing\, who propose a simulation-free optimization method based on the mathematical property of mutual coherence for transform-based wavefield reconstruction, our method involves improving the connectivity of graphs spanned by the binary sampling masks in the midpoint-offset domain. According to @bhojanapalli2014universal\, by maximizing the spectral gap---i.e., the distance between the two first singular values---of the binary sampling mask the connectivity of the graph is improved, which favors reconstruction by matrix completion, an observation recently confirmed by @lopez2022graph for ``2``D and ``3``D seismic wavefield reconstruction. While recent work by @lopez2022graph indeed negates the need to run multiple costly wavefield reconstructions for different candidate sampling masks, this work does not yet provide a constructive method to generate sampling masks that maximize the spectral gap.

Unfortunately, the design of acquisition masks that maximize the spectral gap is an NP-hard problem [@guo2020data] whose solution requires a brute-force search through all combinatorial possibilities [@manohar2018data]\. When the number of sources or receivers becomes "large" [@li2016optimum], this precludes its practical use; for example, there are ``{n \choose m} = \frac{n!}{m!(n-m)!} = 75287520`` possible combinations when selecting ``m = 5`` subsampling positions from a pool of ``n = 100`` candidate sites. For this reason, we propose to obtain an approximate solution by maximizing the spectral gap using simulated annealing [@kirkpatrick1983optimization], a stochastic local search optimization technique that is straightforward to implement, apply, and computationally feasible. The proposed method depends only on binary mask optimization, has a minimal computational cost, and should be adaptable to large-scale survey design. 

We organize this expanded abstract as follows. First, we present the proposed optimization problem to maximize the spectral gap of subsampling masks. Next, we explain how to obtain the approximate acquisition masks  via simulated annealing. We conclude by demonstrating numerical experiments on the ``2``D synthetic Compass dataset [@jones2012building] and show an improvement in recovery quality compared to wavefield reconstruction of data collected with jittered subsampling method[@hennenfent2008simply].

\vspace*{-0.5cm}
## Methodology
\vspace*{-0.25cm}

Successful matrix-completion based seismic wavefield reconstruction [@recht2010guaranteed, @kumar2015efficient, @kumar2017enabling] hinges on three critical factors, namely: (1) an appropriate randomized subsampling scheme, such as uniform random or jittered subsampling [@herrmann2008non, @hennenfent2008simply]\; (2) existence of a transform domain in which the fully sampled seismic data organized as a matrix exhibits low-rank structure; (3) a computationally scalable matrix completion technique, which exploits the property that missing source and/or receivers increases the rank of these matrices. In this study, we propose a new constructive method to automatically generate improved source/receiver sampling masks, which favor seismic wavefield reconstruction via matrix completion in the midpoint-offset domain. We begin by describing our approach to acquisition design.

### Spectral ratio based acquisition design

Following @lopez2022graph\, we constitute the spectral gap by the spectral ratio (SR, the ratio of the first to second singular values), which becomes small for a large spectral gap. While the SR indeed has been shown to be a valuable quantity to predict the quality of wavefield reconstruction with matrix completion [@lopez2022graph], ways to automatically generate acquisition masks with small SRs have so far been lacking. To meet this challenge, we cast the problem of finding optimized acquisition masks with small SRs as a minimization problem. Given ``n_s`` source locations, ``n_r`` receiver locations, and the source subsampling ratio ``r``, we propose to solve a non-convex combinatorial optimization problem with respect to the subsampling mask ``\mathbf M \in \{0,1\}^{n_s \times n_r}``---i.e., we have
```math #eqOptimal
\mathcal L(\vector{M}) = & \minim_{\vector{M}} \quad \frac {\sigma_2(\mathcal{S}( \mathbf{M} ))}{\sigma_1(\mathcal{S}(\mathbf{M}))} \\

& \text{subject to} \\
& \|\mathbf{M}\|_0 = \lfloor n_s \times r\rfloor \times n_r \cap \mathbf{M} \in \mathcal{J} \cap \mathbf{M} \in \{0,1\}^{n_s \times n_r}.
```

In this optimization problem, the objective function consists of the spectral ratio (SR), defined by the ratio of the first, ``\sigma_1 (\cdot)``, and second, ``\sigma_2 (\cdot)``, singular values. ``\mathcal S`` stands for the transformation operator with seismic reciprocity [@fenati1984seismic] from the source-receiver domain to the midpoint-offset domain. We constrain the solution to conserve the subsampling ratio (``\|\mathbf{M}\|_0 = \lfloor n_s \times r\rfloor \times n_r ``) and to stay jittered sampled with ``\mathcal{J} \subseteq \{0,1\}^{n_s \times n_r}`` being the set of all possible jittered subsampling acquisitions. This constraint guarantees that the spread of the survey will not be modified, but only the local source position will be optimized. The symbol ``\lfloor \cdot \rfloor`` denotes the floor operation. As we previously mentioned, in order to solve this combinatorial optimization problem, we implemented a simulated annealing method to obtain a solution in a finite and acceptable time. We now describe this algorithm and link each step to its subsampling mask counterpart.

### Simulated annealing
\vspace*{-0.25cm}

Stochastic local search optimization algorithms are viable approximate methods for solving combinatorial optimization problems (e.g., Equation #eqOptimal)\. Simulated annealing is a global optimization technique that uses local search to find approximate solutions to combinatorial optimization problems given a computational budget [@kirkpatrick1983optimization, @csahin2010simulated].

This optimization method has three main components[@van1987simulated]: (1) an initial state, $\mathbf{M}_0$, representing the initial solution to the optimization problem (Equation #eqOptimal\); (2) a set of neighboring states for any given state, which will be used to update the current state randomly; and (3) a transition probability that determines the probability of moving from one state to another. During optimization, at each given state, $\mathbf{M}_k$, which represents the current solution to the optimization problem #eqOptimal\, a candidate state, $\tilde{\mathbf{M}}_k$, is chosen randomly from the neighboring states. Next, the algorithm transitions from the current state to the candidate state, i.e., from $\mathbf{M}_k$ to $\tilde{\mathbf{M}}_k$, if this transition reduces the objective function, i.e., ``\mathcal{L}(\tilde{\mathbf{M}}_k) < \mathcal{L}(\mathbf{M}_k)``. On the other hand, if the objective function evaluated at the candidate state is larger than the current value, the algorithm makes the transition to the candidate state according to a transition probability, defined as follows [@kirkpatrick1983optimization]:

```math #eqProb
p(\delta{\mathcal L}, k) = \exp \big({\frac{-\delta {\mathcal L}}{T(k)}}\big),
```
where ``k`` is the iteration number, ``\delta\mathcal{L}=\mathcal{L}(\tilde{\mathbf{M}}_k) - \mathcal{L}(\mathbf{M}_k)`` indicates the change in the objective function (Equation #eqOptimal\) by moving to the candidate state, and ``T(k): \mathbb{R} \rightarrow \mathbb{R}^{+}``\, typically called temperature function\, is a monotonically decreasing function that reduces the uphill transition probability towards the end of optimization while allowing uphill movement early in the optimization. We choose the temperature function as ``T(k)= T_0 \times \alpha^k``, following a geometric reduction rule, which is the most commonly used function in the simulated annealing with a start temperature ``T_0`` and the decrease rate ``\alpha`` [@kirkpatrick1983optimization, @abramson1999simulated]. This allows the algorithm to escape from local minima in the initial stages of the optimization while ensuring downhill movement towards the end. Finally, the transition probability is smaller for candidate states that increase the objective function more, i.e., ``\delta\mathcal{L} \gg 0``, minimizing the probability of moving to very bad solutions.

To adapt simulated annealing to the acquisition design optimization problem (cf. Equation #eqOptimal\), we define the states as arbitrary positioning of sources. The algorithm is summarized in Algorithm #alg\. This algorithm is initialized with a subsampling mask ``\mathbf{M}_0`` that is generated by using jittered subsampling method, known to facilitate seismic wavefield recovery [@hennenfent2008simply]. After updating the temperature function ``T(k)`` (line ``1``) [@ma2002simultaneous], we select a source position ``\tilde{\mathbf{M}}_k`` within the neighborhood of the current position (line ``2``). We then update the source position to this updated state according to the loss decrease and probabilistic update rule (lines ``3-8``). After a predetermined number of iterations, the algorithm outputs the source sampling mask ``\mathbf{M}_K`` with smaller SR.

### Algorithm: {#alg}

| **Inputs:**
|		``\mathbf{M}_0``\: Initial source positions using jittered subsampling method.
|		``K``\: Maximum number of iterations.
|		``T(k)``\: Temperature function.
|		``p``\: Transition probability (cf. Equation #eqProb\).
| 0. **for** ``k = 0`` **to** ``K-1`` **do**
| 1.    ``T(k)= T_0 \times \alpha^k``
| 2.	``\tilde{\mathbf {M}}_k \leftarrow`` randomly pick a neighboring state
| 3.	``\delta \mathcal {L} = \mathcal{L}(\tilde{\mathbf{M}}_k) - \mathcal{L}(\mathbf{M}_k)``
| 4. 	**if** ``\delta \mathcal {L} \leq 0``
| 5. 		``\mathbf{M}_{k+1} = \tilde{\mathbf{M}}_k``
| 6. 	**else**
| 7.		``\mathbf{M}_{k+1} = \begin{cases}
                                 &\tilde{\mathbf{M}}_k \text { with probability } p(\delta{\mathcal L}, k) \\
								 &{\mathbf{M}}_k
								 \end{cases}``
| 8. 	**end if**
| 9. **end for**
| **Output:** ``\mathbf{M}_K``\
: Spectral ratio minimization via simulated annealing.

In order to satisfy the constraints introduced in Equation #eqOptimal, we carefully define the neighborhood of acceptable state to prevent sources positions to cluster around a few areas of the survey. We now detail its design.


#### Neighboring states
\vspace*{-0.25cm}

Given the source sampling at an iteration, we randomly select 20% of the source positions in the current state to balance between exploring the search space and avoiding too large change between adjacent iterations [@olorunda2008measuring, @assad2018hybrid]. Using the subsampling factor ``f = \frac{1}{r} ``, the fine grid with all possible source locations is divided into ``f`` equal regions. Each selected source is allowed to randomly shift within the region in which it is located. For clarity, we summarize the perturbation rule in Figure #example_2 with ``n_s = 20`` source positions and a subsampling factor of ``f=5``. We choose the movement range ``R`` to ensure that we remain close to the jittered sampling, which has been shown to result in better wavefield recoveries [@hennenfent2008simply]. We now detailed our synthetic numerical experiment demonstrating the benefits of our method for data reconstruction.

#### Figure: {#example_2}
![](results/select_neighbor1.png){width=100%}
: Random state perturbation rules to satisfy the constraints. White circles (``\circ``) indicate all possible source locations. Five red circles represent an initial state. The blue circle represents the ``20\%`` of sources that we will move within movement range ``R`` to arrive at a neighboring state (purple circle).

\vspace*{-0.5cm}
## Numerical experiments
\vspace*{-0.25cm}

We consider a 2D marine dataset simulated over the realistic Compass model [@jones2012building]. The dataset consists of ``300`` sources and ``150`` receivers sampled at ``12.5\, \mathrm{m}``. The data is recorded at a ``2 \, \mathrm{ms}`` sampling rate for ``2.046 \, \mathrm{s}`` (``1024`` time samples). Based on this dataset, we proceed in two steps. First, we will compare the subsampling mask we obtain with our proposed method against the standard jittered sampling mask. Second, we will show that the recovered data is of better quality as expected from the optimal SR that quantifies the expected quality of recovery. The jittered subsampling method used in this abstract allows neighboring subsamples (non-gap subsamples), which creates favorable conditions for recovery and is defined as optimally-jittered subsampling in @hennenfent2008simply\. We use the weighted matrix completion method [@zhang2020wavefield] to recover the observed data and evaluate the quality of the recovered dataset.

We start with picking a subsampling mask using the jittered method [@herrmann2008non] that includes 20% of the sources. In Figure #OJ_SR, we show the subsampling mask in the source-receiver domain. Under the assumption of the source-receiver reciprocity [@fenati1984seismic], we apply this reciprocity on the subsampling mask to implement a realistic seismic survey design. The subsampling mask in the midpoint-offset domain with seismic reciprocity is depicted in Figure #OJ_MH\. With this mask as an initial guess, we perform ``4000`` iterations of simulated annealing to obtain an optimal subsampling mask. Figure #SA_SR and Figure #SA_MH show the resulting subsampling mask in the source-receiver and midpoint-offset domains, respectively. We observe that the SR was reduced by 30% with a fixed subsampling rate hinting towards well-improved data reconstruction. The proposed method is a simulation-free method that depends only on binary mask optimization.

#### Figure: {#example_1}
![](results/Initial_OJ_SR2.png){#OJ_SR width=40%}
![](results/Initial_OJ_MH2.png){#OJ_MH width=60%}\
![](results/Final_OJ_SR2.png){#SA_SR width=40%}
![](results/Final_OJ_MH2.png){#SA_MH width=60%}
 : Jittered subsampling mask in the *(a)* source-receiver domain and *(b)* midpoint-offset domain (SR ``= 0.331``). Optimized subsampling mask in the *(c)* source-receiver domain and *(d)* midpoint-offset domain (SR ``= 0.242``). 

With this optimized subsampling mask, we now perform data reconstruction via weighted matrix completion [@zhang2020wavefield] and compare the result against reconstructing the data sampled with the initial jittered subsampled mask. In both cases, we use the same algorithm and hyperparameters (e.g., number of iterations, rank) for a fair comparison. We summarize the recovery in Figure #example_4\. 

#### Figure: {#example_4}
![](results/Ground_Truth.png){#GT width=50%}
![](results/Observed_data.png){#Obs width=50%}\
![](results/OJ_Rel.png){#OJR width=50%}
![](results/SA_Rel.png){#SAR width=50%}\
![](results/OJ_Diff.png){#OJD width=50%}
![](results/SA_Diff.png){#SAD width=50%}
 : Wavefield reconstruction results in the time domain. *(a)* Ground truth. *(b)* ``80 \%`` subsampled seismic data with proposed subsampling. Reconstructions from ``80\%`` missing sources: *(c)* jittered subsampling, SNR = ``14.6\, \mathrm{dB}`` and ``12.8\, \mathrm{dB}`` for later arrival events, *(d)* improved subsampling with SNR = ``14.91\, \mathrm{dB}`` and ``12.9\, \mathrm{dB}`` for later arrival events. 

We first show the ground truth in Figure #GT\, where the right plot shows the full shot record and the left one depicts the later arrival events between about ``1\, \mathrm{s}`` to ``2\, \mathrm{s}``. By applying these two masks (jittered mask and proposed mask) individually on the ground truth, we obtain two observed datasets with ``80\%`` of sources missing. The proposed subsampled data is illustrated in Figure #Obs\. Figure #OJR shows the reconstruction from jittered observed data with a signal-to-noise ratio (SNR) of  ``14.6\, \mathrm{dB}`` for the full shot record and ``12.8\, \mathrm{dB}`` for the later arrival events. The reconstruction from the proposed subsampled data is shown in Figure #SAR\, with SNRs of  ``14.91\, \mathrm{dB}`` for the full shot record and ``12.9\, \mathrm{dB}`` for the later arrival events. Figure #OJD illustrates the difference between Figure #OJR and Figure #GT\, whereas Figure #SAD shows the difference between Figure #SAR and Figure #GT\. The wavefield reconstructions demonstrate that the reconstruction from the proposed subsampling mask gives a more accurate data reconstruction for the full shot record and the later arrival events in terms of SNR. The difference is significantly reduced in Figure #SAD in contrast to Figure #OJD\.

To further validate the performance of our proposed method, we show that our proposed mask outperforms on the average standard jittered acquisition and not just for a single experiment. We randomly generate five independent jittered subsampling masks (removing ``80 \%`` of sources) and then utilize the proposed approach to minimize the SR of these five jittered subsampling layouts. These five jittered and proposed masks are then used to perform weighted matrix completion [@zhang2019high, @zhang2020wavefield] and reconstruct the full dataset. The results are summarized in Figure #example_5\. The bar plots in Figure #example_5 lead to the following observations. First, our proposed method consistently reduced the spectral ratio by at least ``11\%`` leading to a similar optimal SR for this given subsampling ratio and acquisition. Second, the recovered data always presents a higher SNR representative of a more accurate wavefield reconstruction. These two results show that despite being a potentially aleatory method, our simulated annealing based SR minimization method consistently provides a subsampling mask best fitted for data reconstruction.

#### Figure: {#example_5}
![](results/SG_boxplot.png){#SG_boxplot width=50%}
![](results/SNR_boxplot.png){#SNR_boxplot width=50%}
: *(a)* SR comparison (lower is better) of subsampling masks using jittered method versus proposed method. *(b)* Reconstruction SNR comparison (higher is better) from observed data using jittered method and proposed method. The results are obtained by five independent experiments. 

\vspace*{-0.5cm}
## Conclusions
\vspace*{-0.25cm}

We proposed a simulation-free method for seismic survey design in this study by minimizing the spectral ratio using the simulated annealing algorithm. Because the proposed method solely relies on a binary mask optimization rather than being data-driven, the computational cost is minimal and should scale to industry-size survey design. Through analysis and experiments, we conclude that the proposed method leads to an optimal subsampling mask best fitted for wavefield reconstruction based on matrix completion. Future work will focus on applying the proposed method to the three-dimensional field data.


\vspace*{-0.5cm}
## Acknowledgement
\vspace*{-0.25cm}

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




```math_def
\def\argmin{\mathop{\rm arg\,min}}
\def\vec{\mbox{$\mathrm{vec}$}}
\def\ivec{\mbox{$\mathrm{vec}^{-1}$}}
%\newcommand{\rr}{{{r}}}
\newcommand{\m}{{\mathsf{m}}}
\newcommand{\PsDO}{\mbox{PsDO\,}}
\newcommand{\Id}{\mbox{$\tensor{I}\,$}}
\newcommand{\R}{\mbox{$\mathbb{R}$}}
\renewcommand{\C}{\mbox{$\mathbb{C}$}}
\newcommand{\E}{\mbox{$\mathbb{E}$}}
\newcommand{\Z}{\mbox{$\mathbb{Z}$}}
\newcommand{\DE}{:=}
\newcommand{\Order}{\mbox{${\mathcal O}$}} \def\bindex#1{{\mathcal{#1}}}
%\newcommand{\Id}{\mbox{$\tensor{I}\vector{d}\,$}}
\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\mathcalsor#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}}}}
%\
\newcommand{\dd}{\ensuremath{\mathrm{d}}}
\newcommand{\blockstack}{\operatorname{blockstack}}
\newcommand{\diag}{\operatorname{diag}}
\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
\newcommand{\rr}[1]{\mbox{\rm #1}}
\usepackage{algorithm2e}
\def\minim{\mathop{\hbox{minimize}}}
```