Time-jittered marine acquisition — low-rank v/s sparsity

Introduction

Simultaneous (or blended) marine acquisition mitigates the challenges posed by conventional marine acquisition in terms of sampling and make large area acquisition inexpensive when compared to conventional acquisition (Beasley et al., 1998; Kok and Gillespie, 2002; Berkhout, 2008; Beasley, 2008; Hampson et al., 2008). 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. (Mansour et al., 2012, Wason and Herrmann (2013)) 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 (Maraschini et al., 2012; Cheng and Sacchi, 2013; Wason et al., 2014). 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 (Aravkin et al., 2014), 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 (Lee et al., 2010) with the Pareto curve approach proposed by (Berg and Friedlander, 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.
(Recht et al., 2010) 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: \[ \begin{equation} \label{BPDN} \min_{X} ||X||_* \quad \text{s.t. } \|\mathcal{A} (X) - b\|_2 \leq \epsilon, \tag{BPDN$_\epsilon$} \end{equation} \] 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, \(\ref{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(Mansour et al., 2012) resulting in a single long 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 \[ \begin{equation*} \min_{X} \sum_{i}^{n_f} ||X_i||_* \quad \text{s.t. } \|\mathcal{A} (X) - b\|_2 \leq \epsilon, \end{equation*} \] 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 \(\ref{BPDN}\), we use an extension of the SPG\(\ell_1\) solver (Berg and Friedlander, 2008) developed for the \(\ref{BPDN}\) problem in (A. Aravkin et al., 2012). The SPG\(\ell_1\) algorithm finds the solution to the \(\ref{BPDN}\) problem by solving a sequence of LASSO subproblems: \[ \begin{equation} \label{LASSO} \min_{X} \|\mathcal{A} (X) - b\|_2 \quad \text{s.t. } \sum_{i}^{n_f}||X_i||_* \leq \tau, \tag{LASSO$_\tau$} \end{equation} \] 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 (Rennie and Srebro, 2005; Lee et al., 2010; Recht and Ré, 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 (Rennie and Srebro, 2005), the nuclear norm obeys the relationship: \[ \begin{equation} \label{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), \end{equation} \] where \(\|\cdot\|_F^2\) is the Frobenius norm of the matrix (sum of the squared entries). Consequently, the LASSO subproblem can be replaced by \[ \begin{equation} \label{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\;, \end{equation} \] 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 (Wason and Herrmann, 2013). 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 (Aravkin et al., 2014). Figure 1 visualizes the different steps of the sampling-transformation operator \(\mathcal{A}\). Figures 2(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 2(g, h)). We further compare the proposed methodology to the sparsity-promoting based techniques (Figures 2(a, b)). As shown in Table 1, 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.

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.