# Geophysics

- © 2008 Society of Exploration Geophysicists. All rights reserved.

## Abstract

We present a new, discrete undersampling scheme designed to favor wavefield reconstruction by sparsity-promoting inversion with transform elements localized in the Fourier domain. The work is motivated by empirical observations in the seismic community, corroborated by results from compressive sampling, that indicate favorable (wavefield) reconstructions from random rather than regular undersampling. Indeed, random undersampling renders coherent aliases into harmless incoherent random noise, effectively turning the interpolation problem into a much simpler denoising problem. A practical requirement of wavefield reconstruction with localized sparsifying transforms is the control on the maximum gap size. Unfortunately, random undersampling does not provide such a control. Thus, we introduce a sampling scheme, termed *jittered undersampling*, that shares the benefits of random sampling and controls the maximum gap size. The contribution of jittered sub-Nyquist sampling is key in formu-lating a versatile wavefield sparsity-promoting recovery scheme that follows the principles of compressive sampling. After the behavior of the jittered-undersampling scheme in the Fourier domain is analyzed, its performance is studied for curvelet recovery by sparsity-promoting inversion (CRSI). The findings on synthetic and real seismic data indicate an improvement of several decibels over recovery from regularly undersampled data for the same amount of data collected.

## INTRODUCTION

The argument has been made that no real theoretical requirement exists for regular spatial sampling of seismic data (Bednar, 1996). However, popular multitrace processing algorithms, e.g., surface-related multiple elimination (SRME) (Verschuur et al., 1992) and wave-equation migration (WEM) (Claerbout, 1971), need a dense and regular coverage of the survey area. Nevertheless, field data sets typically are sampled irregularly and/or coarsely along one or more spatial coordinates and need to be interpolated before being processed.

For regularly undersampled data along one or more spatial coordinates, i.e., data spatially sampled below the Nyquist rate, a wide variety of wavefield-reconstruction techniques exists. Filter-based methods interpolate by convolution with a filter designed so that the error is white noise. The most common of these filters are prediction-error filters (PEFs), which can handle aliased events (Spitz, 1991). Methods based on wavefield operators represent another type of interpolation approach that explicitly includes wave propagation (Canning and Gardner, 1996; Biondi et al., 1998; Stolt, 2002). Transform-based methods also provide efficient algorithms for seismic data regularization (Sacchi et al., 1998; Trad et al., 2003; Herrmann and Hennenfent, 2008; Zwartjes and Sacchi, 2007). However, for irregularly sampled data, e.g., binned data with some of the bins that are empty, or continuous random, undersampled data, the performance of most of these interpolation methods deteriorates.

This article demonstrates that irregular/random undersampling is not a drawback for particular transform-based interpolation methods and for many other advanced processing algorithms, as observed by other workers (Zhou and Schuster, 1995; Sun et al., 1997; Trad and Ulrych, 1999; Xu et al., 2005; Abma and Kabir, 2006; Zwartjes and Sacchi, 2007). We explain why random undersampling is an advantage and how it can be beneficial when designing coarse sampling schemes. To keep the discussion as clear and concise as possible, we focus on regular sampling with randomly missing data points, i.e., discrete random (under)sampling. However, our conclusions extend to continuous random undersampling. Unless otherwise specified, the term *random* is used in the discrete sense.

### Motivation

Recent results in information theory and approximation theory establish that a signal can be recovered exactly from severely undersampled data points if (1) the signal exhibits sparsity in a known transform domain, (2) the artifacts introduced by undersampling look like incoherent random noise in the sparsifying domain, and (3) a data-consistent, sparsity-promoting procedure is used for recovery. It is possible to build an intuitive understanding of these theoretical results, termed *compressive sampling* (CS) (Candès et al., 2006; Candès and Romberg, 2006; Donoho, 2006), by considering a simple example.

Figure 1a shows the superposition of three cosine functions. This signal is sparse in the Fourier domain (condition 1) and is sampled regularly above the Nyquist rate. Its amplitude spectrum is plotted in Figure 1b. When the signal is undersampled randomly threefold according to a discrete uniform distribution as in Figure 1c, its amplitude spectrum, plotted in Figure 1d, is corrupted by artifacts that look like additive incoherent random noise (condition 2). In this case, the significant coefficients of the to-be-recovered signal remain above the noise level. These coefficients can be detected with a denoising technique that promotes sparsity, e.g., nonlinear thresholding (dashed line in Figure 1d and 1f), and recovered exactly by an amplitude-matching procedure to fit the acquired data (condition 3). This experiment illustrates a favorable recovery from severely undersampled data points of a signal that is sparse in the Fourier domain.

When the original signal is undersampled regularly threefold (Figure 1e), the undersampling artifacts coherently interfere, giving rise to well-known aliases that look like the original signal components (Figure 1f). In this case, the sparsity-promoting recovery scheme might fail because the to-be-recovered signal components and the aliases are both sparse in the Fourier domain. This example suggests that random undersampling according to a discrete uniform distribution is more favorable than regular undersampling for reconstruction algorithms that promote sparsity in the Fourier domain.

In general terms, these observations hint at undersampling schemes that lead to more favorable recovery conditions. Significant advances have been made in compressive sampling regarding the main ingredients that go into the design of an undersampling scheme that favors sparsity-promoting recovery. In this article, we draw on these results to design a new coarse spatial sampling scheme for seismic data.

### Main contributions

We propose and analyze a coarse spatial sampling scheme termed *jittered undersampling* (Leneman, 1966; Dippe and Wold, 1992). It creates, under specific conditions, a favorable recovery situation for seismic wavefield-reconstruction methods that impose sparsity in Fourier or Fourier-related domains (see Sacchi et al., 1998; Xu et al., 2005; Herrmann and Hennenfent, 2008; Zwartjes and Sacchi, 2007).

Jittered undersampling differentiates itself from random undersampling according to a discrete uniform distribution — which also creates favorable recovery conditions (Xu et al., 2005; Abma and Kabir, 2006; Zwartjes and Sacchi, 2007) — by controlling the maximum gap in the acquired data. This control makes jittered undersampling very well suited to methods that rely on transforms with localized elements, e.g., windowed Fourier or curvelet transform (Candès et al., 2005a and references therein). These methods are vulnerable to gaps in the data that are larger than the spatiotemporal extent of the transform elements (Trad et al., 2005).

### Outline

After a brief overview of the CS framework and the criteria for a favorable recovery, we study the effects of different undersampling schemes for signals that are sparse in the Fourier domain. Next, we discuss the advantages of random undersampling and design our jittered-undersampling strategy, which increases control on the acquisition grid. The performance of this new scheme for curvelet-based recovery is illustrated on synthetic and real data.

## THEORY

### Basics of compressive sampling

As mentioned before, CS relies on a sparsifying transform for the to-be-recovered signal and uses this sparsity prior to compensate for undersampling during the recovery process. For reconstructing wavefields in the Fourier (Sacchi et al., 1998; Xu et al., 2005; Zwartjes and Sacchi, 2007), Radon (Trad et al., 2003), and curvelet (Hennenfent and Herrmann, 2005; Herrmann and Hennenfent, 2008) domains, sparsity promotion is a well-established technique documented in the geophysical literature. The main contribution of CS is the new light shed on the favorable recovery conditions.

#### Recovery by sparsity-promoting inversion

Consider the following linear forward model for the recovery problem
**S**. equation 1 now can be reformulated as
**A** when estimating **y**.

After sparsity-promoting inversion, the recovered signal is given by **x**.

Minimizing the **x**, and the equality constraint ensures that the solution honors the acquired data. Among all possible solutions of the severely underdetermined system of linear equations

#### Favorable recovery conditions

Following Verdu (1998) and Donoho et al. (2006), we define the matrix **I** is the identity matrix, and the parameter *α* is a scaling factor such that **z** are referred to as multiple-access interference (MAI).

According to the CS theory (Candès et al., 2006; Donoho, 2006), the solution **z** does not contain coherent energy. The first condition of sparsity requires that the energy of **S** in conjunction with the restriction operator **R**.

Intuitively, it requires that the artifacts **z** introduced by undersampling the original signal **S** domain. When this condition on **z** is not met, sparsity alone is no longer an effective prior to solve the recovery problem. Albeit qualitative, the second condition provides a fundamental insight in choosing undersampling schemes that favor recovery by sparsity-promoting inversion.

### Fourier-domain undersampling artifacts

Undersampling artifacts in the Fourier domain are studied for two reasons. First, several interpolation methods are based on the Fourier transform (Sacchi et al., 1998; Xu et al., 2005; Zwartjes and Sacchi, 2007). Second, the curvelet transform, a dyadic-parabolic partition of the Fourier domain, forms the basis of our recovery scheme (Herrmann and Hennenfent, 2008). In many situations, curvelets are preferred over Fourier because of their ability to represent complex seismic data sparsely. For a detailed discussion on this topic, see Candès et al. (2005a) and Hennenfent and Herrmann (2006).

In the coming discussion, the sparsifying transform is defined as the Fourier transform, i.e., **L** are known as spectral leakage (Xu et al., 2005).

#### Regular (under)sampling

When **R** keeps all of the data points of

When **R** corresponds to a regular undersampling scheme,

In the seismic community, difficulties with regularly undersampled data are acknowledged when reconstructing by promoting sparsity in the Fourier domain. For example, Xu et al. (2005, p. V94) write that the antileakage Fourier transform for seismic data regularization “may fail to work when the input data has severe aliasing.”

#### Random undersampling according to a discrete uniform distribution

When **R** corresponds to a random undersampling according to a discrete uniform distribution, the situation is completely different. The matrix **L** is a random matrix (Figure 2f). Consequently, we have
**n**.

For infinitely large systems (Donoho et al., 2006), this approximation becomes an equality. Because of this property, the recovery problem turns into a much simpler denoising problem, followed by a correction for the amplitudes. Remember that the acquired data **y** are noise free (see equation 2) and that the noise **n** in equation 4 comes only from the underdeterminedness of the system. In other words, random undersampling according to a discrete uniform distribution spreads the energy of spectral leakage across the Fourier domain, turning the noise-free underdetermined problem (see equation 2) into a noisy well-determined problem (see equation 4) whose solution can be recovered by solving equation 3. This observation first was reported by Donoho et al. (2006).

#### The practical requirement of maximum gap control

As shown, random undersampling according to a discrete uniform distribution creates favorable recovery conditions for a reconstruction procedure that promotes sparsity in the Fourier domain. However, a global transform such as the Fourier transform typically does not permit a sparse representation for complex seismic wavefields (Hennenfent and Herrmann, 2006). It requires a more local transform, e.g., windowed Fourier (Zwartjes and Sacchi, 2007) or curvelet (Herrmann and Hennenfent, 2008) transform. In this case, problems arise with gaps in the data that are larger than the spatiotemporal extent of the transform elements (Trad et al., 2005). Consequently, undersampling schemes with no control on the size of the maximum gap, e.g., random undersampling according to a discrete uniform distribution, become less attractive.

The term *gap* refers to the interval between two adjacent acquired traces minus the interval associated with the fine interpolation grid, such that adequate sampling has gaps of zero. We present an under-sampling scheme that has an antialiasing effect under some specific conditions yet offering control on the size of the maximum gap.

### Uniform jittered undersampling on a grid

First, the undersampling grid is defined for a discrete uniform jitter. Next, the spectral leakage caused by this scheme is studied.

#### Definition of jittered grid

The basic idea of jittered undersampling is to decimate regularly the interpolation grid and perturb subsequently the coarse-grid sample points on the fine grid. As for random undersampling according to a discrete uniform distribution, in which each location is equally likely to be sampled, a discrete uniform distribution for the perturbation around the coarse-grid points is considered (see Appendix A and Leneman [1966] for more details).

To keep the derivation of our jittered-undersampling scheme succinct, the undersampling factor *γ* is taken to be odd, i.e., *N* of the interpolation grid is a multiple of *γ* so that the number of acquired data points *q*, denoted *q*. The above sampling can be adapted when *γ* is even.

Figure 3 illustrates samplings with increasing randomness. The fine grid of open circles denotes the interpolation grid on which the model

As mentioned earlier, recovery with localized transforms depends on both the maximum gap size and a sufficient sampling randomness to break the coherent aliases. In the next section, we show how the value of the jitter parameter controls these two aspects in our undersampling scheme.

#### Fourier-domain artifacts of the jittered grid

When **R** describes a jittered-undersampling scheme according to a discrete uniform distribution, the stochastic expectation E{⋅} of the first column **a** of the circulant matrix

Equation 6 corresponds to an elementwise multiplication of the periodic Fourier spectrum of the discrete regular sampling vector with a sinc function. This sinc function follows from the Fourier transform of the probability density function for the perturbation with respect to a point of the regularly decimated grid.

In Figure 4, the amplitudes for this Fourier-domain multiplication are plotted for jittered undersampling with

Equation 6 is a special case of the result for jittered undersampling according to an arbitrary distribution introduced by Leneman (1966) and further detailed in Appendix A. Because these results originally were derived for the continuous case, equation 6 is approximate. In practice, however, this formula proves to be accurate, an observation corroborated by numerical results presented below. Consider the following cases for a fixed undersampling factor *γ*.

##### Regular undersampling ( ξ = 0 )

As observed from the first row of Figure 3, there is no jitter in this case and equation 6 becomes
**z** consist of aliased energy.

##### Optimally jittered undersampling ( ξ = γ )

Now the sampling points are perturbed within contiguous windows, as depicted in the third row of Figure 3, and equation 6 reduces to

As with random undersampling according to a discrete uniform distribution, the off-diagonals of the matrix **L** does not contain coherent energy, as observed in Figure 5d, for a fivefold undersampling

##### Jittered undersampling ( 0 < ξ < γ )

In this regime, both coherent aliases and incoherent random undersampling noise are present. Depending on the choice for the jitter parameter, the energy either localizes or randomly spreads across the spectrum. Again, the reduction of the aliases is related to the locations of the zero crossings of the sinc function that move as a function of *ξ*. As *ξ* increases, the zeros move closer to the aliases.

As expected, **L** (Figure 5c) that is a superposition of coherent aliases and incoherent random noise. Although this regime reduces the aliases, coherent energy remains in the undersampling artifacts. This residue creates a situation that is less favorable for recovery. Depending on the relative strength of the aliases compared with the magnitude *n* of the diagonal of

In the next section, a series of controlled experiments is conducted to compare the recovery from regularly, randomly according to a discrete uniform distribution, and optimally jittered undersamplings.

### Controlled recovery experiments for different sampling schemes

With the favorable sampling identified, it remains to be shown that these samplings lead to an improved recovery compared with the unfavorable regular undersampling. In particular, we want to confirm experimentally that jittered undersampling behaves as random undersampling according to a discrete uniform distribution.

For this purpose, we define the sparsifying transform **S** as the Fourier transform **F**, i.e., *k* nonzero entries and of length **y** are obtained by undersampling *k* of nonzero entries of

Each experiment is repeated 100 times for the different undersampling schemes and for varying undersampling factors *γ*, ranging from two to six. The reconstruction error is the number of entries at which the estimated representation

The averaged results for the different experiments are summarized in Figure 6a 6b 6c, which corresponds to regular, random, and optimally jittered undersampling, respectively. The horizontal axes in these plots represent the relative underdeterminedness of the system, i.e., the ratio of the number *k* of nonzero entries in *n* of acquired data points. The vertical axes represent the average reconstruction error. The different curves represent the different undersampling factors. In each plot, the curves from top to bottom correspond to

Figure 6a shows that regardless of the undersampling factor, there is no range of relative underdeterminedness for which

The key observations from these experiments are threefold. First, it is possible, under specific conditions, to recover exactly by sparsity-promoting inversion the original spectrum

## APPLICATION TO SEISMIC DATA

Following recent work on curvelet reconstruction with sparsity-promoting inversion (CRSI), (Herrmann and Hennenfent, 2008), seismic wavefields are reconstructed via **C** is the discrete wrapping-based curvelet transform (Candès et al., 2005a).

Similarly to any other data-independent transforms, curvelets do not provide a sparse representation of seismic data in the strict sense. Instead, the curvelet transform provides a compressible representation—arguably the most compressible (Hennenfent and Herrmann, 2006). Compressibility means that most of the wavefield energy is captured by a few significant coefficients in the sparsifying domain. Because CS guarantees, for sparse enough signal representations, the recovery of a fixed number of largest coefficients for a given undersampling factor (Candès et al., 2005b), a more compressible representation yields a better reconstruction, which explains the success of CRSI.

### Synthetic data example

Figure 7a shows a synthetic data set sampled above the Nyquist rate along both the time and receiver axes. The corresponding amplitude spectrum is plotted in Figure 7b. These two figures serve as references. Comparisons are made between the interpolation results of threefold spatially undersampled data, collected either regularly or optimally jittered.

As expected, the amplitude spectrum (Figure 8c) of the regularly undersampled data (Figure 8a) is severely aliased. Unfortunately, these coherent

Figure 9 shows the CRSI results for these two experiments. Figure 9a and 9b depicts the reconstructions given data sampled regularly (Figure 8a) and optimally jittered (Figure 8b), respectively. Figure 9c and 9d represents the corresponding amplitude spectra. Unlike Figure 9d, which is corrupted only slightly by incoherent errors, Figure 9c still contains substantial energy from the coherent undersampling artifacts. This observation is corroborated by the respective signal-to-reconstruction-error ratios of 6.91 and

It is important to remember that the difference in reconstruction quality is solely because of the difference in spatial sampling; the undersampling factor and the recovery procedure were kept the same. This behavior leads us to conclude that for a given undersampling factor, spatial optimally jittered undersampling is much more favorable for CRSI than regular undersampling.

In addition, Figure 10 shows a recovery experiment given randomly threefold spatially undersampled data. Figure 10a depicts the simulated acquired data, and Figure 10b shows the CRSI result. The signal-to-reconstruction-error ratio is

### Field data example

The far offsets of a regularly sampled shot taken from a real marine data set are considered. Our model consists of 255 traces separated by

## DISCUSSION

### Undersampled data contaminated by noise

Although we focused on a noise-free severely underdetermined system of linear equations, the CS theory and hence our work both extend to the recovery from undersampled data contaminated by noise (Candès et al., 2005b). In this case, the noise **e** that corrupts the data adds to the undersampling artifacts in the sparsifying domain. The quantity that relates to the recoverability is now given by **z** and the imprint of the contaminating noise in the sparsifying domain, i.e.,

### From discrete to continuous spatial undersampling

So far, undersampling schemes based on an underlying fine interpolation grid have been considered. This situation typically occurs when binning continuous randomly sampled seismic data into small bins that define the fine grid used for interpolation. Despite the error introduced in the data, binning presents some computational advantages because it allows for the use of fast implementations of Fourier or Fourier-related transforms, e.g., fastest Fourier transform in the West (FFTW) (Frigo and Johnson, 1998) or fast discrete curvelet transform (FDCT), (Candès et al., 2005a). However, binning can lead at the same time to an unfavorable undersampling scheme, e.g., regular or poorly jittered. In this case, one should consider working on the original data with, for example, an extension to the curvelet transform for irregular grids (Hennenfent and Herrmann, 2006). Despite the extra computational cost for the interpolation, continuous random sampling typically improves interpolation results because it does not create coherent undersampling artifacts (Xu et al., 2005).

### Sparsity-promoting solvers and jittered undersampling

The applicability of CS to the large-scale problems of exploration geophysics relies heavily on the implementation of an efficient

The success and/or efficiency of these approximate solvers depends on the implicit or explicit assumption that the MAI is incoherent. Because optimally jittered undersampling creates such an MAI, these solvers can be used for the sparsity-promoting reconstruction with curvelets or other localized Fourier-based transforms. More important, jittered undersampling can be useful to evaluate the efficiency and robustness of approximate

### Generalization of the concept of undersampling artifacts

Undersampling artifacts are only one particular case of MAI that specifically occurs in the interpolation problem, i.e., **R** can be extended to more general cases (see, e.g., Lustig et al., 2007 in magnetic resonance imaging). For example, when **A** is defined as **M** a modeling/demigration operator (Herrmann et al., 2008; Wang and Sacchi, 2007), **S** domain, and **y** is the incomplete seismic data. The study of the MAI now determines which coarse spatial sampling schemes are more favorable than others in the context of sparsity-promoting migration/inversion. Based on observations in Zhou and Schuster (1995) and in Sun et al. (1997), we believe that discrete random, optimally jittered, and continuous random undersamplings also will play a key role.

## CONCLUSIONS

Successful wavefield recovery depends on three key factors: the existence of a sparsifying transform, a favorable sampling scheme, and a sparsity-promoting recovery method. We have focused on an undersampling scheme designed for localized Fourier-like signal representations such as the curvelet transform. Our scheme builds on the fundamental observation that irregularities in sub-Nyquist sampling are good for nonlinear sparsity-promoting wavefield-reconstruction algorithms because they turn harmful coherent aliases into relatively harmless incoherent random noise. The interpolation problem effectively becomes a much simpler denoising problem in this case.

Undersampling with a discrete random uniform distribution lacks control on the maximum gap size in the acquisition, which causes problems for transforms that consist of localized elements. Our jittered-undersampling scheme remedies this lack of control while preserving the beneficial properties of randomness in the acquisition grid. Our numerical findings on a stylized series of experiments confirm these theoretically predicted benefits.

Curvelet-based wavefield-reconstruction results from jittered-undersampled synthetic and field data sets are better than results obtained from regularly decimated data. In addition, our findings indicate an improved performance compared with traces taken randomly according to a uniform distribution. This is a major result with wide-ranging applications because it entails an increased probability for successful recovery with localized transform elements. In practice, this translates into more robust wavefield reconstruction.

## Acknowledgments

G. H. thanks Ken Bube, Ramesh Neelamani, Warren Ross, Bea-trice Vedel, and Ozgur Yilmaz for constructive discussions about this research. Eric Verschuur and Chevron Energy Technology Company are gratefully thanked for the synthetic and real data sets, respectively. The authors thank the authors of CurveLab (www.curvelet.org) and the authors of SPGL1 (www.cs.ubc.ca/labs/scl/spgl1) for making their codes available. The paper was prepared with Madagascar, a reproducible research package (rsf.sourceforge.net). This work was supported financially in part by NSERC Discovery grant 22R81254 and CRD grant DNOISE 334810-05 of F.J.H. and was carried out as part of the SINBAD project with support, secured through ITF, from BG Group, BP, Chevron, ExxonMobil, and Shell. The authors also appreciate the valuable comments and suggestions from the two reviewers and two associate editors.

## APPENDIX A

### JITTERED UNDERSAMPLING

For the convenience of the reader, we rederive the result originally introduced by Leneman (1966) that leads to equation 6.

Jittered sampling locations *p* on *s* is given by

- Received September 21, 2007.