Fourier methods of seismic data regularisation

Irregular Fourier Transform

The regular unitary discrete Fourier transform can be expressed as follows

$$F(k) = \frac{1}{\sqrt N} \displaystyle\sum_{l=0}^{N-1} f(x_l) \ e^{-2\pi ikx_l}, \quad k = 0, \pm 1, \pm 2,..., \pm N$$ (3)

where $F(k)$ denotes the size of the $k^{th}$ Fourier coefficient, $x_l$ is the position of the input data and $f(x_l)$ denotes the value of the input data series corresponding to this position. To extend this expression to handle an irregular $x_l$ series we must define a new transform, and care must be taken with the transform weights to ensure the forward and inverse transforms are unitary; this is simply done by dividing by $\sqrt N$ in the regular case.

The weights must be related to the relative positioning of $x_l$ on the input axis, with some form of normalisation. This can be done by letting

$$F(k) = \frac{1}{\sqrt N} \displaystyle\sum_{l=0}^{N-1} f(x_l) \ \frac{\Delta x_l}{\Delta x_{av}} \ e^{-2\pi ikx_l},$$ (4)

$$= \frac{1}{\sqrt N} \displaystyle\sum_{l=0}^{N-1} \frac{1}{\Delta x_{av}} \; \displaystyle\sum_{l=0}^{N-1} f(x_l)\Delta x_l \ e^{-2\pi ikx_l},$$ (5)

$$= \frac{\sqrt N}{\Delta X} \; \displaystyle\sum_{l=0}^{N-1} f(x_l)\Delta x_l \ e^{-2\pi ikx_l},$$ (6)

$\Delta X$ is the range of the input axis, $x_N - x_1$ . The weighting here has taken into account the data range, size, and relative irregularity. Upon multiple tests it is easy to show that this transform is indeed unitary, and corretly preserves amplitudes after hundreds of domain transformations. The $\Delta x_l$ values are calculated using a symmetric differencing star for stability, such that

$$\Delta x_l \ = \frac{1}{2} (x_{l+1} - x_{l-1})$$ (7)

a weighting scheme for applying the transform in multiple dimensions is more difficult to estimate, but the problem can be broken down and each dimension considered separately for this application.

However designing an appropriate weight does not alleviate spectral leakage. One first method used was setting up the problem as a least squares inverse, using the Cauchy norm of $F(k)$ as a regularisation term within the cost function (Saachi and Ulrych, 1995), however a key problem with this approach was that the data is not sufficiently honoured.

Fourier methods of seismic data regularisation