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.
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.
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,
Approximation of covariances
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.
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.
To estimate , 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.
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:
Computing the NDFT takes 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).
For each iteration, the CG scheme takes 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.
APPLICATIONS AND RESULTS
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.