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.
 Part I of this paper is adapted from Randomized sampling and sparsity: getting more information from fewer samples. Geophysics 75, WB173 (2010); doi:10.1190/1.350614.
 Seismic Laboratory for Imaging and Modeling, Department of Earth and Ocean Sciences, the University of British Columbia, 6339 Stores Road, Vancouver, V6T 1Z4, BC, Canada. Email: email@example.com
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 (Fs), followed by uniformly drawn random phase rotations ϑ ∈ [0, π], an inverse Fourier transform (Fs*), 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 / ns . An example of this type of sampling, resulting in a seismic line consisting of n’s / ns supershots, is included in Figure 1. In this Figure, ns 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···ns] discrete source indices and from discrete random time indices, selected uniformly from [0···(n*s – 1) x nt)] 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 SDRLand = (ns x n’s) / ns = n’s and for Marine data SDRMarine = (n*s x ns). 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 STRLand = n’s / ns in the first and STRMarine = ns x T0 / T with T0 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 ns and for `Marine’ acquisition proportional to n*s x T0 / 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 (ns – 1) x T0, 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 l1 optimization problem using SPGl1 (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.
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.
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, SPGl1, 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.
About the Author(s)
Felix J. Herrmann received in 1997 a Ph.D. degree in Engineering Physics from the Delft University of Technology, the Netherlands. He was a visiting scholar at Stanford in 1998, a post-doctoral fellow at MIT's Earth Resources Laboratory from 1999 to 2002, and was a senior Fellow at the UCLA's Institute for Pure and Applied Mathematics. He is currently associate professor at the Department of Earth & Ocean Sciences at the University of British Columbia where he founded the Seismic Laboratory for Imaging and Modeling (SLIM). He is interested in theoretical and applied aspects of (exploration) reflection seismology and runs the industry/government-supported research programs SINBAD & DNOISE. He serves on the editorial boards of the Journal of Applied Geophysics and Geophysical Prospecting. He is a member of the Institute of Applied mathematics (IAM) and an associate member of the Institute for Computing, Information, and Cognitive Science (ICICS). He is also a member of the EAGE, CSEG, SEG, SIAM, and AGU.
Haneet Wason is a graduate from the University of Calgary with a B.Sc. degree in Geophysics, and is currently pursuing a M.Sc. degree in Geophysics at the University of British Columbia.
Tim T.Y. Lin is a Ph.D. Geophysics student at the University of British Columbia's Department of Earth and Ocean Sciences. He is part of the Seismic Laboratory for Imaging and Modelling with a special interest in the confluence of physics and computational optimization algorithms, and its application to reflection seismic signal processing. He holds a combined B.Sc. degree in Physics and Computer Science from the University of British Columbia.
Abma, R., and N. Kabir, 2006, 3D interpolation of irregular data with a POCS algorithm: Geophysics, 71, E91–E97.
Berg, E. v., and M. P. Friedlander, 2008, Probing the Pareto frontier for basis pursuit solutions: Technical Report 2, Department of Computer Science, University of British Columbia.
Berkhout, A. J., 2008, Changing the mindset in seismic data acquisition: The Leading Edge, 27, 924–938.
Bleistein, N., J. Cohen, and J. Stockwell, 2001, Mathematics of Multidimensional Seismic Imaging, Migration and Inversion: Springer.
Bobin, J., J. Starck, and R. Ottensamer, 2008, Compressed sensing in astronomy: IEEE J. Selected Topics in Signal Process, 2, 718–726.
Candès, E., J. Romberg, and T. Tao, 2006, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223.
Candès, E. J., and L. Demanet, 2005, The curvelet representation of wave propagators is optimally sparse: Comm. Pure Appl. Math, 58, 1472–1528.
Candès, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2006a, Fast discrete curvelet transforms: Multiscale Modeling and Simulation, 5, 861–899.
Candès, E. J., J. Romberg, and T. Tao, 2006b, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information: IEEE Transactions on Information Theory, 52, 489 – 509.
Canning, A., and G. H. F. Gardner, 1996, g61, 1103–1114.
Demanet, L., and L. Ying, 2007, Wave atoms and sparsity of oscillatory patterns: Applied and Computational Harmonic Analysis, 23, 368–387.
Donoho, D., A. Maleki, and A. Montanari, 2009, Message passing algorithms for compressed sensing: Proceedings of the National Academy of Sciences.
Donoho, D., and J. Tanner, 2009, Observed universality of phase transitions in highdimensional geometry, with implications for modern data analysis and signal processing: Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367, 4273.
Donoho, D. L., 2006a, Compressed sensing: IEEE Trans. Inform. Theory, 52, 1289–1306.
——–, 2006b, Compressed sensing: IEEE Transactions on Information Theory, 52, 1289 – 1306.
Donoho, D. L., and X. Huo, 2001, Uncertainty principles and ideal atomic decomposition: IEEE Transactions on Information Theory, 47, 2845–2862.
Donoho, D. L., and B. F. Logan, 1992, Signal Recovery and the Large Sieve: SIAM Journal on Applied Mathematics, 52, 577–591.
Donoho, P., R. Ergas, and R. Polzer, 1999, Development of seismic data compression methods for reliable, low-noise performance: SEG International Exposition and 69th Annual Meeting, 1903–1906.
Fomel, S., 2003, Theory of differenctial offset continuation: 68, 718–732.
Fomel, S., J. G. Berryman, R. G. Clapp, and M. Prucha, 2002, Iterative resolution estimation in leastsquares Kirchoff migration: Geophys. Pros., 50, 577–588.
Güntürk, S., A. Powell, R. Saab, and Ö. Yılmaz, 2010, Sobolev duals for random frames and Sigma- Delta quantization of compressed sensing measurements: Arxiv preprint arXiv:1002.0182.
Hale, D., 1995, DMO processing: Geophysics Reprint Series.
Harlan, W. S., J. F. Claerbout, and F. Rocca, 1984, Signal/noise separation and velocity estimation: Geophysics, 49, 1869–1880.
Hennenfent, G., and F. J. Herrmann, 2008, Simply denoise: wavefield reconstruction via jittered undersampling: Geophysics, 73.
Herrmann, F. J., 2009, Compressive imaging by wavefield inversion with group sparsity: SEG Technical Program Expanded Abstracts, SEG, SEG, 2337–2341.
Herrmann, F. J., U. Boeniger, and D. J. Verschuur, 2007, Non-linear primary-multiple separation with directional curvelet frames: Geophysical Journal International, 170, 781–799.
Herrmann, F. J., Y. A. Erlangga, and T. Lin, 2009a, Compressive sensing applied to full-waveform inversion.
——–, 2009b, Compressive simultaneous full-waveform simulation: Geophysics, 74, A35.
Herrmann, F. J., and G. Hennenfent, 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, 233–248.
Herrmann, F. J., P. P. Moghaddam, and C. C. Stolk, 2008, Sparsity- and continuitypromoting seismic imaging with curvelet frames: Journal of Applied and Computational Harmonic Analysis, 24, 150–173. (doi:10.1016/j.acha.2007.06.007).
Krohn, C., and R. Neelamani, 2008, Simultaneous sourcing without compromise: Rome 2008, 70th EAGE Conference & Exhibition, B008.
Levy, S., D. Oldenburg, and J. Wang, 1988, Subsurface imaging using magnetotelluric data: Geophysics, 53, 104–117.
Lin, T., and F. J. Herrmann, 2009a, Designing simultaneous acquisitions with compressive sensing: Presented at the EAGE Technical Program Expanded Abstracts, EAGE.
——–, 2009b, Unified compressive sensing framework for simultaneous acquisition with primary estimation: SEG Technical Program Expanded Abstracts, SEG, 3113–3117.
Malcolm, A. E., M. V. de Hoop, and J. A. Rousseau, 2005, The applicability of dmo/amo in the presence of caustics: Geophysics, 70, S1–S17.
Mallat, S. G., 2009, A Wavelet Tour of Signal Processing: the Sparse Way: Academic Press.
Moldoveanu, N., 2010, Random sampling: A new strategy for marine acquisition: SEG, Expanded Abstracts, 29, 51–55.
Oldenburg, D. W., S. Levy, and K. P. Whittall, 1981, Wavelet estimation and deconvolution: Geophysics, 46, 1528–1542.
Rauhut, H., K. Schnass, and P. Vandergheynst, 2008, Compressed sensing and redundant dictionaries: IEEE Trans. Inform. Theory, 54, 2210–2219.
Romberg, J., 2009, Compressive sensing by random convolution: SIAM Journal on Imaging Sciences, 2, 1098–1128.
Sacchi, M. D., T. J. Ulrych, and C. Walker, 1998, Interpolation and extrapolation using a high-resolution discrete Fourier transform: IEEE Trans. Signal Process., 46, 31–38.
Sacchi, M. D., D. R. Velis, and A. H. Cominguez, 1994, Minimum entropy deconvolution with frequency-domain constraints: Geophysics, 59, 938–945.
Santosa, F., and W. Symes, 1986, Linear inversion of band-limited reflection seismogram: SIAM J. of Sci. Comput., 7.
Smith, H. F., 1998, A Hardy space for Fourier integral operators: J. Geom. Anal., 8, 629–653.
Spitz, S., 1999, Pattern recognition, spatial predictability, and subtraction of multiple events: The Leading Edge, 18, 55–58.
Tang, G., J. Ma, and F. Herrmann, 2009, Design of two-dimensional randomized sampling schemes for curvelet-based sparsity-promoting seismic data recovery: Technical Report TR-2009-03, the University of British Columbia.
Taylor, H. L., S. Banks, and J. McCoy, 1979, Deconvolution with the 1 norm: Geophysics, 44, 39.
Trad, D., T. Ulrych, and M. Sacchi, 2003, Latest views of the sparse radon transform: Geophysics, 68, 386–399.
Trad, D. O., 2003, Interpolation and multiple attenuation with migration operators: Geophysics, 68, 2043–2054.
Ulrych, T. J., and C. Walker, 1982, Analytic minimum entropy deconvolution: Geophysics, 47, 1295–1302.
Xu, S., Y. Zhang, D. Pham, and G. Lambare, 2005, Antileakage Fourier transform for seismic data regularization: Geophysics, 70, V87 – V95.
Zwartjes, P. M., and M. D. Sacchi, 2007a, Fourier reconstruction of nonuniformly sampled, aliased data: Geophysics, 72, V21–V32.
——–, 2007b, Fourier reconstruction of nonuniformly sampled, aliased seismic data: Geophysics, 72, V21–V32.