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.
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, 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.
Two-dimensional Fourier transform design.
The two-dimensional Fourier transform is defined as:
formula_1
where formula_2 is the spatial frequency (also known as wavenumber) and formula_3 is the temporal frequency. The two-dimensional equivalent of the frequency domain is also referred to as the formula_4 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, formula_5. 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.
formula_6
Each value on the trace "p" is the sum of all the samples along the line
formula_7
The transform is defined by:
formula_8
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 "τ". 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. 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:
formula_9
where formula_10 represents the source wavelet, formula_11 is the reflectivity function, formula_12 is the space vector and formula_13 is the time variable. The iterative equation for deconvolution is of the form:
formula_14
and
formula_15, 
where
formula_16
Taking the Fourier transform of the iterative equation gives:
formula_17
This is a first-order one-dimensional difference equation with index formula_18, input formula_19, and coefficients that are functions of formula_20. The impulse response is formula_21, where formula_22 represents the one-dimensional unit step function. The output then becomes:
formula_23
The above equation can be approximated as
formula_24,
if formula_25 and formula_26
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.
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 formula_27 travelling upwards towards the surface, where formula_28 is the position on the surface and formula_29 is the depth. The wave's propagation is described by:
formula_30
Migration refers to this wave's backward propagation. The two-dimensional Fourier transform of the wave at depth formula_31 is given by:
formula_32
To obtain the wave profile at formula_33, the wave field formula_27 can be extrapolated to formula_35 using a linear filter with an ideal response given by:
formula_36
where formula_37 is the x component of the wavenumber, formula_38, formula_39 is the temporal frequency formula_40 and
formula_41
For implementation, a complex fan filter is used to approximate the ideal filter described above. It must allow propagation in the region formula_42 (called the propagating region) and attenuate waves in the region formula_43 (called the evanescent region). The ideal frequency response is shown in the figure.