Multidimensional seismic data processing

Multichannel filters may be applied to each individual record or to the final seismic profile. This may be done to separate different types of waves and to improve the signal-to-noise ratio. There are two well-known methods of designing velocity filters for seismic data processing applications.[2]

The two-dimensional Fourier transform is defined as:

${\displaystyle F({\underline {k}},\omega )=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f({\underline {x}},t)e^{-j(\omega t-{\underline {k}}{\underline {x}})}d{\underline {x}}dt}$

where ${\displaystyle {\underline {k}}}$ is the spatial frequency (also known as wavenumber) and ${\displaystyle \omega }$ is the temporal frequency. The two-dimensional equivalent of the frequency domain is also referred to as the ${\displaystyle {\underline {k}}-\omega }$ domain. There are various techniques to design two-dimensional filters based on the Fourier transform, such as the minimax design method and design by transformation. One disadvantage of Fourier transform design is its global nature; it may filter out some desired components as well.

The τ-p transform is a special case of the Radon transform, and is simpler to apply than the Fourier transform. It allows one to study different wave modes as a function of their slowness values, ${\displaystyle p}$.[3] Application of this transform involves summing (stacking) all traces in a record along a slope (slant), which results in a single trace (called the p value, slowness or the ray parameter). It transforms the input data from the space-time domain to intercept time-slowness domain.

${\displaystyle p={\frac {1}{v}}={\frac {dt}{dx}}}$

Each value on the trace p is the sum of all the samples along the line

${\displaystyle t=\tau +px}$

The transform is defined by:

${\displaystyle F(p,\tau )=\int _{-\infty }^{\infty }f(x,\tau +px)dx=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f(x,t)\delta (t-\tau -px)dxdt}$

The τ-p transform converts seismic records into a domain where all these events are separated. Simply put, each point in the τ-p domain is the sum of all the points in the x-t plane lying across a straight line with a slope p and intercept τ.[4] That also means a point in the x-t domain transforms into a line in the τ-p domain, hyperbolae transform into ellipses and so on. Similar to the Fourier transform, a signal in the τ-p domain can also be transformed back into the x-t domain.

During data acquisition, various effects have to be accounted for, such as near-surface structure around the source, noise, wavefront divergence and reverbations. It has to be ensured that a change in the seismic trace reflects a change in the geology and not one of the effects mentioned above. Deconvolution negates these effects to an extent and thus increases the resolution of the seismic data.

Seismic data, or a seismogram, may be considered as a convolution of the source wavelet, the reflectivity and noise.[5] Its deconvolution is usually implemented as a convolution with an inverse filter. Various well-known deconvolution techniques already exist for one dimension, such as predictive deconvolution, Kalman filtering and deterministic deconvolution. In multiple dimensions, however, the deconvolution process is iterative due to the difficulty of defining an inverse operator. The output data sample may be represented as:

${\displaystyle y({\underline {x}},t)=f({\underline {x}},t)**r({\underline {x}},t)}$

where ${\displaystyle f({\underline {x}},t)}$ represents the source wavelet, ${\displaystyle r({\underline {x}},t)}$ is the reflectivity function, ${\displaystyle {\underline {x}}}$ is the space vector and ${\displaystyle t}$ is the time variable. The iterative equation for deconvolution is of the form:

${\displaystyle f_{0}({\underline {x}},t)=\lambda y({\underline {x}},t)}$ and

${\displaystyle f_{n+1}({\underline {x}},t)=\lambda y({\underline {x}},t)+q(x,t)**f_{n}({\underline {x}},t)}$, where ${\displaystyle q(x,t)=\delta ({\underline {x}},t)-r({\underline {x}},t)}$

Taking the Fourier transform of the iterative equation gives:

${\displaystyle F_{n+1}({\underline {k}},\omega )=\lambda Y({\underline {k}},\omega )+F_{n}({\underline {k}},\omega )-\lambda F_{n}({\underline {k}},\omega )R({\underline {k}},\omega )}$

This is a first-order one-dimensional difference equation with index ${\displaystyle n}$, input ${\displaystyle \lambda Y({\underline {k}},\omega )u(n)}$, and coefficients that are functions of ${\displaystyle ({\underline {k}},\omega )}$. The impulse response is ${\displaystyle [1-\lambda R({\underline {k}},\omega )]^{n}u(n)}$, where ${\displaystyle u(n)}$ represents the one-dimensional unit step function. The output then becomes:

${\displaystyle F_{n}({\underline {k}},\omega )={\frac {Y({\underline {k}},\omega )}{R({\underline {k}},\omega )}}\lbrace 1-[1-\lambda R({\underline {k}},\omega )]^{n+1}\rbrace u(n)}$

The above equation can be approximated as

${\displaystyle F_{n}({\underline {k}},\omega )={\frac {Y({\underline {k}},\omega )}{R({\underline {k}},\omega )}}}$, if ${\displaystyle n\rightarrow \infty }$ and ${\displaystyle |1-\lambda R({\underline {k}},\omega )|<1}$

Note that the output is the same as the output of an inverse filter. An inverse filter does not actually have to be realized and the iterative procedure can be easily implemented on a computer.[6]

Stacking is another process used to improve the signal-to-noise ratio of the seismic profile. This involves gathering seismic traces from points at the same depth and summing them. This is referred to as "Common depth-point stacking" or "Common midpoint stacking". Simply speaking, when these traces are merged, the background noise cancels itself out and the seismic signal add up, thus improving the SNR.

Assuming a seismic wave ${\displaystyle s(x,z,t)}$ travelling upwards towards the surface, where ${\displaystyle x}$ is the position on the surface and ${\displaystyle z}$ is the depth. The wave's propagation is described by:

Evanescent and propagation regions for migration filter

${\displaystyle {\partial ^{2}s \over \partial x^{2}}+{\partial ^{2}s \over \partial z^{2}}={\frac {1}{c^{2}}}{\partial ^{2}s \over \partial t^{2}}}$

Migration refers to this wave's backward propagation. The two-dimensional Fourier transform of the wave at depth ${\displaystyle z_{0}}$ is given by:

${\displaystyle S(k_{x},z_{0},\Omega )=\int \int s(x,z_{0},t)e^{-j(\Omega t-k_{x}x)}dxdt}$

To obtain the wave profile at ${\displaystyle z=z_{0}}$, the wave field ${\displaystyle s(x,z,t)}$ can be extrapolated to ${\displaystyle s(x,z_{0}+\Delta z,t)}$ using a linear filter with an ideal response given by:

${\displaystyle H(\omega _{1},\omega _{2})={\begin{cases}e^{j{\sqrt {\alpha ^{2}{\omega _{2}}^{2}-{\omega _{1}}^{2}}}},&{\text{ for }}|\omega _{1}|<|\alpha \omega _{2}|\\0,&{\text{ else}}\end{cases}}}$

where ${\displaystyle \omega _{1}}$ is the x component of the wavenumber, ${\displaystyle k_{x}}$, ${\displaystyle \omega _{2}}$ is the temporal frequency ${\displaystyle \Omega }$ and

${\displaystyle \alpha ={\frac {1}{c}}{\frac {\Delta z}{\Delta t}}}$

For implementation, a complex fan filter is used to approximate the ideal filter described above. It must allow propagation in the region ${\displaystyle |\alpha \omega _{2}|>|\omega _{1}|}$ (called the propagating region) and attenuate waves in the region ${\displaystyle |\alpha \omega _{2}|<|\omega _{1}|}$ (called the evanescent region). The ideal frequency response is shown in the figure.[7]