Compressive sensing in seismic exploration: an outlook on a new paradigm

Abstract

Many seismic exploration techniques rely on the collection of massive data volumes that are subsequently mined for information during processing. While this approach has been extremely successful in the past, current efforts toward higher resolution images in increasingly complicated regions of the Earth continue to reveal fundamental shortcomings in our workflows. Chiefly amongst these is the so-called “curse of dimensionality” exemplified by Nyquist’s sampling criterion, which disproportionately strains current acquisition and processing systems as the size and desired resolution of our survey areas continues to increase.

We offer an alternative sampling method leveraging recent insights from compressive sensing towards seismic acquisition and processing for data that, from a traditional point of view, are considered to be undersampled. The main outcome of this approach is a new technology where acquisition and processing related costs are decoupled from the stringent Nyquist sampling criterion.

At the heart of our approach lies randomized incoherent sampling that breaks subsampling-related interferences by turning them into harmless noise, which we subsequently remove by promoting sparsity in a transform-domain. Acquisition schemes designed to fit into this regime no longer grow significantly in cost with increasing resolution and dimensionality of the survey area, but instead its cost ideally only depends on transform-domain sparsity of the expected data. Our contribution is twofold. First, we demonstrate by means of carefully designed numerical experiments that ideas from compressive sensing can be adapted to seismic acquisition. Second, we leverage the property that seismic data volumes are well approximated by a small percentage of curvelet coefficients. Thus curvelet-domain sparsity allows us to recover conventionally-sampled seismic data volumes from compressively-sampled data volumes whose size exceeds this percentage by only a small factor. Because compressive sensing combines transformation and encoding by a single linear encoding step, this technology is directly applicable to seismic acquisition and therefore constitutes a new paradigm where acquisitions costs scale with transform-domain sparsity instead of the gridsize. We illustrate this principle by showcasing recovery of a real seismic line from simulated compressively sampled acquisitions.

1 Outline

In Part I, we showed by means of carefully designed experiments that the theory of compressive sensing can be successfully adapted to seismic-data acquisition. The experiments confirmed that accurate recoveries are possible from compressively sampled data volumes that exceed the size of conventionally compressed data volumes by only a small factor. In this part II, we adapt these findings to the recovery of a real seismic line from compressively-sampled ‘Land’ and ‘Marine’ data. Because compressive sensing combines transformation and encoding by a single linear encoding step, this technology is directly applicable to seismic acquisition and therefore constitutes a new paradigm where acquisitions costs scale with transform-domain sparsity instead of the grid size. We illustrate this principle by showcasing recovery of a real seismic line from simulated compressively sampled acquisitions.

2 An academic case study

Now that we established that high SNR’s are achievable with modest oversampling ratios, we study the performance of our recovery algorithm on a seismic line from the Gulf of Suez by comparing two simultaneous-source scenarios with coincident source-receiver positions:

Land acquisition with random amplitude encoding: Here, sequential impulsive sources are replaced by impulsive simultaneous `phase-encoded’ sources. Mathematically, simultaneous measurements are obtained by replacing the sampling matrix for the sources—normally given by identity matrix—by a measurement matrix obtained by phase encoding along the source coordinate. Following (Romberg, 2009) and (Herrmann et al., 2009b), we define the measurement matrix by the following Kronecker product

In this expression, conventional sampling, which corresponds to the action of the identity matrix I, is replaced by a ‘random phase encoding’ consisting of applying a Fourier transform along the source coordinate (F_s), followed by uniformly drawn random phase rotations ϑ ∈ [0, π], an inverse Fourier transform (F_s*), and multiplication by a random sign vector (i.e., multiplication by with diag (η) with η ∈ ± 1 with equal probability). As shown by (Romberg, 2009), the combined action of these operations corresponds to the action of a Gaussian matrix at reduced computational costs (see also Herrmann et al., 2009b). Application of this matrix to a conventionally-sampled seismic line turns sequential impulsive source into a simultaneous “supershot” where all sources fire simultaneously with weights drawn from a single Gaussian distribution. As before, the restriction operator selects a subset of n’_s “supershots” generated by different randomly-weighted simultaneous sources. After restriction along the source coordinate, the sampling matrix has an aspect (or undersampling) ratio of δ = n’_s / n_s . An example of this type of sampling, resulting in a seismic line consisting of n’_s / n_s supershots, is included in Figure 1. In this Figure, n_s single impulsive-source experiments (1(a)) become n’_s simultaneoussource experiments (juxtapose Figure 1(a) and 1(c)). While this sort of sampling is perhaps physically unrealizable—i.e., we typically do not have large numbers of vibroseis trucks available— it gives us the most favorable recovery conditions from the compressive-sensing perspective. Therefore, our `Land’ acquisition will serve as a benchmark with which we can compare alternative and physically more realistic acquisition scenarios.

Marine acquisition with random-time dithering: Here, sequential acquisition with a single airgun is replaced by continuous acquisition with multiple airguns that fire continuously at random times and at random locations. In this scenario, a seismic line is mapped into a single long `supershot’. Mathematically, this type of acquisition is represented by the following sampling operator

In this expression, the linear operator T turns sequential recordings (1(b)) with synchronized impulsive shots (Figure 1(a)) into continuous recordings with n*_s impulsive sources firing at random positions (Figure 1(e)), selected uniformly-random from [1···n_s] discrete source indices and from discrete random time indices, selected uniformly from [0···(n*_s – 1) x n_t)] time indices. Note that T acts both on the shot and the time coordinate. The resulting data is one long ‘supershot’ that contains a superposition of n*_s impulsive shots. For plotting reasons, we reshaped in Figure 1(f) this long record into multiple shorter records. Notice that this type of Marine acquisition is physically realizable as long as the number of simultaneous sources involved is limited.

Aside from mathematical factors, such as the mutual coherence (discussed in part 1) that determines the recovery quality, there are also economical factors to consider. For this purpose, (Berkhout, 2008) proposed two performance indicators, which quantify the cost savings associated with simultaneous and continuous acquisition. The first measure compares the number of sources involved in conventional and simultaneous acquisition and is expressed in terms of the source-density ratio

For Land data acquisition, this quantity equals SDR_Land = (n_s x n’_s) / n_s = n’_s and for Marine data SDR_Marine = (n*_s x n_s). Remember that the number of sources refers the number of sources firing and not the number of source experiments. Clearly, Land acquisition has a significant higher SDR.

Aside from the number of sources, the cost of acquisition is also determined by survey-time ratio

Ignoring overhead in sequential shooting, this quantity equals STR_Land = n’_s / n_s in the first and STR_Marine = n_s x T₀ / T with T₀ the time of a single sequential experiment and T the total survey time of the continuous recording. The overall economic performance is measured by the product of these two ratios. For Land acquisition this product is proportional to n_s and for `Marine’ acquisition proportional to n*_s x T₀ / T.

As we have seen from our discussion on compressive sensing, recovery depends on the mutual coherence of the sampling matrix. So, the challenge really lies in the design of acquisition scenarios that obtain the lowest mutual coherence while maximizing the above two economic performance indicators. To get a better insight in how these factors determine the quality of recovered data, we conduct a series of experiments by simulating possible alternative acquisition strategies on a perviously traditionally recorded real seismic line.

First, we simulate `Land’ data for δ - 0.5 ( 64 simultaneous source experiments with all sources firing) and study the recovery based on 2-D and 3-D curvelets. The former is based on a 2-D discrete curvelet transform along the source and receiver coordinates, and the discrete wavelet transform along the remaining time coordinate:

We conduct a similar experiment for the `Marine case’. In this case, we randomly select 128 shots from the total survey time T = δ x (n_s – 1) x T₀, yielding the same aspect ratio for the sampling matrix.

Figures 2 and 3 summarize the results for Land and Marine acquisition using recoveries based on the 2-D and 3-D curvelet transform. The following observations can be made. First, it is clear that accurate recovery is possible by solving an l₁ optimization problem using SPGl₁ (Berg and Friedlander, 2008) while limiting the number of iterations for the 2-D case to 500 and the 3-D case to 200. Second, the recovery results for 3-D recovery of Land data show an improvement of 1.3 dB by exploiting 3-D structure of the wavefronts. Similarly, we find an improvement of 3.9 dB for the Marine case. Both observations can be explained by the fact that the 3-D curvelet transforms attains higher sparsity because it explores continuity of the wavefield along all three coordinate axes. Second, Land acquisition clearly favors recovery by curvelet-domain sparsity promotion compared to Marine acquisition. This is true despite the fact that the subsampling ratio, i.e., the aspect ratio of the sampling matrices, are the same. Clearly this difference lies in the mutual coherence of the sampling matrix. The columns of the sampling matrix for Land acquisition are more incoherent and hence more independent and this favors recovery. These observations are confirmed by the SNRs, which for Land acquisition equal 10.3 dB and 11.6 dB, for the 2-D/3-D recovery, respectively, and 7.2 dB and 11.1 dB, for Marine acquisition.

Unfortunately, recovery quality is not the only consideration. The economics expresse so play a role. In the above setting, the Land acquisition has a SDR = 64 and STR = 2 while the Marine acquisition has SDR = 1 and STR = 2. Clearly, the SDR for land acquisition may not be realistic.

3 Discussion

The presented results illustrate that we are at the cusp of exciting new developments where acquisition workflows are no longer impeded by subsampling related artifacts. Instead, we arrive at acquisition schemes that control these artifacts. We accomplish by applying the following design principles: (i) randomize— break coherent aliases by introducing randomness, e.g. by designing randomly perturbed acquisition grids, or by designing randomized simultaneous sources; and (ii) sparsify—utilize sparsifying transforms in conjunction with sparsity-promoting programs that separate signal and subsampling artifacts and that restore amplitudes. The implications of randomized incoherent sampling go far beyond the examples presented here. For instance, our approach is applicable to land acquisition for physically realizable sources (Krohn and Neelamani, 2008; Romberg, 2009) and can be used to compute solutions to wavefield simulations (Herrmann et al., 2009b) and to compute full waveform inversion (Herrmann et al., 2009b) faster. Because randomized sampling is linear (Bobin et al., 2008), wavefield reconstructions and processing can be carried out incrementally as more compressive data becomes available.

Indeed, compressive sensing offers enticing perspectives towards the design of future Land and Marine acquisition systems. In order for this technology to become successful the following issues need to be addressed, namely the performance of recovery

• from field data including all its idiosyncrasies. This will require a concerted effort from practitioners in the field and theoreticians. For Marine acquisition, recent work by (Moldoveanu, 2010) has shown early indications that randomized jittered sampling leads to improved imaging.
• from discrete data with quantization errors. Addressing this issue calls for integration of digital-to-analog conversion into compressive and recent progress has been made in this area (see e.g. Gunturk et al., 2010);
• from Land data that has the imprint of statics. Addressing this issue will be essential because severe static effects may adversely affect transform-domain sparsity on which recovery from compressive-sampled data relies.

At the UBC Seismic Laboratory for Imaging and Modelling (http://slim.eos.ubc.ca/SLIM), we hope to report progress on these important topics in future publications.

4 Conclusions

We successfully demonstrated that compressive sensing can be used to recover seismic lines from simultaneous and continuously acquired data. This observation is a direct consequence of the fact that compressive sensing combines sampling and compression in a single linear encoding step. Because simultaneous and continuous acquisition can be identified as an instance of compressive sampling, this has profound implications for exploration seismology that include: a new randomized sampling paradigm, where the cost of acquisition are no longer dominated by resolution and size of the acquisition area, but by the desired reconstruction error and transform domain sparsity of the wavefield, and a new paradigm for randomized processing and inversion, where dimensionality reductions will allow us to mine high-dimensional data volumes for information in ways, which previously, would have been computationally infeasible.

Acknowledgments This paper was written as a follow up to the first author’s presentation during the ``Recent advances and Road Ahead’’ session at the 79th annual meeting of the Society of Exploration Geophysicists. I would like to thank the organizers of this session Dave Wilkinson and Masoud Nikravesh for their invitation. I also would like to thank Dries Gisolf and Eric Verschuur of the Delft University of Technology for their hospitality during my sabbatical and Gilles Hennenfent for preparing some of the figures. This publication was prepared using CurveLab, a toolbox implementing the Fast Discrete Curvelet Transform, WaveAtom a toolbox implementing the Wave Atom transform, Madagascar, a package for reproducible computational experiments, SPGl₁, SPOT, a suite of linear operators and problems for testing algorithms for sparse signal reconstruction, and pSPOT, SLIM’s parallel extension of SPOT. FJH was in part financially supported by NSERC Discovery Grant (22R81254) and by CRD Grant. The industrial sponsors of the Seismic Laboratory for Imaging and Modelling, BG Group, BP, Chevron, ConocoPhillips, Petrobras, Total SA, WesternGeco are also gratefully acknowledged.