The aim of exploration seismology is to obtain an image of the subsurface by probing it with seismic waves at various locations. These waves are generally generated by using airguns in marine, and vibroseis or dynamite in land. They propagate downwards through the subsurface, and are reflected at interfaces between geological layers or refracted within layers. Parts of these waves subsequently propagate upwards to the surface, where they are detected and recorded.
In exploration seismology, although the time coordinate is regularly sampled, spatial coordinates are typically irregularly sampled due to the presence of obstacles in land and strong currents in marine. But even for example receivers placed within a marine survey cable or streamer may not be always equidistant. Hence, the inline sampling can be quite irregular.
The regularization of seismic data is very important especially in time-lapse survey matching, multiple suppression and imaging. If the irregular nature of the sampling grid is ignored or handled poorly, notable errors are introduced the severity of which may be further amplified at later stages of the seismic processing chain. The problem of signal reconstruction from uniformly spaced data has been investigated in depth. The Whittaker-Kotel'nikov-Shannon sampling theorem states that any signal f(x) can be reconstructed from its uniformly spaced samples if the sampling interval is less than half the period of the highest frequency component in that signal. Thus if f(x) is bandlimited to the wavenumber s/2, which is known as the Nyquist wavenumber, then the sampling theorem provides the following formula to interpolate any function value from uniformly spaced values f(m/s):
                              f          ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                                    f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢            sin            ⁢                                                  ⁢                          c              ⁡                              (                                                      σ                    ⁢                                                                                  ⁢                    x                                    -                  m                                )                                                                        [        1        ]            where sinc(x)=sin(πx)/πx. Thus, when the sampling rate is sufficient and there is no aliasing, the sampling theorem provides a way to reconstruct the signal “exactly” from its uniformly spaced samples. To satisfy requirements of the sampling theorem, the signal should be sampled at a rate greater than twice the Nyquist rate, i.e., s. The seismic signal may be given by:
                                          f            ⁡                          (              x              )                                ≈                                    f              L                        ⁡                          (              x              )                                      =                  σ          ⁢                                    ∑                              m                =                0                                            L                -                1                                      ⁢                          Δ              ⁢                                                          ⁢                              x                m                            ⁢                              f                ⁡                                  (                                      x                    m                                    )                                            ⁢              sin              ⁢                                                          ⁢                              c                ⁡                                  (                                      σ                    ⁡                                          (                                              x                        -                                                  x                          m                                                                    )                                                        )                                                                                        [        2        ]            where Δxm is the Jacobian weight, i.e., Δxm=xm+1−xm and f(xm) the value of the seismic data at irregular offset xm. It is important to note that, when Δxm=1/σ the sinc interpolator of eq. [2] is exact since
                              f          ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                                    f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢            sin            ⁢                                                  ⁢                          c              ⁡                              (                                  σ                  ⁡                                      (                                          x                      -                      τ                                        )                                                  )                                                                        [        3        ]            by Whittaker-Kotel'nikov-Shannon theorem. On the other hand, when Δxm is not equal 1/σ, the sinc interpolator provides only a crude approximation to the continuous signal.
For such cases, it has been found that a better approach is to invert equation [1] for the desired uniformly spaced signal values f(m/σ). This inversion when written in matrix notation yields
                    h        ≡                                                                                        f                  ⁡                                      (                                          x                      1                                        )                                                                                                                        f                  ⁡                                      (                                          x                      2                                        )                                                                                                      M                                                                                      f                  ⁡                                      (                                          x                      L                                        )                                                                                                ≈                                                                                                              s                    11                                                                                        s                    11                                                                    Λ                                                                      s                                          1                      ⁢                                                                                          ⁢                      L                                                                                                                                        s                    21                                                                                        s                    11                                                                    Λ                                                                      s                                          2                      ⁢                                                                                          ⁢                      L                                                                                                                    M                                                  M                                                  O                                                  M                                                                                                  s                    11                                                                                        s                    11                                                                    Λ                                                                      s                    LL                                                                                            ·                                                                                                      f                    ⁡                                          (                      0                      )                                                                                                                                        f                    ⁡                                          (                                              1                        /                        σ                                            )                                                                                                                                                                                                                                                            f                    ⁡                                          (                                                                        (                                                      L                            -                            1                                                    )                                                /                        σ                                            )                                                                                                                      ≡        Sg                            [        4        ]            where σ/2 is the bandwidth of the signal f(x) and S is the sinc matrix with entries sij=sinc(σ(xi−j/σ). If the matrix S is well conditioned than the seismic data at regular offsets can be computed by standard matrix inversion:g=S−1h  [5]Otherwise, a least squares minimum norm inversion can be used:g=(STW1S+W2)−1STW1h  [6]where W1 is usually chosen as a diagonal matrix whose mth diagonal entry is the Jacobian weight Δxm=xm+1−xm and W2 is usually chosen as a small multiple of identity matrix, i.e., W2=ε2I.
This necessitates the use of robust and efficient techniques for seismic data regularization.
An interpolator in accordance with equation [6] has been described in: Yen J. L., 1956, On nonuniform sampling of bandwidth-limited signals, IRE Trans. Circuit Theory, CT-3, 251-257 (1956). It is therefore sometimes referred to as Yen's Interpolator of Type 1.
Many interpolators used in seismics are variations of the Yen's interpolator for example those described in: Duijndam, A. J. W., Schonewille, M. A., Hindriks, C. O. H., Irregular and sparse sampling in exploration seismology (chapter of L. Zhang: Nonuniform sampling: theory and practice, Kluwer Academic/Plenum Publishers, New York, USA (2001)). The regularization in Duijndam is formulated as a spectral domain problem. There it is proposes to estimate the spectrum of the signal by using non-uniform Fourier transform of the irregular samples, and then reconstruct the regular samples by using inverse discrete Fourier transformation. It can be shown that, this regularization technique is exactly equivalent to Yen's Interpolator of Type 1.
A variant of this method is described in: Zwartjes, P. M. and M. D. Sacchi, 2004, Fourier reconstruction of non-uniformly sampled, aliased data: 74th Ann. Internat. Mtg.: Soc. of Expl. Geophys., 1997-2000. However in Zwartjes there is proposed least squares inversion of the Fourier transformation instead of using an inverse discrete Fourier transform. To this purposes a cost function is defined, which also involves a non-quadratic penalty term to obtain a parsimonious model.
Another example is: Hale, I. D., 1980, Resampling irregularly sampled data, Stanford Exploration Project, SEP-25, 39-58. Hale's method is based on the more general version of the Yen's Interpolator, where a space limited signal assumption is not used. In that case, the uniform samples f(m/σ) can still be computed by solving a matrix equation similar to [6]. What Hale suggests is to replace the entries in the inverse matrix by their locally computed approximations.
The interpolators based on Yen's 1st theorem usually provide satisfactory results on non-aliased signals with little high-wavenumber content. However their performance degrades significantly when the interpolated signal has a substantial amount of high wavenumber spectral content. Another shortcoming of the interpolators based on Yen's 1st theorem is that in order to solve [4], at least as many irregular sampling positions as regular sampling positions are required. Hence, if some seismic traces are dropped out, traces which reside at further locations must be used to solve the system of equations given by [4]. Usually this degrades the accuracy of the interpolated sample values.
Further, although Yen's 1st interpolator is exact for infinite length signals, it is an approximation when only a finite extent of the signal is available for interpolation.
Given the problems of the existing interpolators it remains an object to find improved interpolates capable of interpolating data received by receivers at irregular locations to regulars sampling locations.