On sampling functions and Fourier reconstruction methods - CiteSeerX

2 downloads 2735 Views 2MB Size Report
versely, regular sampling often hampers the Fourier data recovery methods. ... Data reconstruction methods based on signal processing concepts have been ...
On sampling functions and Fourier reconstruction methods Mostafa Naghizadeh



and Mauricio D. Sacchi



ABSTRACT Random sampling can lead to algorithms where the Fourier reconstruction is almost perfect when the underlying spectrum of the signal is sparse or band-limited. Conversely, regular sampling often hampers the Fourier data recovery methods. However, two-dimensional (2D) signals which are band-limited in one spatial dimension can be recovered by designing a regular acquisition grid that minimizes the mixing between the unknown spectrum of the well-sampled signal and aliasing artifacts. This concept can be easily extended to higher dimensions and used to define potential strategies for acquisition-guided Fourier reconstruction. In this paper we derive the wavenumber response of various sampling operators and investigate sampling conditions for optimal Fourier reconstruction using synthetic and real data examples.

INTRODUCTION The sampling theorem is an interesting topic of chief importance in the physical sciences and engineering. Sampling methods can be classified into two categories, uniform and nonuniform sampling methods. In uniform sampling methods, a signal is sampled periodically at a constant rate. This leads to recovery conditions based on the well-known Nyquist sampling theorem. An overview of uniform sampling and its properties is given by Unser (2000). For nonuniform sampling methods, a distinction is made between methods that assume known sampling positions, and those where the position at which observations were taken are not precisely known (Vandewalle et al., 2007). In this paper, we will concentrate on methods with known sampling locations. Multi-channel sampling is another form of sampling that leads to the so called superresolution reconstruction techniques (Bertero and Boccacci, 2003; Park et al., 2003). A geophysical example of multi-channel sampling was provided by Ronen (1987) who proposed an algorithm to estimate an alias-free zero-offset section from aliased common offset seismic sections. This algorithm was later developed into inversion to zero-offset (Ronen et al., 1991, 1995). Data reconstruction methods based on signal processing concepts have been proposed to reconstruct seismic data. Two important categories are prediction filtering methods, and transform-based methods. Prediction filtering approaches have been proposed by Spitz (1991); Porsani (1999). The frequency-wavenumber (f-k ) equivalence of prediction filtering methods is introduced by Gulunay (2003). Their techniques are capable of interpolating regularly sampled aliased data. These type of methods have been modified by Naghizadeh and Sacchi (2007) to cope with irregularly sampled data. Transform-based methods map the signal to a new domain and a synthesis of the signal from the transformed domain is used to generate data at un-recorded spatial locations. Transform-based seismic data Fourier reconstruction

2

reconstruction have been proposed by, to cite a few, Bardan (1987), Darche (1990), and Trad (2008). Transform-based methods which operate with Fourier basis have been also adapted for seismic data reconstruction. In this case the seismic signal is represented via a superposition of complex exponentials. Fourier reconstruction methods rely on two basic assumptions. In general, we assume band-limited signals in the wavenumber domain and/or signals that can be represented by a parsimonious distribution of Fourier coefficients. Examples of the aforementioned assumptions abound in the signal and imaging processing literature. For instance, reconstruction of band-limited signals via Fourier methods have been studied by Strohmer (1997); Feuer and Goodwin (2005); Eldar (2006). In the geophysical literature Duijndam et al. (1999) and Schonewille et al. (2003) proposed a band-limited signal reconstruction method for seismic data that depends on two and three spatial dimensions. Methods that used the sparsity assumption to reconstruct data were proposed by Sacchi and Ulrych (1996) and Sacchi et al. (1998) and expanded to the multidimensional case by Zwartjes and Gisolf (2006) and Zwartjes and Sacchi (2007). The intention of this paper is to provide an understanding of the importance of sampling at the time of defining a Fourier reconstruction strategy. In particular we investigate different type of samplings and the impact they have on the accuracy of the general form of band-limited Fourier reconstruction introduced by Liu (2004). The latter, also called Minimum Weighted Norm Interpolation (MWNI) method, permits to retrieve the wavenumber spectrum of desired data from a decimated data. It must be mentioned that this article only investigates the MWNI method. The conclusions arising from our analysis, however, are valid for other reconstruction methods that adopt a simple or sparse spectral representation. This includes the anti-leakage Fourier transform (Xu et al., 2005), projection onto convex sets (Abma and Kabir, 2006) and other sparsity promoting solutions to the interpolation problem (Sacchi et al., 1998; Hennenfent, 2008). Also, our general conclusions may apply to methods that use other transforms such as wavelets (Beylkin et al., 1991), curvelets (Candes, 2006), seislets (Fomel, 2006) and etc. Reconstruction methods that adopt sparsity promoting strategies are part of the general theory of compressive sensing (Donoho, 2006; Tsaig and Donoho, 2006). The latter deals with the optimal recovery conditions for signals that admit sparse representation in surrogate domain. In our paper, random sampling means randomly picked samples from a regular grid of data. A regular grid of data is a requirement for the MWNI method (Liu, 2004). The latter permits to design an computationally efficient interpolator that heavily relies on spatial fast Fourier transforms (FFT). The anti-leakage Fourier transform reconstruction (Xu et al., 2005), on the other hand, uses the exact irregular spatial positions and consequently requires operators based on non-uniform discrete Fourier transforms. It must be mentioned that in real practical surveys the pure irregularity of spatial sampling might completely eliminate the alias problem. However, in most of seismic surveys, the spatial sampling operators are near regular and therefore they might suffer from the same problems of regular grids. This paper is organized as follows. First, we review one-dimensional and two-dimensional sampling operators and their frequency or wavenumber responses. We examine the wavenumber domain footprint of the sampling operators and study cases where MWNI can properly reconstruct data from under-sampled observations. We then provide examples of synthetic and real data, highlighting our finding on the relationship between sampling operators and Fourier reconstruction.

3

SAMPLING THEORY One-dimensional sampling functions We start our discussion by considering a discrete signal x of length N x = {x(0), x(1), x(2), . . . , x(N − 1)}.

(1)

In addition, we will replace N − M arbitrary samples of x by zeros to simulate missing observations. The location of the available samples are indicated by the set H = {h(0), h(1), h(2), . . . , h(M − 1)}. The signal with missing samples is represented via the following expression  x(h(n)) h(n) ∈ H xs (n) = . (2) 0 h(n) ∈ /H The discrete Fourier transform (DFT) of x is defined by X(k) =

N −1 −i2πnk 1 X x(n)e N N

k = 0, 1, 2, . . . , N − 1.

(3)

n=0

We now we compute the DFT of signal with missing samples Xs (k) =

M −1 −i2πh(n)k 1 X x(h(n))e N , N

k = 0, 1, 2, . . . , N − 1,

(4)

n=0

It is easy to show (Appendix) that the DFT of the sampled signal can be represented in vector form as follows Xs = X ∗ Q,

(5)

where the symbol ∗ indicates discrete periodic convolution (Oppenheim and Schafer, 1974) and Q denotes the Fourier response of sampling function with elements given by Q(k) =

M −1 1 X −i2πh(m)k N e N

k = 0, 1, 2, . . . , N − 1.

(6)

m=0

The sampling function in equation 6 convolves the spectrum of the original signal to yield the spectrum of the signal with missing samples. This equation is fundamental to understanding Fourier reconstruction methods. Let’s use examples to study the footprint of the sampling operator in the wavenumber domain. Figure 1 shows sampling functions (left side) and their correspondent frequency responses (right side). In this example we assume that the data has been regularly decimated. The decimation factors are 1, 2, 3, 4 and 5 from top to bottom, respectively. For the decimation factor equal to 1 (no decimation), the wavenumber response of the sampling function has only one impulse at the normalized wavenumber k = 0.0. This means that the DFT vectors Xs and X are equal. When the data are decimated by a factor of 2, every other sample is replaces by zeros and the wavenumber response contains two impulses. In other words, the DFT of the sampled data will contain replicas of the DFT of

4

1 0

1

5

10 15 20 Sample number

25

30

0 −0.5

0 Normalized wavenumber

0.5

Figure 1: Sampling operators (left) and their associated wavenumber spectra (right). The decimation factors are 1, 2, 3, 4 and 5 from top to bottom.

the data prior to decimation. Increasing the decimation factor results in more replications of the original spectrum. It is important to notice that the size of replicated impulses is preserved. If KN denotes the Nyquist wavenumber, the decimation process by an integer factor L of a periodic sequence can also be described as the lateral cyclic rotation of the original wavenumber spectrum by amounts KN /L, 2KN /L, . . . , (L − 1)KN /L followed by summation (Gulunay, 2003). Figure 2 depicts 5 different random sampling scenarios (left side) and their wavenumber responses (right side). For all of the examples, 50% of the samples are missing. The wavenumber response of a nonuniform sampling operator has a different structure than the uniform one. It has nonzero values for the all of its components inside the fundamental interval of the wavenumber spectrum. Interestingly, it has only one large value which appears at the exact location k = 0. The rest of the components are small. As the number of missing samples increases, more components of the sampling function are dominant. Figure 3 shows 5 random sampling operators with 85% missing samples. As the number of missing samples increases, the amplitude of the impulses in the wavenumber response of the random sampling operator increases. To continue with the analysis we now examine the case where the data contain gaps. Figure 4 shows the plot of various sampling functions (left side) and their associated wavenumber responses (right side). Similarly, Figure 5 depicts the plot of different sampling functions (left side) for the extrapolation case and their corespondent wavenumber responses (right side). Both cases produce wavenumber responses with a dominant amplitude responses at the negative and positive normalized wavenumbers near k = 0. The spectrum of the data that has undergone this type of sampling will be blurred by the convolution with the wavenumber response of the sampling function. On the other hand, in the regularly sam-

5

1 0

1

5

10 15 20 Sample number

25

30

0 −0.5

0 Normalized wavenumber

0.5

Figure 2: Sampling operators (left) and their associated wavenumber spectra (right). For each example 50% of data were randomly removed.

1 0

1

5

10 15 20 Sample number

25

30

0 −0.5

0 Normalized wavenumber

0.5

Figure 3: Sampling operators (left) and their associated wavenumber spectra (right). For each example 85% of data were randomly eliminated.

6

1 0

1

5

10 15 20 Sample number

25

0 −0.5

30

0 Normalized wavenumber

0.5

Figure 4: Sampling operators (left) and their associated wavenumber spectra (right) for a gap inside the data. The size of the gap decreases from top to bottom.

pled case the spectrum of the sampled data will exhibit, depending on wavenumber support, organized replicas (alias), which can only be removed when there is no overlap between the spectrum of the signal and the aliased components (band-limited signal). In other words, in order to interpolate the data, such overlaps must be suppressed. (Gulunay, 2003).

Multidimensional sampling functions The formulas introduced in the previous section can be easily generalized to multidimensional cases. We will concentrate on the 2D case but bear in mind that a generalization to N-D case is straightforward. Let us consider the lattice u, with size Nx × Ny represented by  nx ∈ 0 : N x − 1 u(nx , ny ) = . (7) ny ∈ 0 : N y − 1 Assume that we arbitrarily reserve M samples of data and set the rest to zero to obtain us , the signal with missing samples. Let us indicate the location of the available samples by H = {(hx (m), hy (m)) | m ∈ 0 : M − 1},

(8)

where, hx and hy indicate the location of the available samples on the axis X and Y , respectively. The DFT of the original data, U, is given by: Ny −1 1 X U (kx , ky ) = Ny ny =0

Nx −1 kx nx 1 X u(nx , ny )e−i2π Nx Nx nx =0

!

−i2π

e

ky ny Ny

.

(9)

7

1 0

1

5

10 15 20 Sample number

25

30

0 −0.5

0 Normalized wavenumber

0.5

Figure 5: Sampling operators (left) and their associated wavenumber spectra (right) for the extrapolation case. The number of missing samples is decreasing from top to bottom.

Now, the DFT of the data with missing samples, Us , is obtained by Us = U ~ Q,

(10)

where, ~ represents the 2D convolution operator and Q is the 2D Fourier response given by M −1 1 X −i2π kx hNx (m) −i2π ky hNyy(m) x Q(kx , ky ) = e e . (11) Nx Ny m=0

In this paper the absolute value (or amplitude) of the wavenumber response of the 2D sampling function, Q, is plotted for various sampling scenarios. To start the analysis, Figure 6b shows the wavenumber response of a sampling function with no missing samples (Figure 6a). The wavenumber response is a single impulse in the middle of the 2D spectrum. Figure 6c shows a 2D sampling operator that eliminates every other slice of data in the Y direction. Filled circles represent the available samples and crosses show the missing samples. Figure 6d shows the wavenumber response of Figure 6c. The wavenumber response has two impulses located at normalized wavenumbers (kx , ky) = {(0, 0), (0, −0.5)}. Figure 6e shows another 2D decimation function that samples the data in a chessboard pattern. Figure 6f shows the wavenumber response of Figure 6e. In comparison to Figure 6d, the impulse at the normalized wavenumber (0, −0.5) is moved to (−0.5, −0.5). This is an interesting phenomenon and can be efficiently utilized for data reconstruction purposes. Assume that the original 2D signal in Figure 6b is band-limited in the kx axis in the interval [−0.2, 0.2] and full band in the ky axis. Therefore, in Figure 6d the replication impulse produce by the sampling operator is inside the interval where the signal lives. In Figure 6f, however, the impulse resides outside the region containing the signal. Hence,

8

in the case of Figure 6f one use a band-limiting constraint on the kx axis to remove the sampling artifacts. Conversely, in Figure 6d band-limitation cannot be used to eliminate the sampling artifacts. Figure 7a shows a sampling function which has eliminated two Y slices between each available Y slice. Figure 7c shows another sampling function which also has eliminated two samples between available samples and the available samples are located in all of the Y and X slices. Figures 7b and 7d show the wavenumber responses of the Figures 7a and 7c. In Figure 7d the band-limitation information on the kx axis can be used to eliminate the artifacts. This is not possible in the situation depicted in Figure 7b. To continue with our analysis, we now examine the case where samples from the grid are randomly eliminated. Figures 8a and 8c show two 2D random sampling operators with 50% and 85% missing samples, respectively. Figures 8b and 8d depict the wavenumber responses of 8a and 8c, respectively. It is interesting that in both wavenumber responses the amplitude of impulses are very small compared to the impulse of the original data located in the center of the plots. Therefore, for randomly sampled multidimensional data with a high percentage of missing samples, the spectral distortion could be minimal. Figure 8e shows a 2D sampling operator for a gap inside the data. The wavenumber response of 2D gap sampling (Figure 8f) contains relatively high amplitude impulses only in the vicinity of the main original impulse.

Band-limited Minimum Weighted Norm Interpolation In the previous section we have reviewed basic concepts in sampling theory. In this section we examine how one can reconstruct seismic data using MWNI for different sampling scenarios. In particular, we will highlight the importance of understanding the interplay between sampling and a given reconstruction algorithm for a successful reconstruction of seismic data. The MWNI method is described in detail in Liu (2004), Liu and Sacchi (2004) and Trad (2009). Interpolation of band-limited data with missing samples using the MWNI method is summarized in the following paragraphs. The observed data dobs is related to the unknown data in the desired grid, d, via the sampling operator G. The desired data can be retrieved by minimizing the following cost function J = ||Gd − dobs ||22 + µ kdk2W ,

(12)

where k.k2W indicates a weighted norm and G is the sampling matrix which maps desired data samples d = d(xn , f ) to available samples dobs = d(xh , f ) at a given frequency f . Its transpose, GT , fills the position of missing samples with zeros (Liu and Sacchi, 2004). The norm is expressed as follows kdk2W =

X D∗ Dk k , Wk2

(13)

k∈K

where Dk indicates the coefficients of the Fourier transform of the vector of spatial data d.

9

X sample number 5 10

a)

15

b)

-0.4

-0.2

Kx 0

0.2

0.4

-0.4

-0.2

Kx 0

0.2

0.4

-0.4

-0.2

Kx 0

0.2

0.4

5

-0.2 Ky

Y sample number

-0.4

0

10 0.2 15

0.4

X sample number 5 10

c)

15

d)

5

-0.2 Ky

Y sample number

-0.4

0

10 0.2 15

0.4

X sample number 5 10

e)

15

f)

5

-0.2 Ky

Y sample number

-0.4

0

10 0.2 15

0.4

Figure 6: a) Full 2D sampling operator. c) 2D regular sampling operator with a decimation of factor 2 in the Y direction. e) 2D regular sampling operator with chessboard pattern decimation. b), d), and f) are 2D wavenumber spectra of (a), (c), and (e), respectively.

10

X sample number 5 10

a)

15

b)

-0.4

-0.2

Kx 0

0.2

0.4

-0.4

-0.2

Kx 0

0.2

0.4

-0.2

5 Ky

Y sample number

-0.4

0

10 0.2

0.4

15

X sample number 5 10

c)

15

d)

-0.2

5 Ky

Y sample number

-0.4

0

10 0.2

15

0.4

Figure 7: a) and c) Two different 2D regular sampling functions with a decimation factor of 3. b) and d) are the wavenumber spectra of (a) and (c), respectively.

11

X sample number 5 10

a)

15

b)

-0.4

-0.2

Kx 0

0.2

0.4 1.0 0.8

5

-0.2 0.6 Ky

Y sample number

-0.4

0

10

0.4

15

X sample number 5 10

c)

15

0.2

0.2

0.4

0

d)

-0.4

-0.2

Kx 0

0.2

0.4 1.0 0.8

5

-0.2 0.6 Ky

Y sample number

-0.4

0

10

0.4

15

X sample number 5 10

e)

15

0.2

0.2

0.4

0

f)

-0.4

Kx 0

0.2

0.4 1.0 0.8

5

-0.2 0.6 Ky

Y sample number

-0.4

-0.2

0 0.4

10 0.2 15

0.4

0.2 0

Figure 8: a) and c) Two 2D random sampling operators with 50% and 85% of missing samples, respectively. e) 2D gap sampling operator. b), d), and f) are 2D wavenumber spectra of (a), (c), and (e), respectively. The second largest amplitudes in (b) and (d) are 0.15 and 0.33, respectively.

12

For the MWNI method the diagonal matrix is defined as  Wk2 k ∈ K Υk = , 0 k∈ /K

(14)

where K indicates the region of support of the Fourier transform. The pseudoinverse of Υ is defined as:  Wk−2 k ∈ K, † Υk = (15) 0 k∈ / K. For the band-limited MWNI, the values of Wk are equal to one, while for the MWNI method, their values must be iteratively updated to find an optimal reconstruction. The minimizer of the cost function in expression 12 is given by ˆ = FH ΥFGT (GFH ΥFGT + αI)−1 dobs , d

(16)

where F is the Fourier transform operator in matrix form, α is trade-off parameter, I is the identity matrix, T and H stand for transpose and Hermitian transpose operators, respectively. For further details see Liu (2004) and Liu and Sacchi (2004). Tables 1 3 demonstrate the Conjugate Gradients (CG) algorithm used to implement the MWNI method. The symbols ~ and represent the circular convolution and element by element multiplication, respectively.

EXAMPLES Synthetic 1D examples Random sampling operator In order to examine the performance of the MWNI method with various sampling operators, a synthetic and real-valued 1D signal with two harmonics was created (Figure 9a). The data contains 120 samples and the harmonics are at 0.06 and 0.13 normalized wavenumbers with amplitudes equal to one. For our 1D examples, the left panels will show the data in the spatial domain (where we have assumed ∆x = 1) and the right panels will represent their correspondent normalized wavenumber domains. To proceed with the first example, 80 percent of the data was eliminated randomly. The missing samples (Figure 9c) were reconstructed using the MWNI method. Figure 9e shows the result of the reconstruction. Since available samples were picked by a random nonuniform sampling operator the reconstruction was successful. This is due to the fact that the wavenumber response for the nonuniform sampling operator has only one high amplitude impulse (Figure 2). Hence, its convolution with the spectrum of the original data (which is simple as well) will create a simple signal (with impulses at the same location of the impulses of the original spectrum) suitable to be reconstructed using the MWNI method. Figures 9b, 9d, and 9f show the Fourier panels of Figures 9a, 9c, and 9e, respectively. There is a common misinterpretation pertaining random sampling, which needs to be clarified. One might think that the shortest distance between consecutive samples (or average sampling rate) determines the success of a reconstruction method based on sparsity or MWNI. Figure 10a shows an example of randomly picked samples in which the shortest

13

Table 1: Conjugate Gradients algorithm for the MWNI method Input: Total number of data samples: nt Indices of the available samples: hns×1 Values of the available samples: yns×1 Hanning window: Hnh×1 Number of external and internal iterations: niter1 , niter2 Error reduction toleration: tol1 , tol2 Output: ˆ Reconstructed data: x begin  1 Desired wavenumbers Initiate band-limiting function : Υnt×1 = 0 Elsewhere for i = 1 to niter1 X = 0nt×1 (Initiation with zero vector) s=y R = Forward Operator(s, h, Υ, nt) (Table 2) P=R 1 = R · R∗ δ1 = s · s∗ for j = 1 to niter2 q = Adjoint Operator(P, h, Υ) (Table 3) α = 1 /(q · q∗ ) X = X + αP s = s − αq R = Forward Operator(s, h, Υ, nt) (Table 2) 2 = 1 1 = R · R∗ β = 1 /2 P = R + βP δ2 = s · s∗ if [(δ1 − δ2 ) < tol1 ] or [δ2 < tol2 ] then Break end δ1 = δ2 end Update and smooth band-limiting function:Υ = abs(X) ~ H end ˆ = IF F T (X) x end

14

Table 2: Forward operator of CG algorithm for the MWNI method. begin t = 0nt×1 t(h) = s T = F F T (t) R=Υ T end

Table 3: Adjoint operator of CG algorithm for the MWNI method. begin T=Υ P t = IF F T (T) qns×1 = t(h) end

a)

b) Amplitude

Amplitude

1 0 -1 10

40

70

5 -0.4

-0.2

0

0.2

0.4

-0.4

-0.2

0

0.2

0.4

d) Amplitude

1 Amplitude

15

100

c) 0 -1 10

40

70

25 15 5

100

e)

f) Amplitude

1 Amplitude

25

0 -1 10

40 70 Sample number

100

25 15 5 -0.4

-0.2 0 0.2 0.4 Normalized wavenumber

Figure 9: a) Original data. c) Data with 80% randomly missing samples. e) Reconstructed data via MWNI. b), d), and f) are the wavenumber spectra of (a), (c), and (e), respectively.

15

a)

b) Amplitude

Amplitude

1 0 -1 10

40

70

15 5

100

c)

-0.4

-0.2

0

0.2

0.4

d) Amplitude

1 Amplitude

25

0 -1 10

40 70 Sample number

100

25 15 5 -0.4

-0.2 0 0.2 0.4 Normalized wavenumber

Figure 10: a) Randomly sampled data with the restriction of a minimum 4 samples distance between consecutive available samples. c) Reconstructed data via MWNI. b) and d) are the wavenumber spectra of (a) and (c), respectively.

distance between consecutive samples was forced to be more than 4 samples. The reconstructed data using the MWNI method (Figure 10c) shows that the signal was reconstructed successfully even though we did not have samples with less than 4 samples apart. This reflects the fact that it is the main underlying grid of the sampling operator that determines the success of the MWNI method. Figures 10b and 10d show the Fourier panels of Figures 10a and 10c, respectively. It is important to mention that for random sampling even relatively high amplitude artifacts can be eliminated using the MWNI method. This is due to the fact that the MWNI algorithm is a non-linear interpolation method which iteratively updates its weighting function in the wavenumber domain by comparing the latest result of the reconstruction to the available samples. Hence, if there were high amplitude artifacts in the spectrum that do not match the available samples, they are eliminated by iterative updating. Unfortunately, artifacts in the spectrum due to regular decimation in the spatial domain will completely match the available samples and therefore, non-linear updating can not eliminate them. Figure 11 shows the relationship between the percentage of available samples and the total reconstruction error for the MWNI method when random sampling is adopted to the original data in Figure 9a. For each percentage of missing samples, 20 different random sampling realizations were performed and the mean of reconstruction error (normalized RMS difference) and its standard deviation were plotted. As the number of available samples decreases, the reconstruction error and its variance increases.

16

0.4

Error

0.3

0.2

0.1

0 0

10

20

30

40 50 60 70 Percentage of avialable samples

80

90

100

Figure 11: Reconstruction error vs percentage of available samples for 1D MWNI of randomly sampled data. For each case 20 different random realizations were carried out. The variance is shown by the error bars.

Regular decimation For the next example (Figure 12a), the signal was decimated by a decimation factor of 2 (every other sample). Due to the regularity of the samples, the MWNI method was not able to reconstruct the missing samples (Figure 12c). While MWNI fails to recover the missing data, one can utilize band-limitation in the Fourier domain for reconstruction purposes. Figure 12e shows the successful reconstruction of the data using the band-limited MWNI method. The band-limiting function in the Fourier domain is designed to eliminate any event out of the normalized frequencies [−0.15, 0.15]. The interval is chosen based on the a priori information regarding the spectrum of the original data. Figures 12b, 12d, and 12f show the Fourier panels of Figures 12a, 12c, and 12e, respectively. While MWNI with the addition of band-limitation can eliminate high amplitude artifacts caused by regular sampling operators, it will fail to do so if the artifacts are mixed with the original spectrum of the data. Figure 13a show a regularly decimated signal with decimation factor of 6. Figure 13b shows the Fourier panel of Figure 13a. It is evident that parts of the artificial events introduced by the regular sampling operator are interfering with the original spectrum of the data. Figure 13c and 13d show the reconstructed data using the band-limited MWNI method and its Fourier panel, respectively. The reconstructed data are different from the original data due to the presence of high amplitude artifacts in the band selected for reconstruction. Combination of regular and random sampling operators A sampling operator can be a combination of regular and random sampling functions. This means that first, a discrete signal is decimated using a regular sampling operator and then the resultant sampled function is sampled using a random sampling operator. Figure 14a shows an example where the available samples have been randomly picked from an already decimated original signal by a factor of 3. Figure 14c shows the reconstructed data using the MWNI method. Since the random sampling operator is applied on the already decimated signal, the MWNI method was able to reconstruct the decimated signal, not the original one. Figure 14e shows the reconstructed data using the MWNI method when

17

a)

b) Amplitude

Amplitude

1 0 -1 10

40

70

5 -0.4

-0.2

0

0.2

0.4

-0.4

-0.2

0

0.2

0.4

d) Amplitude

1 Amplitude

15

100

c) 0 -1 10

40

70

25 15 5

100

e)

f) Amplitude

1 Amplitude

25

0 -1 10

40 70 Sample number

25 15 5

100

-0.4

-0.2 0 0.2 0.4 Normalized wavenumber

Figure 12: a) Regularly decimated data (decimation factor of 2). c) Reconstructed data via MWNI. e) Reconstructed data using band-limited MWNI in the normalized frequency interval of [-0.15,0.15]. b), d), and f) are the wavenumber spectra of (a), (c), and (e), respectively.

a)

b) Amplitude

Amplitude

1 0 -1 10

40

70

15 5

100

c)

-0.4

-0.2

0

0.2

0.4

d) Amplitude

1 Amplitude

25

0 -1 10

40 70 Sample number

100

25 15 5 -0.4

-0.2 0 0.2 0.4 Normalized wavenumber

Figure 13: a) Regularly decimated data (decimation factor of 6). c) Reconstructed data using bandlimited MWNI in the normalized frequency interval of [-0.15,0.15]. b) and d) are the wavenumber spectra of (a) and (c), respectively.

18

a)

b) Amplitude

Amplitude

1 0 -1 10

40

70

5 -0.4

-0.2

0

0.2

0.4

-0.4

-0.2

0

0.2

0.4

d) Amplitude

1 Amplitude

15

100

c) 0 -1 10

40

70

25 15 5

100

e)

f) Amplitude

1 Amplitude

25

0 -1 10

40 70 Sample number

100

25 15 5 -0.4

-0.2 0 0.2 0.4 Normalized wavenumber

Figure 14: a) Random sampling from a regularly decimated data (decimation factor of 3). c) Reconstructed data via MWNI. e) Reconstructed data using band-limited MWNI in the normalized frequency interval of [-0.15,0.15]. b), d), and f) are the wavenumber spectra of (a), (c), and (e), respectively.

band-limitation is included. Notice that we have simultaneously utilized band-limitation and the minimum weighted norm constraint to remove the effects of both regular and random sampling operators. Figures 14b, 14d, and 14f show the Fourier panels of Figures 14a, 14c, and 14e, respectively. Gap operator Figure 15a shows an example with a gap inside the data. The reconstructed data using the MWNI method (Figure 15c) shows a successful recovery of the samples in the gap. It is clear that the artifacts caused by the gap in the frequency domain (15b) are all in the vicinity of the main spectrum of the original signal (9b). This was expected since the Fourier response of the gap sampling operator (Figure 4) has high value impulses only in the vicinity of the main impulse of the original spectrum. Therefore, as the Fourier panel of the reconstructed data (Figure 15d) shows, the MWNI method was capable of eliminating the artifacts caused by a gap inside the regularly sampled grid of data.

19

a)

b) Amplitude

Amplitude

1 0 -1 10

40

70

15 5

100

c)

-0.4

-0.2

0

0.2

0.4

d) Amplitude

1 Amplitude

25

0 -1 10

40 70 Sample number

100

25 15 5 -0.4

-0.2 0 0.2 0.4 Normalized wavenumber

Figure 15: a) Data with erasures due to the presence of a gap. c) Reconstructed data via MWNI. b) and d) are the wavenumber spectra of (a) and (c), respectively.

Synthetic 2D examples Random sampling In this section we utilize the MWNI method to reconstruct 2D data. In fact, the 2D MWNI method is carried out in the f-x domain and, therefore, it is implemented as a 1D spatial reconstruction process for each frequency. Figure 16a shows an original synthetic 2D section composed of three linear events. The original data are severely aliased and the band-limited MWNI method can not be used for the frequencies above the normalized frequency 0.12. More than 60 percent of traces are randomly eliminated creating a section with missing traces shown in Figure 16b. Figure 16c shows the reconstructed data using the MWNI method at each frequency. Figures 16d, 16e and 16f show the f-k spectrum of Figures 16a, 16b and 16c, respectively. The random sampling in the spatial domain creates low amplitude artifacts around the spectrum of the original data. This is the behavior expected from the random sampling operator as shown in Figure 2. In the case of random sampling of traces, the MWNI method is capable of interpolating the missing traces (Figure 16c). However, there are some isolated high amplitude artifacts left in the f-k spectrum of the reconstructed data. Regular sampling The same 2D synthetic section was also decimated by a factor of 2 (Figure 17a) in order to be reconstructed by the MWNI method. Since the sampling operator was regular, the MWNI method was not able to recover the missing traces (Figure 17b). Figures 17c and 17d show the f-k spectrum representation of the figures 17a and 17b, respectively. The effect of regular sampling function on the spectrum of the data is the superposition of the

20

10

Trace number 20

30

b)

Trace number 20

30

c)

0.2

0.4

0.4

0.4

Time (s)

0.2

0.6

0.6

0.6

0.8

0.8

0.8

e) -0.4

Normalized wavenumber -0.2 0 0.2 0.4

-0.4

0.2

0.3

0.4

Normalized wavenumber -0.2 0 0.2 0.4

-0.4

0.1

Normalized frequency

0.1

0.2

0.3

0.4

10

Trace number 20

30

f) Normalized wavenumber -0.2 0 0.2 0.4

0.1

Normalized frequency

d)

Normalized frequency

10

0.2

Time (s)

Time (s)

a)

0.2

0.3

0.4

Figure 16: 2D data in the t-x domain. a) Original data. b) Randomly sampled data. c) Reconstructed data via MWNI. d), e), and f) are the f-k spectra of (a), (b), and (c), respectively.

21

10

Trace number 20

30

b)

0.2

0.2

0.4

0.4

Time (s)

Time (s)

a)

0.6

0.6

0.8

0.8

c)

-0.4

Normalized wavenumber -0.2 0 0.2 0.4

d)

0.2

0.3

0.4

-0.4

Trace number 20

30

Normalized wavenumber -0.2 0 0.2 0.4

0.1

Normalized frequency

Normalized frequency

0.1

10

0.2

0.3

0.4

Figure 17: a) Regular decimation of data in Figure 17a. b) Reconstructed data via MWNI. c) and d) are the f-k spectra of (a) and (b), respectively.

cyclic shifted spectrum with the original one creating the wrapping of the signal in the f-k domain (Gulunay, 2003). Therefore, since the original spectrum of data was spatially aliased, the resulted replicated spectrum completely mixes with the original spectrum. Applying the MWNI method separately for each frequency can not discriminate between the original and the replicated spectrum. The replicated spectrum is well-separated in the low frequency portion of the data and a band-limiting (low-pass) function can be deployed in the wavenumber axis to eliminate them. However, high frequency spatial data, the bandlimitation is no longer effective. To overcome this problem, Gulunay (2003) introduced an f-k reconstruction method where a mask function is built from low frequencies to eliminate the aliasing artifacts in the high frequencies.

22

10

Trace number 20

30

b)

0.2

0.2

0.4

0.4

Time (s)

Time (s)

a)

0.6

0.6

0.8

0.8

c)

-0.4

Normalized wavenumber -0.2 0 0.2 0.4

d)

0.2

0.3

0.4

-0.4

Trace number 20

30

Normalized wavenumber -0.2 0 0.2 0.4

0.1

Normalized frequency

Normalized frequency

0.1

10

0.2

0.3

0.4

Figure 18: a) Data with gap inside the original data in Figure 17a. b) Reconstructed data via MWNI. c) and d) are the f-k spectra of (a) and (b), respectively.

Gap sampling Next, we investigate the spatial gap of traces inside a 2D synthetic example (Figure 18a). As it was expected and shown in Figure 18c, the gap causes low amplitude repetition of the spectrum in the vicinity of the original spectrum. Figures 18b and 18d show the reconstructed data using the MWNI method in the t-x and f-k domains, respectively. The results show that the MWNI method is a good candidate for reconstruction of a gap inside the data.

23

3D examples Synthetic linear events In order to examine the performance of the MWNI method for multidimensional data, a 3D cube of synthetic data was created. The original data are composed of three dipping planes with slopes (αx , αy ) = {(2.5, 3.75), (1.25, −1.25), (−3.5, −3.0)}. The data set contains 2400 traces in a spatial grid of 40 × 60 and slightly contaminated with random noise. Figure 19a shows a perspective view of the cube containing the synthetic data. The trajectory lines of the three synthetic events can be seen on the boundaries of the cube. In order to get an in-depth view of the data one can pick slices of data in different locations and project them to the sides of the cube. Figure 19b shows a cube of the original data where the top view is the time slice at 0.65 (s), the front view is the 21st slice in the Y direction, and the side view is 17th slice in the X direction. Figure 19c shows the cube with missing traces on which about 95% of traces have been eliminated randomly. Figure 19d shows the reconstructed data using the MWNI method. Despite the large number of missing traces, the reconstruction was successful. Figures 20a, 20b, and 20c show the f-k panels of the front views in Figures 19b, 19c, and 19d, respectively. Figure 21a shows a cube with missing traces which was created from Figure 19b by first decimating every other slice in the X direction and then eliminating 50% of the remaining traces. Overall, in Figure 21a, about 75% of the original traces were eliminated. Figure 21b shows the reconstructed cube using the MWNI method. It is interesting to note that only every other slice in the X direction is successfully reconstructed. This was expected since the original decimation of slices in X direction results in a 2D Fourier response similar to the one shown in Figure 6d. This means that the replicated artifacts in the X direction are as large as the amplitudes in the original spectrum and cannot be removed by the MWNI method. Since the original data was aliased in both spatial directions we could not use the band-limited MWNI method to eliminate the artifacts. Figures 22a and 22b show the f-k panels of the front views in Figures 22a and 22b, respectively. Synthetic curved events It is interesting to examine the performance of Fourier reconstruction for non-linear (curved) events. Events with curvatures in the t-x domain do not have a simple or sparse representation in the frequency-wavenumber domain. New representation methods like curvelets (Candes et al., 2005) and seislets (Fomel, 2006) were recently introduced to obtain a sparse and local representation of seismic events. This paper, however, is focused on reconstruction methods based on Fourier transforms. Interested readers are referred to Hennenfent and Herrmann (2008); Naghizadeh and Sacchi (2010) for curvelet-based reconstruction and Liu and Fomel (2010) for seislet interpolation. Figure 23a shows a synthetic cube of curved events which are aliased in the X direction and alias-free in the Y direction. Figure 23b shows slice views of the original data at time 0.35 (s) (top view), X slice 17 (side view), and Y slice 21 (front view). Figure 23c shows data after eliminating every other slice of data in the X direction. The sampling operator is equivalent to the one depicted in Figures 6c. Figure 23d shows the reconstructed data using the band-limited MWNI method by applying band-limitations on the normalized

24

b) r

X trace number 15 25 35

X trace number 17

m be

5

m be

r

a)

nu tra ce Y

Y

tra ce

nu

55 45 35 25 15 5

21

Time (s)

Time (s)

0.2 0.5

0.65

0.8

r

d) m tra

ce

nu

m nu ce tra

Y

21

Time (s)

Y

Time (s)

0.65

X trace number 17

be

X trace number 17

be

r

c)

21

0.65

Figure 19: a) 3D synthetic cube of data. b) Slice views from inside the original data. c) Data with 95% randomly missing traces. d) Reconstructed data using MWNI method.

25

Normalized frequency

a)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

b)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

c)

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

Figure 20: a), b), and c) are the f-k spectra of the front views of the data in Figures 19b, 19c, and 19d.

r

b) m tra

ce

nu

m nu ce tra

Y

21

Time (s)

Y

Time (s)

0.65

X trace number 17

be

X trace number 17

be

r

a)

21

0.65

Figure 21: a) Cube of missing traces created from the original data shown in Figure 19b by first eliminating every other X slice and then randomly elimination 50% of the remaining traces. b) Reconstructed data using MWNI.

26

Normalized frequency

a)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

b)

0.1

0.1

0.2

0.2

0.3

0.3

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

Figure 22: a) and b) are the f-k spectra of the front views of the data in Figures 21a and 21b.

wavenumber in the Y direction (ky ∈ [−0.15, 0.15]). The reconstruction was not successful because the sampling artifacts were inside the band-width used by the inversion. Figure 23e shows the cube of data with missing traces which has been sampled in a chessboard pattern. The sampling operator is equivalent to the one in Figure 6e. Figure 23f shows the reconstruction of the data in Figure 23e using the band-limited MWNI method in the bandwidth ky ∈ [−0.15, 0.15]. For the chessboard pattern, the artifacts fall outside the chosen band-width and as a result the band-limited MWNI method can recover the missing traces successfully. Figures 24a, 24b, and 24c show the f-k panel of front views in Figures 23b, 23e, and 23f, respectively.

Interpolation of 3D real data 2D spatial random sampling A real data set from the Gulf of Mexico was selected to analyze the performance of the MWNI method under different sampling scenarios. A data cube of nineteen consecutive shots, each with 91 traces, each trace having 900 time samples, is used for this specific test. Figure 25a shows a perspective view of the cube containing the selected real data with trajectories of the events on the boundaries of the cube. Figure 25b shows slice views of Figure 25a at time 3.5 (s) (top view), X slice 46 (side view), and Y slice 10 (front view). Next 50% of the original traces were eliminated randomly to create a cube with missing traces (Figure 25c). Figure 25d shows the reconstructed data using the MWNI method. Figures 26a, 26b, and 26c show the f-k panels of the front views in Figures 26b, 26c, and 26d, respectively. Because of the presence of curved events there is some differences between the original data and the reconstructed data. Also, notice that the reconstructed data have been muted above the first arrival. By choosing a proper windowing strategy in the spatial and time directions, one can further reduce the differences between the original

27

5

X trace number 15 25 35

b) be r

be r

a)

nu m tra ce Y

Y

tra ce

nu m

55 45 35 25 15 5

X trace number 17

21

Time (s)

Time (s)

0.2 0.5

0.35

0.8

r

d) nu tra

ce

nu ce tra

Y

Time (s)

f) r

X trace number 17

0.35

X trace number 17

m nu ce tra Y

21

Time (s)

Time (s)

Y

tra

ce

nu

m

be

r

e)

21

be

Y

Time (s)

21

0.35

0.35

X trace number 17

m be

X trace number 17

m be

r

c)

21

0.35

Figure 23: a) 3D synthetic cube of data with mild curvature in the Y direction and alias in the X direction . b) Slice view from inside the original data. c) Data after elimination of every other slice in the X direction. d) Reconstruction of (c) using band-limited MWNI. e) Data after chessboard elimination of traces. f) Reconstruction of (e) using band-limited MWNI.

28

Normalized frequency

a)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

b)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

c)

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

Figure 24: a), b), and c) are the f-k spectra of the front views of the data in Figures 23b, 23e, and 23f.

and reconstructed data. 2D spatial regular sampling Figure 27a shows the cube of data from the Gulf of Mexico on which the original traces are decimated by a chessboard pattern. The 2D sampling operator for this example is similar to the one in Figure 6e. Figure 28a shows the f-k spectrum of the front view of Figure 27a. It is clear that regular sampling resulted in an f-k spectrum which is the sum of wrapped (shifted) spectrum with the original spectrum. Also there is a prominent overlap between the original and shifted spectrum. Figure 27b shows the reconstructed data using the MWNI method. Figure 28b shows the f-k spectrum of the front view of Figure 28b. The MWNI method was not able to eliminate the artifacts. Figure 27c shows the reconstructed data using the band-limited MWNI method. The band-limitation was applied in the interval of ky = [−0.2, 0.2] normalized wavenumbers and full-band on kx axis. The band-limitation on the Y direction (common-offset direction) makes sense since one would not expect strong dips as well as sharp changes in this direction. The reconstructed data still contains some leftover artifacts. This could be due to the presence of noise and curvature in the original data. Figure 28c shows the f-k panel of the front view of Figure 28c.

CONCLUSIONS In this paper we examined the effects of random sampling and regular sampling operators on the performance of the MWNI method. The spectrum of sampling operators and their influence on the original spectrum of the data were mathematically derived. The spectrum of the random sampling operators has only one impulse and the rest of its components are

29

b)

ce tra

X trace number 46

10

Y

Y

tra

ce

15 10 5

r

X trace number 40 70

m be

10

nu

nu

m be

r

a)

2.5

Time (s)

Time (s)

1.5

3.5

3.5

4.5

d)

3.5

ce

nu

m

be r

X trace number 46

X trace number 46

10

Y

tra

10

Time (s)

Time (s)

Y

tra

ce

nu

m

be r

c)

3.5

Figure 25: a) original 3D cube of data from Gulf of Mexico. b) Slices view from inside the original data. c) Data with 50% randomly missing traces. d) Reconstructed data using MWNI.

30

Normalized frequency

a)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

b)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

c)

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

Figure 26: a), b), and c) are the f-k spectra of the front views of the data in Figures 25b, 25c, and 25d.

small. Conversely, for regularly sampled data, depending on the decimation factor, several equally sized impulses will be present in the spectral domain. The latter is responsible for spectral mixing in the form of alias. The spectrum of the sampling operator is cyclicly convolved with the spectrum of original data to give the spectrum of sampled data. Provided that the original data had a simple spectrum, its convolution with the spectrum of the random sampling operator will still create a simple signal with small amplitude artifacts. Therefore, randomly sampled signals are recoverable by assuming simplicity in the Fourier domain. The latter can be in the form of sparsity or via the minimum weighted norm constraint studied in this paper. The spectrum of the regular sampling operator contains replications and its cyclic convolution with the spectrum of original data will produce unwanted repetitions of the spectrum of original data. Hence, eliminating the artifacts caused by the regular sampling operator using the MWNI method will be impossible unless the original data were bandlimited. For the multidimensional case (more than one spatial dimensions), the success of the band-limited MWNI method will depend on the type of regular sampling operator and band-limited properties of the original data. For a chessboard pattern decimation, the sampling artifacts in the Fourier domain will have the optimal possible separation from the spectrum of original data. This property helps one to use the band-limitation information in one spatial direction to recover the signal in other, aliased, spatial directions.

ACKNOWLEDGMENTS This research was funded by the generous support of the industrial members of the Signal Analysis and Imaging Group (SAIG) at the University of Alberta and by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada to M.D. Sacchi. We

31

X trace number 46

10

Time (s)

Y

tra

ce

nu

m be

r

a)

3.5

c)

3.5

ce

nu

m

be r

X trace number 46

X trace number 46

10

Y

tra

10

Time (s)

Time (s)

Y

tra

ce

nu

m

be r

b)

3.5

Figure 27: a) Decimation of data in Figure 25b with a chessboard pattern. b) Reconstructed data using MWNI. c) Reconstructed data using the band-limited MWNI method in the interval of ky = [−0.15, 0.15].

32

Normalized frequency

a)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

b)

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

c)

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

Normalized wavenumber Kx -0.4 -0.2 0 0.2 0.4

Figure 28: a), b), and c) are the f-k spectra of the front views of the data in Figures 27a, 27b, and 27c.

also thank inspiring interpolation discussions with Ray Abma, Liu Bin, Mike Perz, Daniel Trad and Juefu Wang.

APPENDIX We call Q the wavenumber response of the sampling operator. It can be proven that Xs = X ∗ Q, where ∗ represents the convolution operator. Writing the convolution operator explicitly, we have: Xs (k) =

N −1 X

X(j)Q(k − j)

j=0

=

−1 N −1 N X X

[

x(n)e

−i2πnj N

m=0

j=0 n=0

=

N −1 N −1 X X

[

x(n)e

j=0 n=0

=

−1 N −1 M X X n=0 m=0

M −1 1 X −i2πh(m)(k−j) N ][ e ] N

x(n)e

−i2πnj N

][

1 N

−i2πh(m)k N

M −1 X

e

−i2πh(m)k N

e

i2πh(m)j N

]

m=0

N −1 1 X −i2πnj i2πh(m)j [ e N e N ]. N j=0

(17)

33

Now due to the orthogonality of Fourier series, the expression inside the bracket in the last line of (17) is nonzero if n = h(m). This means that expression (17) can be reduced to: Xs (k) =

=

M −1 X m=0 M −1 X

x(h(m))e

−i2πh(m)k N

x(h(m))e

−i2πh(m)k N

[

1 N] N k = 0, 1, 2, . . . , N − 1,

(18)

m=0

REFERENCES Abma, R. and N. Kabir, 2006, 3D interpolation of irregular data with a POCS algorithm: Geophysics, 71, E91–E97. Bardan, V., 1987, Trace interpolation in seismic data processing: Geophysical Prospecting, 35, 342–358. Bertero, M. and P. Boccacci, 2003, Super-resolution in computational imaging: Micron, 34, 65273. Beylkin, G., R. Coifman, and V. Rokhlin, 1991, Fast wavelet transforms and numerical algorithms: Comm. on Pure and Appl. Math., 44, 141–183. Candes, E. J., 2006, Modern statistical estimation via oracle inequalities: Acta Numerica, 15, 257–325. Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005, Fast discrete curvelet transforms: Multiscale Modeling and Simulation, 5, 861–899. Darche, G., 1990, Spatial interpolation using a fast parabolic transform: 60th Annual International Meeting, SEG, Expanded Abstarcts, 1647–1650. Donoho, D. L., 2006, Compressed sensing: IEEE Transactions on Information Theory, 52, 1289–1306. Duijndam, A. J. W., M. A. Schonewille, and C. O. H. Hindriks, 1999, Reconstruction of band-limited signals, irregularly sampled along one spatial direction: Geophysics, 64, 524–538. Eldar, Y. C., 2006, Mean-squared error sampling and reconstruction in the presence of noise: IEEE Transactions on Signal Processing, 54, 4619–4633. Feuer, A. and G. C. Goodwin, 2005, Reconstruction of multidimensional bandlimited signals from nonuniform and generalized samples: IEEE Transactions on Signal Processing, 53, 4273–4282. Fomel, S., 2006, Towards the seislet transform: SEG Technical Program Expanded Abstracts, 25, 2847–2851. Gulunay, N., 2003, Seismic trace interpolation in the Fourier transform domain: Geophysics, 68, 355–369. Hennenfent, G., 2008, Sampling and reconstruction of seismic wavefields in the curvelet domain: PhD thesis, The University of British Columbia, Vancouver, BC Canada. Hennenfent, G. and F. J. Herrmann, 2008, Simply denoise: Wavefield reconstruction via jittered undersampling: Geophysics, 73, V19–V28. Liu, B., 2004, Multi-dimensional reconstruction of seismic data: PhD thesis, University of Alberta. Liu, B. and M. D. Sacchi, 2004, Minimum weighted norm interpolation of seismic records: Geophysics, 69, 1560–1568.

34

Liu, Y. and S. Fomel, 2010, Oc-seislet: seislet transform construction with differential offset continuation: Geophysics, Accepted for special issue of sampling. Naghizadeh, M. and M. D. Sacchi, 2007, Multistep autoregressive reconstruction of seismic records: Geophysics, 72, V111–V118. ——– 2010, Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data: Geophysics, Accepted for special issue of sampling. Oppenheim, A. V. and R. W. Schafer, eds. 1974, Digital signal processing: Prentice-Hall Inc. Park, S. C., M. K. Park, and M. G. Kang, 2003, Super-resolution image reconstruction: a technical overview: IEEE Signal Processing Magazine, 20, 21–32. Porsani, M., 1999, Seismic trace interpolation using half-step prediction filters: Geophysics, 64, 1461–1467. Ronen, J., 1987, Wave-equation trace interpolation: Geophysics, 52, 973–984. Ronen, S., D. Nichols, R. Bale, and R. Ferber, 1995, Dealiasing dmo: Good-pass, bad-pass, and unconstrained: SEG Technical Program Expanded Abstracts, 14, 743–746. Ronen, S., V. Sorin, and R. Bale, 1991, Spatial dealiasing of 3-d seismic reflection data: Geophysical Journal International, 503–511. Sacchi, M. D. and T. J. Ulrych, 1996, Estimation of the Discrete Fourier Transform, a linear inversion approach: Geophysics, 61, 1128–1136. Sacchi, M. D., T. J. Ulrych, and C. J. Walker, 1998, Interpolation and extrapolation using a high-resolution discrete fourier transform: IEEE Transaction on Signal Processing, 46, 31–38. Schonewille, M. A., R. Romijn, A. J. W. Duijndam, and L. Ongkiehong, 2003, A general reconstruction scheme for dominant azimuth 3D seismic data: Geophysics, 68, 2092–2105. Spitz, S., 1991, Seismic trace interpolation in the F-X domain: Geophysics, 56, 785–794. Strohmer, T., 1997, Computationally attractive reconstruction of bandlimited images from irregular samples: IEEE Transactions on Image Processing, 6, 540–548. Trad, D., 2008, Five dimensional seismic data interpolation: 78th Annual International Meeting, SEG, Expanded Abstarcts, 978–981. ——– 2009, Five-dimensional interpolation: Recovering from acquisition constraints: Geophysics, 74, V123–V132. Tsaig, Y. and D. L. Donoho, 2006, Extensions of compressed sensing: Signal Processing, 86, 549–571. Unser, M., 2000, Sampling - 50 years after Shannon: Proceedings of the IEEE, 88, 569–587. Vandewalle, P., L. Sbaiz, J. Vandewalle, and M. Vetterli, 2007, Super-resolution from unregistered and totally aliased signals using subspace methods: IEEE Transactions on Signal Processing, 55, 3687–3703. Xu, S., Y. Zhang, D. Pham, and G. Lambare, 2005, Antileakage Fourier transform for seismic data regularization: Geophysics, 70, V87–V95. Zwartjes, P. and A. Gisolf, 2006, Fourier reconstruction of marine-streamer data in four spatial coordinates: Geophysics, 71, V171–V186. Zwartjes, P. and M. D. Sacchi, 2007, Fourier reconstruction of nonuniformly sampled, aliased seismic data: Geophysics, 72, V21–V32.