Seismic signals are often irregularly sampled along spatial coordinates, leading to suboptimal processing and imaging results. Least-squares estimation of Fourier components is used for the reconstruction of band-limited seismic signals that are irregularly sampled along two spatial coordinates. A simple and efficient diagonal weighting scheme, based on the areas surrounding the spatial samples, takes the properties of the noise (signal outside the bandwidth) into account in an approximate sense. Diagonal stabilization based on the energies of the signal and the noise ensures robust estimation.

Reconstruction by temporal frequency component allows the specification of varying bandwidth in two dimensions, depending on the minimum apparent velocity. This parameterization improves the reconstruction capability for lower temporal frequencies. The shape of the spatial aperture affects the method of sampling the Fourier domain. Taking into account this property, a larger bandwidth can be recovered.

The properties of the least-squares estimator allow a very efficient implementation which, when using a conjugate gradient algorithm, requires a modest number of 2-D fast Fourier transforms per temporal frequency. The method shows signicant improvement over the conventionally used binning and stacking method on both synthetic and real data. The method can be applied to any subset of seismic data with two varying spatial coordinates.


Seismic data are typically irregularly sampled along spatial coordinates. Among other reasons, this is caused by the presence of obstacles, inaccessible areas, feathering of streamers, dead channels and economics. Most multitrace data processing algorithms cannot adequately handle irregular sampling. For the situation with two spatial coordinates, as considered in this paper, a popular technique is binning, where seismic traces are simply assigned to the nearest bin centers with no amplitude and phase corrections. Equivalent to introducing a significant positioning error, this method is clearly not optimal. A more accurate procedure is needed, in which the signal is reconstructed from the irregularly spaced samples, either in the spatial domain or in a domain associated with a suitable transformation (e.g, a Fourier transformation).

Reconstruction of data irregularly sampled along more than one coordinate is still an area to be investigated for seismic data. Most methods described in the literature make some assumptions which are not applicable to seismic data, such as infinite sampling and high sampling rates in regions where higher wavenumber components have contributing energy. Strohmer (1993) proposed a reconstruction algorithm for two dimensions that only assumes the data are band limited, which makes the algorithm able as a base for reconstruction of seismic signals.

In this paper, we extend the 1-D reconstruction method proposed by Duijndam et al. (1999) to two spatial dimensions. The method estimates Fourier coefficients using a least-squares approach. It does not require any geological or geophysical assumptions other than that the signal is band limited in two dimensions. A similar algorithm has been proposed by Rauth and Strohmer (1998) for potential field data. Here, we discuss the special aspects related to seismic data. The computational steps of the algorithm, which can be carried out very efficiently, are discussed.

The method can be applied to any subset of seismic data with two varying spatial coordinates. Specific applications are the construction of a pseudo zero-offset data set from NMO-corrected data (alternative to binning and stacking), the reconstruction of prestack data along in-lines or cross-lines using midpoint and offset coordinates, the reconstruction of 2-D shot or receiver gathers (if sampled densely enough), and the reconstruction of 3-D vertical seismic profiling (VSP) data sets.


We define the continuous 2-D forward spatial Fourier transformation as
where x and y are the spatial variables, kx and ky are the corresponding wavenumbers, and ω is the temporal frequency. The inverse transformation is given by
For band-limited data, regularly sampled along x and y, the well-known discrete Fourier transform (DFT) related to equation (1) is used:
with ΔS = Δx × Δy. The spatial sampling intervals Δx and Δy should be chosen small enough to avoid aliasing in the spatial Fourier domain.
In the presence of irregular sampling, a straightforward approach to obtain the data in the Fourier domain is to use the Riemann sum, where the integral in equation (1) is replaced by a sum corresponding to the actual sample locations graphic:
where ΔSn is the area corresponding to the sample (xn, yn). In contrast to the 1-D problem the assignment of this area is not trivial. Triangulation algorithms can be used to determine it. We use a Delauney triangulation (Ahuja and Schacter, 1983). After partitioning the spatial sampling area in triangles, each sample is weighted with a value ΔSn equal to the sum of one-third part of the area of each adjacent triangle (Figure 1). The samples defining the edges of the sampling area are assigned smaller weights than deserved, which has the effect of a small taper. We will refer to the transform of equation (4) as the 2-D nonuniform discrete Fourier transform (NDFT). The 2-D NDFT does not reproduce the Fourier spectrum exactly, as will be demonstrated later.

Forward model

Extending the 1-D reconstruction algorithm, as described by Duijndam et al. (1999), the next approach is proposed in order to improve upon the NDFT.

For practical purposes, let us consider a regular sampling of the Fourier domain where the data in the Fourier domain is restricted to a certain limited area (band limited in the two dimensions). The discrete inverse Fourier transform for any spatial location (x, y) is given by
where the precise specification of (kx,m, ky,m) and the elementary area ΔSF in the Fourier domain depend on the specific parameterization chosen (see below).

The discrete inverse Fourier transform is exact, except for a possible incomplete coverage of the so-called “region of support” (ROS), i.e., the region in the Fourier domain where the signal is located,

and apart from spatial aliasing, which depends on the chosen sampling of the Fourier domain. The combination of N versions of equation (5) for the irregularly spaced sample locations graphic can be written in vector notation as
In practice the data will not be exactly band limited, and some spatial frequency components exist beyond the bandwidth that we specify the data to have. These components constitute the errors or “noise" in the forward model and should be incorporated in a noise term n:

MAP estimator

We now recognize a standard linear inverse problem similar to the 1-D irregular sampling problem (Duijndam et al., 1999), in which the unknown vector graphic, containing the data in the regularly sampled Fourier domain, is to be estimated from the data vector y containing the irregularly sampled data in the spatial domain. Any desired parameter estimation technique can be used. Within the Bayesian formalism (Bard, 1974; Tarantola, 1987; Duijndam, 1988) under Gaussian noise assumptions and zero a priori model, the maximum of the a posteriori probability density function (MAP) estimator is given by
where AH denotes the complex conjugate transpose of A, Cn is the covariance matrix of the noise, and Cgraphic is the a priori covariance matrix of the signal.

Approximation of covariances

An accurate representation of the noise in Cn is expensive and typically not realistic in the sense that detailed information on the noise is not available. Duijndam et al. (1999) have shown that a practical approximation in case of 1-D irregular sampling is obtained by taking
with W a diagonal matrix containing the areas surrounding the samples,
Equation (12) can be shown to hold for 2-D irregular sampling as well. However, determination of the weights is not trivial in case of 2-D irregular sampling; similar to the NDFT of equation (4), we used a Delauney triangulation algorithm. Due to the use of the matrix W, closely spaced samples have less weight than widely spaced samples. Duijndam et al. (1999) already gave the argument that for band-limited noise, closely spaced samples contain almost the same signal and noise values; hence, they do not represent independent information and should instead approximately be regarded as one piece of information or sample, which is achieved by using the weights.
In general, we have no reason to assume correlations between the spatial frequencies and Cgraphic will be diagonal. By varying the diagonal, it is possible to specify a preference for certain spatial frequencies to others (e.g., when at events are most likely to occur). Usually, however, Cgraphic will take the form
where σ2graphic is the expectation of the signal. By substituting the approximated covariance matrix of the noise, equation (12), and the covariance matrix of the signal, equation (14), in the estimator of equation (11), the estimator becomes
with k2 = c22graphic. If one is prepared to specify the ratio F between the expected energies of the noise and the signal in the input data, then it can be shown ( Appendix A) that the stabilization constant can be approximated by
where M is the total number of Fourier coefficients. Otherwise, a rule of thumb such as taking 1% of the value of the diagonal or a trial and error procedures can be used.
When we evaluate the term AHWy on the right side of equation (15), we obtain
and we observe that, apart from the constant factor ΔSF/4π2, this first term represents the Riemann sum (NDFT) of equation (4) and is a first approximation to the answer. The operator graphic in equation (15) “deconvolves” or improves the output of the NDFT. Without weights the first term would be the poorer approximation AHy, and the corresponding correction process described by the operator graphic would show slower convergence when using a conjugate-gradient scheme.
The estimator of equation (15) for the general situation of irregular sampling reduces to a standard result for the special case of regular, orthogonal sampling, with rectangular region of support. This means that the sample positions are given by (nxΔx, nyΔy). With Nx and Ny sample points in x and y, we have a total number of sample points N = NxNy. It can be shown that for the choices Δkx = 2π/NxΔx and Δky = 2 π/NyΔy the estimator becomes
Setting the stabilization term to zero, substitution of the definitions of graphic, A, and y in the estimator of equation (18), and using ΔSF = ΔkxΔky indeed yields the standard DFT [equation (3)]:
Apart from the weighting matrix W and a proper stabilization, the methodology obtained so far is comparable with those derived by others. We already mentioned the work done by Duijndam et al. (1999), who investigated 1-D reconstruction. In papers by Thorson and Claerbout (1985) and Kabir and Verschuur (1995), the equations for the least-squares Radon transform, used for rconstruction, were set up in a similar way. Rauth and Strohmer (1998) proposed a similar method for reconstructing 2-D irregularly sampled potential field data. However, in this paper, we are dealing with a problem concerning 2-D irregularities within a 3-D seismic experiment. Regularization has to be carried out for each time or temporal frequency component, the parameterization is less trivial compared to 1-D irregular sampling and, since we are dealing with large systems to be solved, emphasis should be put on efficiency as well.


The choice of parameterization is crucial, as in most if not all inverse problems. A bad choice of parameterization can destabilize the system to be solved and introduce aliasing. Compared to the 1-D case we have more freedom in sampling the Fourier domain.

Petersen and Middleton (1962) discussed multidimensional regular sampling. Nonorthogonal regular sampling in the Fourier domain can be described by two base vectors, each with two components. The position of the mth Fourier domain sample on an arbitrary sampling grid determined by the real base-vectors v1 and v2 is then described by
The vectors v1 and v2 determine how the inverse Fourier transformation of equation (5) is discretized. Due to this discretization, ΔSF is equal to the area surrounding each Fourier domain sample. As this area is spanned by v1 and v2, ΔSF is determined by the length of the crossproduct of v1 and v2,
and the forward model of equation (5) can be rewritten as,
By discretizing the inverse Fourier transformation, the signal will become periodic in the spatial domain. The spatial periodicity is described by two base vectors u1 and u2 related to v1 and v2 by
The orientation and length of the vectors v1 and v2 should be chosen such that overlap, or aliasing, will not occur in the spatial domain. For a rectangular ROS in the spatial domain with lengths of X and Y along, respectively, the x- and y-axis, we may allow a periodicity described by the vectors u1 and u2 of
Through equation (24), this corresponds to
In case the sampling area has an elliptical-shaped region of support with axes parallel to the x- and y-axes, the optimal periodicity is determined by a closest packing of ellipses (Figure 2). The base vectors describing the repetition in the spatial domain are then defined as
which leads to the following base vectors describing the sampling in the wave-number domain:
These definitions for the base vectors, describing the sampling positions, lead to a less dense sampling in the wavenumber domain compared to orthogonal sampling as shown in Figure 2. A closest packing of repeated regions of support in the spatial domain leads to the least amount of parameters covering the 2-D ROS in the Fourier domain, resulting in a better determination of the estimator of equation (15).

Apart from the sampling grid, we haveto determine the size and the shape of the ROS in the Fourier domain. For many cases in seismic analysis, it will be reasonable to assume a minimum apparent velocity cmin that will not vary with azimuth.

Through Fourier theory, it follows that the ROS is then limited by this minimum velocity and temporal frequency through the condition
The volume, in which the relevant data resides, has a cone shape, determined by the minimum velocity (see Figure 3). For low temporal frequencies, fewer parameters have to be estimated, which results in a better determined system of equations. Because the stabilization factor depends on the number of parameters M [equation (16)], it will increase along with the temporal frequency. Hence, when the method is applied to severely irregularly and sparsely sampled signals, data are suppressed within gaps for high temporal frequencies only, while for low temporal frequencies, hardly any stabilization is needed, resulting in correct reconstruction. When reconstruction is done in the time domain, the number of parameters to be estimated must be chosen corresponding to the highest frequency. This leads to a less stable and efficient method. Therefore, the temporal frequency domain approach is preferred. In contradiction to Rauth and Strohmer (1998), who proposed a similar method to be used on potential field data, the parametrization we propose is justified by the physical properties of the seismic experiment.


To estimate graphic, the signal in the Fourier domain, from its spatial samples y by means of equation (15), two major computational efforts have to be made for each frequency component. First, the computation of the NDFT, AHWy, is required. For this purpose nonuniform fast Fourier transform (NFFT) schemes have been developed which only take a number of operations proportional to the number of operations for the fast Fourier transformation (Duijndam and Schonewille, 1999). The second computational effort is the inversion of (AHWA + k2I).

Because we are dealing with a number of parameters on the same order as the number of input samples, which can be on the order of tens of thousands, standard inversion will be very expensive. However, (AHWA + k2I) turns out to have special properties, which can be used in deriving a very efficient inversion scheme.

Let graphic be ordered in the following way:
with M1 × M2 being the total number of parameters to be estimated. The matrix element [AHWA]pq is determined by the inner product of the pth row of AH, associated with (kx, ky)p, and the qth column of A, associated with (kx, ky)q. Through equation (20), let (kx, ky)p be defined by some scalars p1 and p2, and (kx, ky)q by some scalars q1 and q2:
Then, by combining equations (9), (13), and (32), [AHWA]pq can be written as
Through the ordering as defined by equations (30) and (31), we can draw the following conclusions. (1) As long as the parameters p1 and q1 are running, p2 and q2 are constant. Therefore, the matrix AHWA is made of blocks with constant p2 and q2:
with Bm = BHm and
and with m = p2q2. (2) Within each block, the factor (p1q1) is constant on the diagonals. Therefore, the blocks are Toeplitz. (3) The same holds for the factor m = (p2q2) with respect to distinctive blocks. Therefore, AHWA is block Toeplitz as well. (4) AHWA is Hermitian. The operator (AHWA + k2I) has a so-called block-Toeplitz-Toeplitz-block (BTTB) structure and is independent of the temporal frequency, apart from k2I.
A Toeplitz matrix T can be augmented to a circulant matrix C as follows:
Then, within a matrix operation on an arbitrary vector x, T and C are related through
in which y represents residual elements we are not directly interested in. With the definition of C [equation (36)], the left side of equation (37) resembles a discrete convolution of x with the first row of C, which can also be computed using fast Fourier transforms (FFTs). In a similar way, relations can be found to augment a BTTB matrix to a block-circulant-circulant-block (BCCB) matrix, representing a 2-D convolution operator. The 2-D convolution that is defined in this way can be carried out by using a 2-D FFT. After augmenting (AHWA + k2I) to a BCCB matrix, all distinctive elements are collected in a matrix H, and we have
which is the 2-D equivalent of the first row of C of equation (37). Matrix-vector multiplication (AHWA + k2I)x, with arbitrary x, is replaced by an elementwise multiplication with the 2-D FFTs of H and x, followed by an inverse 2-D FFT (Strang, 1986). These aspects can be advantageously incorporated within iterative conjugate-gradient schemes.
The estimator of equation (15) can be written as follows:
where L is the bounded matrix (AHWA + k2I) and f = AHWy. The matrix L is Hermitian: LH = L.
Through standard derivations of conjugate gradient (CG) schemes (Kleinmann and Van den Berg, 1991), the following iterative scheme is found, which solves the linear set of equations [equation (40)] by minimizing the residual ‖Lgraphicf‖. Starting at n = 1,
The only computationally intensive element in each iteration is the computation of Lwn which, as shown above, can be carried out by using a 2-D FFT.

Recalling that the BTTB operator is data independent and needs to be computed only once, the total number of operations for each frequency component is derived from the following contributions:

  1. Computing the NDFT takes graphic operations, with N the number of input samples and M the total number of parameters to estimate. Typical values for q and f are q = 7 and f = 2 and (Duijndam and Schonewille, 1999).

  2. For each iteration, the CG scheme takes graphic operations, with niter the number of required iterations.

In total, the estimation of the wavenumber parameters per temporal frequency takes the same amount of operations as a modest number of 2-D FFTs does.


The 2-D reconstruction is illustrated on synthetic data with a random sampling geometry. The sampling geometry and the model are shown in Figure 4. The model consists of six plane-dipping events. The sixth event is undersampled with regard to the spatial bandwidth. The ROS in the Fourier domain in this example is chosen such that the sixth event is not included and therefore serves as noise. Its energy is 5% of the total energy. Regarding this event, the factor F in equation (16) is set at 0.05.

The data has a temporal frequency content of 125 Hz and is spatially 20% over-sampled with respect to the five properly sampled events. The input data is tapered towards the edges to avoid edge effects. The signal is “estimated” in the Fourier domain with three methods: (1) conventional binning and stacking, (2) the NDFT, and (3) the 2-D reconstruction technique. All figures are plotted in absolute values on the same amplitude scale.

Figure 5 shows the wavenumber spectra for a frequency slice (f = 58.5 Hz). Figure 5a shows the reference data. Figure 5b shows the result after applying conventional binning and stacking. The result is distorted throughout the spectrum. Furthermore, the spectrum shows attenuation towards the edges of the spectrum. The NDFT (Figure 5c), computed within the bandwidth we specify the data to have, shows the same amount of “background noise”; however, the events do not show attenuation towards the edges, resulting in a better signal-to-noise ratio. Figure 5d shows the result after reconstruction (or improvement of the NDFT). The observed distortion is due only to the undersampled event, which is not within the specified bandwidth. Clearly the reconstruction technique has the best definition, the least distortion, and shows no attenuation towards the edges of the spectrum.

Figure 6 shows the amplitude spectra as a function of kx and ω for y = 0. In Figure 6c (result after applying the NDFT) and Figure 6d (result after reconstruction), clearly the bandwidth we specified the data to have is visible according to a defined minimum apparent velocity. Conclusions similar to those for Figure 5 can be drawn.

Figure 7 shows the reconstructed, or “regularized,” data in the spatial domain as a function of x and t for y = 0. Binning and stacking the data (Figure 7b) renders no data in the gaps, displayed as zero traces. The NDFT (Figure 7c) partly fills the gaps but still shows ringing of the signal. The reconstruction (Figure 7d) shows excellent reconstructed, continuous events. The event specified to be outside the bandwidth is suppressed to a limited extent. The higher the number of input samples, the better the anti-alias protection capability of the reconstruction.

3-D VSP walkaway experiment

To illustrate the 2-D reconstruction on field data, a 3-D VSP experiment is performed (courtesy of Norsk Hydro). This marine 3-D walkaway VSP experiment consists of a downhole receiver line and 9930 source positions.

The source locations are irregularly distributed (see Figure 8), and reconstruction is carried out over the source coordinates. The acquisition geometry shows two large gaps, one in the center of the survey (caused by the maneuvre limits of the seismic vessel) and one in the lower right corner (due to a platform, which formed an obstacle in the survey).

The orthogonal sampling used in the Fourier domain corresponds with a square area of 4 × 4 km in the spatial domain. Since this area is not fully covered by the acquisition, no useful results can be expected in the corners. We will illustrate the method on one receiver gather.

The data has a temporal frequency content of 100 Hz. The average sampling rate in both the x- and y-directions is approximately 40 m. With regard to these sampling rates, aliasing is observed above a temporal frequency of 30 Hz. The data is not tapered towards the edges, and no further preprocessing steps have been applied.

Figure 9 shows the amplitude spectra after binning and stacking the data (left figure) and after applying the reconstruction method (right figure) at a frequency of 7 Hz. It is clearly visible that, except for the noise, the bandwidth of the signal is bounded by the presence of a minimum apparent velocity. We defined the ROS to cover a larger area (shear-wave energy will be present as well). Binning and stacking (left figure) has more distortion and shows attenuation within the bandwidth in comparison with the reconstruction technique.

Figure 10 shows the amplitude spectra in the kx-ω domain along the line y = 0 up to 12 Hz, clearly visualizing the cone-shaped ROS we have chosen and the presence of a minimum acoustic velocity. Again similar observations can be made.

Figure 11 shows the results of the binning and stacking method (left figure) and the reconstruction method (right figure) in the space-time domain. The upper part of the cubes are cut along the line AB (Figure 8) and t = 2 s. The binning and stacking method renders zero data in the gaps and a poor definition of the events, whereas the reconstruction method shows good continuation of the events. Even the very large gaps are partly filled.


We have derived an algorithm for the reconstruction of band-limited signals that are irregularly sampled along two spatial coordinates. It is generally applicable since it is based on Fourier theory and discrete inverse theory only and requires no further geophysical assumptions. Moreover, it can be used for arbitrary sampling geometries. Reconstruction per temporal frequency component allows a varying spatial bandwidth, depending on a minimum apparent velocity to be specified, resulting in greater reconstruction capability than reconstruction per time slice.

A diagonal weighting scheme based on the area surrounding the input samples is an approximate and very efficient way of taking the noise properties into account. With a slight average oversampling, gaps up to three times the Nyquist interval in both spatial directions can be reconstructed. For larger gaps, stabilization is required. The stabilization constant can be based on the energies of the signal and the noise (signal outside the bandwidth).

The algorithm can be applied to any subset of (preprocessed) traces where two spatial coordinate vary. Specific applications are the replacement of conventional binning and stacking and the regularization of 3-D VSP data. The improvement over the bin stack consists of both an improved anti-alias protection and the reduction of signal distortion within the bandwidth.

A very fast implementation of the algorithm is possible, which takes, apart from FFTs to the temporal Fourier domain and back, a modest number times the 2-D FFT corresponding to the number of spatial output samples.

It is important to note that through the mild assumptions chosen for this algorithm, it is not aimed at “interpolation beyond aliasing.” Severely undersampled events cannot be reconstructed, and an uncertainty analysis will show this to the user. Amplitude or static time shift perturbations of the data set will broaden the spatial bandwidth and can best be removed before application of the reconstruction algorithm.


We thank Norsk Hydro for the provision of their 3-D VSP data and the permission to publish the results based on their data.



The stabilization factor k2 = c22graphic depends on the ratio of second-order moments (or expectations of energies) of the noise and the signal. When the energy of the noise is believed to be a factor F times the energy of the signal within the bandwidth considered, we can write in the continuous domain, using Parseval's relation,
In order to obtain practical results, we have to approximate equation (A-1) by discretizing both sides. The noise is discretized irregularly, the right side is discretized in a regular fashion. Furthermore, we take the expectations of both sides:
with ni = N(xi) and graphic. Using
and graphic we obtain:

First Page Preview

First page PDF preview