Multidimensional seismic data processing

Last updated October 30, 2021

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

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:

$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 ${\underline {k}}$ is the spatial frequency (also known as wavenumber) and $\omega$ is the temporal frequency. The two-dimensional equivalent of the frequency domain is also referred to as the ${\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, $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.

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

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

$t=\tau +px$

The transform is defined by:

$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

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:

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

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

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

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

Taking the Fourier transform of the iterative equation gives:

$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 $n$ , input $\lambda Y({\underline {k}},\omega )u(n)$ , and coefficients that are functions of $({\underline {k}},\omega )$ . The impulse response is $[1-\lambda R({\underline {k}},\omega )]^{n}u(n)$ , where $u(n)$ represents the one-dimensional unit step function. The output then becomes:

$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

$F_{n}({\underline {k}},\omega )={\frac {Y({\underline {k}},\omega )}{R({\underline {k}},\omega )}}$ , if $n\rightarrow \infty$ and $|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

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

${\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 $z_{0}$ is given by:

$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 $z=z_{0}$ , the wave field $s(x,z,t)$ can be extrapolated to $s(x,z_{0}+\Delta z,t)$ using a linear filter with an ideal response given by:

$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 $\omega _{1}$ is the x component of the wavenumber, $k_{x}$ , $\omega _{2}$ is the temporal frequency $\Omega$ and

$\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 $|\alpha \omega _{2}|>|\omega _{1}|$ (called the propagating region) and attenuate waves in the region $|\alpha \omega _{2}|<|\omega _{1}|$ (called the evanescent region). The ideal frequency response is shown in the figure.^[7]

Related Research Articles

In engineering, a transfer function of a system, sub-system, or component is a mathematical function which theoretically models the system's output for each possible input. They are widely used in electronics and control systems. In some simple cases, this function is a two-dimensional graph of an independent scalar input versus the dependent scalar output, called a transfer curve or characteristic curve. Transfer functions for components are used to design and analyze systems assembled from components, particularly using the block diagram technique, in electronics and control theory.

In physics, a standing wave, also known as a stationary wave, is a wave which oscillates in time but whose peak amplitude profile does not move in space. The peak amplitude of the wave oscillations at any point in space is constant with respect to time, and the oscillations at different points throughout the wave are in phase. The locations at which the absolute value of the amplitude is minimum are called nodes, and the locations where the absolute value of the amplitude is maximum are called antinodes.

In mathematics, a Fourier transform (FT) is a mathematical transform that decomposes functions depending on space or time into functions depending on spatial or temporal frequency, such as the expression of a musical chord in terms of the volumes and frequencies of its constituent notes. The term Fourier transform refers to both the frequency domain representation and the mathematical operation that associates the frequency domain representation to a function of space or time.

In mathematics, deconvolution is the operation inverse to convolution. Both operations are used in signal processing and image processing. For example, convolution can be used to apply a filter, and it may be possible to recover the original signal using deconvolution.

Fourier optics is the study of classical optics using Fourier transforms (FTs), in which the waveform being considered is regarded as made up of a combination, or superposition, of plane waves. It has some parallels to the Huygens–Fresnel principle, in which the wavefront is regarded as being made up of a combination of spherical wavefronts whose sum is the wavefront being studied. A key difference is that Fourier optics considers the plane waves to be natural modes of the propagation medium, as opposed to Huygens–Fresnel, where the spherical waves originate in the physical medium.

The Short-time Fourier transform (STFT), is a Fourier-related transform used to determine the sinusoidal frequency and phase content of local sections of a signal as it changes over time. In practice, the procedure for computing STFTs is to divide a longer time signal into shorter segments of equal length and then compute the Fourier transform separately on each shorter segment. This reveals the Fourier spectrum on each shorter segment. One then usually plots the changing spectra as a function of time, known as a spectrogram or waterfall plot, such as commonly used in Software Defined Radio (SDR) based spectrum displays. Full bandwidth displays covering the whole range of an SDR commonly use Fast Fourier Transforms (FFTs) with 2^24 points on desktop computers.

In mathematics, the Weierstrass elliptic functions are elliptic functions that take a particularly simple form. They are named for Karl Weierstrass. This class of functions are also referred to as ℘-functions and they are usually denoted by the symbol ℘, a uniquely fancy script p. They play an important role in the theory of elliptic functions. A ℘-function together with its derivative can be used to parameterize elliptic curves and they generate the field of elliptic functions with respect to a given period lattice.

In mathematics, the Poisson summation formula is an equation that relates the Fourier series coefficients of the periodic summation of a function to values of the function's continuous Fourier transform. Consequently, the periodic summation of a function is completely defined by discrete samples of the original function's Fourier transform. And conversely, the periodic summation of a function's Fourier transform is completely defined by discrete samples of the original function. The Poisson summation formula was discovered by Siméon Denis Poisson and is sometimes called Poisson resummation.

In mathematics and in signal processing, the Hilbert transform is a specific linear operator that takes a function, $u (t)$ of a real variable and produces another function of a real variable $H(u)(t)$ . This linear operator is given by convolution with the function $. The Hilbert transform has a particularly simple representation in the frequency domain: It imparts a phase shift of \pm90° to every frequency component of a function, the sign of the shift depending on the sign of the frequency. The Hilbert transform is important in signal processing, where it is a component of the analytic representation of a real-valued signal u (t) . The Hilbert transform was first introduced by David Hilbert in this setting, to solve a special case of the Riemann-Hilbert problem for analytic functions.$

In the theory of stochastic processes, the Karhunen–Loève theorem, also known as the Kosambi–Karhunen–Loève theorem is a representation of a stochastic process as an infinite linear combination of orthogonal functions, analogous to a Fourier series representation of a function on a bounded interval. The transformation is also known as Hotelling transform and eigenvector transform, and is closely related to principal component analysis (PCA) technique widely used in image processing and in data analysis in many fields.

In system analysis, among other fields of study, a linear time-invariant system is a system that produces an output signal from any input signal subject to the constraints of linearity and time-invariance; these terms are briefly defined below. These properties apply to many important physical systems, in which case the response y(t) of the system to an arbitrary input x(t) can be found directly using convolution: y(t) = x(t) ∗ h(t) where h(t) is called the system's impulse response and ∗ represents convolution. What's more, there are systematic methods for solving any such system, whereas systems not meeting both properties are generally more difficult to solve analytically. A good example of an LTI system is any electrical circuit consisting of resistors, capacitors, inductors and linear amplifiers.

Geophysical survey is the systematic collection of geophysical data for spatial studies. Detection and analysis of the geophysical signals forms the core of Geophysical signal processing. The magnetic and gravitational fields emanating from the Earth's interior hold essential information concerning seismic activities and the internal structure. Hence, detection and analysis of the electric and Magnetic fields is very crucial. As the Electromagnetic and gravitational waves are multi-dimensional signals, all the 1-D transformation techniques can be extended for the analysis of these signals as well. Hence this article also discusses multi-dimensional signal processing techniques.

The Frank–Tamm formula yields the amount of Cherenkov radiation emitted on a given frequency as a charged particle moves through a medium at superluminal velocity. It is named for Russian physicists Ilya Frank and Igor Tamm who developed the theory of the Cherenkov effect in 1937, for which they were awarded a Nobel Prize in Physics in 1958.

The prolate spheroidal wave functions are eigenfunctions of the Laplacian in prolate spheroidal coordinates, adapted to boundary conditions on certain ellipsoids of revolution. Related are the oblate spheroidal wave functions.

Time-domain thermoreflectance is a method by which the thermal properties of a material can be measured, most importantly thermal conductivity. This method can be applied most notably to thin film materials, which have properties that vary greatly when compared to the same materials in bulk. The idea behind this technique is that once a material is heated up, the change in the reflectance of the surface can be utilized to derive the thermal properties. The reflectivity is measured with respect to time, and the data received can be matched to a model with coefficients that correspond to thermal properties.

In optics, the Fraunhofer diffraction equation is used to model the diffraction of waves when the diffraction pattern is viewed at a long distance from the diffracting object, and also when it is viewed at the focal plane of an imaging lens.

In digital signal processing, multidimensional sampling is the process of converting a function of a multidimensional variable into a discrete collection of values of the function measured on a discrete set of points. This article presents the basic result due to Petersen and Middleton on conditions for perfectly reconstructing a wavenumber-limited function from its measurements on a discrete lattice of points. This result, also known as the Petersen–Middleton theorem, is a generalization of the Nyquist–Shannon sampling theorem for sampling one-dimensional band-limited functions to higher-dimensional Euclidean spaces.

In signal processing, nonlinear multidimensional signal processing (NMSP) covers all signal processing using nonlinear multidimensional signals and systems. Nonlinear multidimensional signal processing is a subset of signal processing. Nonlinear multi-dimensional systems can be used in a broad range such as imaging, teletraffic, communications, hydrology, geology, and economics. Nonlinear systems cannot be treated as linear systems, using Fourier transformation and wavelet analysis. Nonlinear systems will have chaotic behavior, limit cycle, steady state, bifurcation, multi-stability and so on. Nonlinear systems do not have a canonical representation, like impulse response for linear systems. But there are some efforts to characterize nonlinear systems, such as Volterra and Wiener series using polynomial integrals as the use of those methods naturally extend the signal into multi-dimensions. Another example is the Empirical mode decomposition method using Hilbert transform instead of Fourier Transform for nonlinear multi-dimensional systems. This method is an empirical method and can be directly applied to data sets. Multi-dimensional nonlinear filters (MDNF) are also an important part of NMSP, MDNF are mainly used to filter noise in real data. There are nonlinear-type hybrid filters used in color image processing, nonlinear edge-preserving filters use in magnetic resonance image restoration. Those filters use both temporal and spatial information and combine the maximum likelihood estimate with the spatial smoothing algorithm.

The Fokas method, or unified transform, is an algorithmic procedure for analysing boundary value problems for linear partial differential equations and for an important class of nonlinear PDEs belonging to the so-called integrable systems. It is named after Greek mathematician Athanassios S. Fokas.

Tau functions are an important ingredient in the modern theory of integrable systems, and have numerous applications in a variety of other domains. They were originally introduced by Ryogo Hirota in his direct method approach to soliton equations, based on expressing them in an equivalent bilinear form. The term Tau function, or $-function$ , was first used systematically by Mikio Sato and his students in the specific context of the Kadomtsev–Petviashvili equation, and related integrable hierarchies. It is a central ingredient in the theory of solitons. Tau functions also appear as matrix model partition functions in the spectral theory of Random Matrices, and may also serve as generating functions, in the sense of combinatorics and enumerative geometry, especially in relation to moduli spaces of Riemann surfaces, and enumeration of branched coverings, or so-called Hurwitz numbers.

References

↑ Rector, James; Mangriotis, M. D. (2010). "Vertical Seismic Profiling". Encyclopedia of Solid Earth Geophysics. Springer. pp. 430–433. ISBN 978-90-481-8702-7.
↑ 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.
↑ Donati, Maria (1995). "Seismic reconstruction using a 3D tau-p transform" (PDF). CREWES Research Report. 7.
↑ 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.
↑ 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.
↑ Mersereau, Russell; Dudgeon, Dan. Multidimensional Digital Signal Processing. Prentice-Hall. pp. 350–352.
↑ Mersereau, Russell; Dudgeon, Dan. Multidimensional Digital Signal Processing. Prentice-Hall. pp. 359–363.

External links

This page is based on this Wikipedia article
Text is available under the CC BY-SA 4.0 license; additional terms may apply.
Images, videos and audio are available under their respective licenses.

[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.

[1]

[2]

[3]

[4]

[5]

[6]

[7]