Multidimensional seismic data processing

Multidimensional seismic data processing forms a major component of seismic profiling, a technique used in geophysical exploration. The technique itself has various applications, including mapping ocean floors, determining the structure of sediments, mapping subsurface currents and hydrocarbon exploration. Since geophysical data obtained in such techniques is a function of both space and time, multidimensional signal processing techniques may be better suited for processing such data.

Data acquisition

Offset VSP

There are a number of data acquisition techniques used to generate seismic profiles, all of which involve measuring acoustic waves by means of a source and receivers. These techniques may be further classified into various categories,[1] depending on the configuration and type of sources and receivers used. For example, zero-offset vertical seismic profiling (ZVSP), walk-away VSP etc.

The source (which is typically on the surface) produces a wave travelling downwards. The receivers are positioned in an appropriate configuration at known depths. For example, in case of vertical seismic profiling, the receivers are aligned vertically, spaced approximately 15 meters apart. The vertical travel time of the wave to each of the receivers is measured and each such measurement is referred to as a “check-shot” record. Multiple sources may be added or a single source may be moved along predetermined paths, generating seismic waves periodically in order to sample different points in the sub-surface. The result is a series of check-shot records, where each check-shot is typically a two or three-dimensional array representing a spatial dimension (the source-receiver offset) and a temporal dimension (the vertical travel time).

Data processing

The acquired data has to be rearranged and processed to generate a meaningful seismic profile: a two-dimensional picture of the cross section along a vertical plane passing through the source and receivers. This consists of a series of processes: filtering, deconvolution, stacking and migration.

Multichannel filtering

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]

Two-dimensional Fourier transform design

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.

τ-p transform design

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.

Deconvolution

See also: Deconvolution

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

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.

Migration

See also: Seismic migration

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]

References

1. Rector, James; Mangriotis, M. D. (2010). "Vertical Seismic Profiling". Encyclopedia of Solid Earth Geophysics. Springer. pp. 430–433. ISBN 978-90-481-8702-7.
2. Tatham, R; Mangriotis, M (Oct 1984). "Multidimensional Filtering of Seismic Data". Proceedings of the IEEE. 72 (10): 1357–1369. doi:10.1109/PROC.1984.13023.
3. Donati, Maria (1995). "Seismic reconstruction using a 3D tau-p transform" (PDF). CREWES Research Report. 7.
4. McMechan, G. A.; Clayton, R. W.; Mooney, W. D. (10 February 1982). "Application of Wave Field Continuation to the Inversion of Refraction Data" (PDF). Journal of Geophysical Research. 87: 927–935. doi:10.1029/JB087iB02p00927.
5. Arya, V (April 1984). "Deconvolution of Seismic Data - An Overview". IEEE Transactions on Geoscience Electronics. 16 (2): 95–98. doi:10.1109/TGE.1978.294570.
6. Mersereau, Russell; Dudgeon, Dan. Multidimensional Digital Signal Processing. Prentice-Hall. pp. 350–352.
7. Mersereau, Russell; Dudgeon, Dan. Multidimensional Digital Signal Processing. Prentice-Hall. pp. 359–363.

External links

• Tau-P Processing of Seismic Refraction Data
• Reflections on the Deconvolution of Land Seismic Data
• Seismic profiling
• COMMON-MIDPOINT STACKING