Note: Descriptions are shown in the official language in which they were submitted.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
Processing Seismic T)ata
The present invention relates to a method of processing seismic data, in
particular to a
method of interpolating from acquired seismic data traces to produce an
estimated trace
for a location fox which no actual data trace is available.
In seismic exploration, acoustic signals produced by a seismic source travel
downwards
into the earth's interiox. The acoustic energy is partially reflected by
features within the
earth that include an acoustic impedance boundary, and the reflected seismic
energy is
detected by seismic receivers which contain ane or more sensors for detecting
acoustic
energy. The reflected acoustic energy is processed in order to obtain
information about
the nature of the earth's interior at the survey location. For example, the
received
acoustic signals carry information relating to the structure of reflective
layers within the
earth's interior, such as faults or boundaries between different types of
rocks.
A seismic receiver generally acquires a digital record of the acoustic energy
incident on
the receiver, by sampling the incident acoustic energy at regular sampling
intervals.
The digital signals acquired by a receiver are normally referred to as
"traces", and
provide a record of, for example, the amplitude of the acquired acoustic
energy or a
component thereof.
Figure 1 is a schematic illustration of a seismic arrangement. A seismic
energy source
(acoustic source) 1 and a receiver 2 are located on the surface 3 of the earth
4. When
the source 1 is actuated it emits acoustic energy downwards into the earth,
and one path
of acoustic energy is shown by reference 5.
Figure 1 shows one geological feature 7 that acts as a partial reflector of
acoustic
energy. Acoustic energy that propagates along the downwards path 5 is
reflected
upwards along path 6 by the reflector 7, and is incident on the receiver 2. In
figure 1, it
is assumed that the reflector 7 is horizontal, so that the point of reflection
S lies directly
below the mid point 9 on the earth's surface between the source 1 and the
receiver 2.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
2
The time between emission of a pulse of acoustic energy from the source l and
the
arrival of its reflection at the receiver 2 is measured directly. This time is
the travel time
of acoustic energy along the downwards path 5 and the upwards path 6.
A seismic survey will generally contain an array of source positions and an
array of
receivers. It may occur that a point on a geological structure acts as a
reflection point
for more than one source-receiver pair, and this is illustrated in figure 1.
Figure 1
shows a second acoustic source 1' and a second receiver 2', and also shows
paths 5',6'
of acoustic energy from the source 1' to the receiver 2' via a reflection at
the geological
feature 7. The second source 1' and the second receiver 2' are disposed such
that the
two sources 1, 1' and the two receivers 2, 2' lie on a straight line, and so
that the
midpoint on the earth's surface between the second source 1' and the second
receiver 2'
is coincident with the midpoint 9 between the first source 1 and the first
receiver 2. The
midpoint 9 is therefore referred to as the "common mid point" or CMP for the
two
source-receiver pairs l, 2 and 1', 2'. It will also be seen that acoustic
energy that travels
from the second acoustic source 1' to the second receiver 2' along the paths
5',6' will
again undergo reflection at the reflection point 8 - the point 8 is a common
reflection
point for acoustic energy travelling from the source 1 to the receiver 2 and
for acoustic
energy travelling from the second source 1' to the second receiver 2'.
In figure 1, the geological feature 7 that acts as a reflector of acoustic
energy is assumed
to be substantially horizontal so that the point of reflection 8 Iies
vertically below the
mid point 9. However, for situations in which the inclination of the reflector
is
significant, this assumption is inadequate. This is illustrated in figure 2,
which
corresponds to figure 1 but shows a "dipping" reflector (that is, a reflector
that is
inclined to the horizontal).
As shown in figure 2, for a dipped reflector 7, when the acoustic source 1 is
actuated
acoustic energy will again be reflected towards the receiver 2. However, the
length of
the corresponding normal travel path 10 in figure 2 will be different from the
length of
the normal travel path 10 in figure 1. Furthermore, the surface location 9
corresponding
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
3
to the normal travel path 10 is not located at the mid point between the
source l and the
receiver 2, as a result of the dip of the reflector.
In a typical three-dimensional seismic survey, sources 1 and receivers 2 are
located so
as to cover an area of interest, for example as shown in figure 3. The sources
1 are
arranged as an array, which may be regular if the surface enviromnent at the
survey
location permits this. The receivers 2 may similarly be arranged as a surface
array,
which may again be regular if this is possible. The sources and receivers are
arranged
such that the target sub-surface is adequately illuminated by acoustic energy
so that an
accurate three-dimensional subsurface image may be obtained for the survey
location.
A survey location is usually represented by Cartesian co-ordinates x and y in
metric
unit. To facilitate further processing, a survey location is also divided into
grid cells or
"bins", and this enables the survey location to be represented by grid
numbers. Each
trace acquired in the survey is allocated to one of the bins, and this
allocation process is
known as "binning".
Figure 4 shows a portion of a typical grid of bins 11 (in practice, a bin
array will contain
many more bins than shown in Figure 4). The bins 11 are denoted as squares in
figure
4, and in practice are generally square (although in principle the bins do not
need to be
square). The length and width of a bin are typically of the order of tens of
metres.
Figure 4 also illustrates the location of acquired traces. In Figure 4 the
traces have been
allocated to bins according to the mid-point of each trace, so Figure 4 shows
a grid of
common mid-point (CMP) bins 11. (Traces may be binned according to parameters
other than their respective mid-points.) Each solid circle 12 in figure 4
represents the
mid-point of an acquired data trace. The mid-point location of a trace is
defined as the
midpoint of the source and the receiver used to acquire that trace. In the
case of a flat
reflector the reflection point associated with a trace is directly below the
mid-point of
the trace, although this is not true for a dipping reflector.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
4
A trace is described by "trace attributes". Trace attributes may contain
information, for
example, about the midpoint co-ordinates (x,y) of the trace, about its bin
information
(e.g.. bin ordinal, bin centre co-ordinate etc). Other trace attributes
include, for
example, the source location and the receiver location, the source-receiver
separation
("offset"), and the source receiver azimuth.
The locations of the mid-points of the traces acquired at the survey location
are
determined by the arrangement of sources and receivers. Ideally, each mid-
point bin
would contain at least one trace, but this does not always happen - in figure
4, it will be
seen that the central bin does not contain a trace.
Ideally, each bin would contain the same number of traces, but it will be seen
that this is
also not the case for figure 4. While most bins contain one trace, a few bins
contain two
traces (and, as noted above, one bin does not contain a trace). The number of
traces in a
bin is referred to as the "fold" of the bin.
The conventional approach to the situation of figure 4 in which a bin does not
contain a
trace is the "copy and move" technique. In this prior art method, a trace from
a
neighbouring bin is copied to the vacant bin. This process is unsatisfactory,
however,
since it necessarily results in a loss of resolution of the subsurface image
finally
obtained from the seismic data.
The present invention provides a method of processing seismic data comprising
the step
of: obtaining an estimated data trace for a pre-determined target location by
interpolation from a set of spatially irregularly sampled seismic data traces;
wherein the
interpolation step comprises applying a weighting function dependent on the
target
location.
The weighting function may be dependent on the value at the target location of
the local
input density in a co-ordinate direction. It may be proportional to a sinc
function having
a locally varying sinc frequency, where sinc (x,f) = sin (fx)/(fx), and it may
be
proportional to a sinc function of the general form sinc (x,sf), where s (s <
1) is a scalar.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
The invention may be used to estimate a trace for a vacant bin, in which case
the target
location identifies the vacant bin for example by its centre co-ordinates.
Alternatively,
the method may comprise obtaining a plurality of estimated seismic data
traces, each for
a respective one of a plurality of pre-determined spatial co-ordinates, from
the first set
of data traces. This embodiment may be used to convert the data to a new bin
array.
The method may further comprise the step of estimating at least one attribute,
for
example such as the source-receiver offset or the source-receiver azimuth, or
an
attribute relating to the interpolation process used, for the or each
estimated data trace.
The at least one attribute may be estimated for a selected one of the
estimated data
traces using the same weighting function as used in the step of estimating the
selected
estimated data trace.
The method may further comprise obtaining a plurality of estimated seismic
data traces,
each for one of the plurality of pre-determined spatial co-ordinates, from
another set of
spatially irregularly sampled seismic data traces. This enables two sets of
data traces to
be regularised to a common bin array, thereby allowing the two sets of data
traces to be
combined and/or compared. This is of use in a 4-D survey.
A second aspect of the present invention provides a method of seismic
surveying
comprising acquiring at least one set of spatially irregularly sampled seismic
data traces,
and processing the data traces according to a method as defined above.
A third aspect of the present invention provides an apparatus for processing
seismic data
comprising: means for obtaining an estimated data trace for a pre-determined
target
location by interpolation from a set of spatially irregularly sampled seismic
data traces
using a weighting function dependent on the target location.
The apparatus may comprise a programmable data processor.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
6
A fourth aspect of the invention provides a storage medium comprising a
program for a
data processor of an apparatus as defined in above.
A fifth aspect of the invention provides a storage medium containing a program
for
controlling a data processor to perform a method as defined above.
Further aspects and preferred features of the invention are defined in the
appended
claims.
Preferred embodiments of the present invention will now be described by way of
illustrative example with reference to the accompanying figures in which:
Figure 1 is a schematic illustration of a seismic survey involving a
horizontal reflector;
Figure 2 is a schematic illustration of a seismic survey involving a dipped
reflector;
Figure 3 is a schematic plan view of a multi-source, multi-receiver seismic
survey;
Figure 4 illustrates typical results for distribution of traces within CMP
bins;
Figure 5(a) illustrates a first embodiment of the present invention;
Figure 5(b) illustrates a second embodiment of the present invention;
Figure 6 illustrates the sinc function;
Figure 7 is a schematic flow diagram illustrating a method of the present
invention;
Figures 8(a), 8(b) and 8(c) illustrate three possible ways of defining the
neighbourhood
over which interpolation is carried out; and
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
7
Figure 9 is a block schematic diagram of an apparatus according to the present
invention.
Preferred embodiments of the invention will be described with reference to mid-
point
bins as shown in Figure 4, although the invention is not limited to use with
mid-point
bins.
In the present invention, a seismic data trace for pre-determined coordinates
is estimated
by interpolation from a plurality of acquired seismic data traces. The
acquired seismic
data traces represent an irregular spatial sampling of the survey location, so
that the
traces are irregularly distributed among the mid-point bins as shown in figure
4.
The present invention may be used to regularise a set of acquired seismic data
traces to
a pre-selected bin arrangement, and this embodiment is shown schematically in
figure
5(a). The left hand side of Figure 5(a) shows the original set of acquired
data traces
and, as noted above, the original set of traces have an irregular spatial
distribution of
their mid points through the bin array.
In this embodiment, the user defines a grid - that is, an array of mid-point
bins - and this
will normally be a two-dimensional array of bins. The user-defined mid-point
bin array
may be the same as, or different to, the original mid-point bin array of the
acquired data
traces. The invention is then employed to estimate a trace for bins of the
user defined
bin array. In general, a user will desire to estimate a trace for a particular
output
location, such as the bin centre, for each bin of the user-defined bin array
(although if
there is a bin of the user-defined bin array for which there is already a
trace at the exact
output location it would be possible to retain the existing trace rather than
estimate a
trace for that bin).
The result, as shown in the right hand view of figure 5(a) is one estimated
tarace 13 for
each bin 11, with each trace preferably being centared in its respective bin.
This
embodiment is particularly useful where it is desired to compare andlor merge
the
results of two seismic surveys carried out at the same survey location.
Different bin
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
grids may have been defined for the two surveys. Furthermore, if the two
surveys were
carried out using different source and/or receiver arrays, the two surveys may
have a
different distribution of trace mid-points even if they have the same defined
mid-point
bin grid. In either of these cases it would be difficult simply to combine or
compare the
two sets of data traces. The invention, however, makes it possible to
regularise the mid-
point bin array of one survey to the bin array of the other survey.
Alternatively, each of
the surveys may be regularised to a common pre-determined mid-point bin array
which
is different from the original mid-point bin array of either survey.
Each trace is preferably centred in its associated mid-point bin. This fully
regularises
the data, as required by some processing software:
The embodiment of 5(a) is also useful for data comparison in 4D studies.
Furthermore,
complete regularisation of traces to regularly-spaced CMP bin grid centres,
also fulfils
the requirements of some forms of pre-stack migration.
Figure 5(b) illustrates a second embodiment of the invention, in which the
acquired
seismic traces are partially regularised. The left hand side of Figure 5(b) is
the same as
the left hand sides of Figure 5(a), and shows the original set of acquired
data traces as
binned into a mid-point bin array. In this embodiment, the original acquired
traces 12
are retained, and the invention is simply used to obtain an estimated trace 13
for any
vacant bin. The result is that each bin 11 either contains an estimated trace
13 or
contains one or more actual data traces 12. This embodiment of the invention
may be
useful where it is desired to apply Kirchhoff pre-stack migration.
It is generally preferable fox each bin to contain a trace, so that the fold
of the data is
regularised. It is however possible for one or more bins to be empty, for
example for a
part of the survey location where the input data density is too low. In
principle, it would
be possible to estimate more than one trace for a bin, by specifying two or
more target
locations within a bin.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
9
In a variant of the embodiment of figure 5(b), if a bin contains two or more
acquired
data traces then one of the traces is retained. For example, the trace with a
mid-point
closest to the bin centre may be retained, and other traces in that bin
discarded.
Alternatively, a stack trace derived from all traces in a bin may be output
rather than the
trace closest to the bin centre. In these embodiments, the result is a grid in
which each
bin contains one data trace (which will be an estimated data trace for a bin
that was
originally a vacant bin). .
According to the invention, the interpolation to obtain a trace for a target
location is
carried out using a weighting function that depends on the target location.
The
weighting function is dependent on the local input cluster density - that is,
the density
of traces - at the target location. The invention thus provides an adaptive
interpolation
method.
In a particularly preferred embodiment, the interpolation is carried out
according to the
following equation:
~D(x~,y~,t-~t(x~ -x;,y~ -y~,dip~))~W~x~ -x~l,ly~ -ytlafxt~f,,t)
D~~x~~y~~t~= J
~~ x~ -xil~ly~ -y~I~fxa,fyr
(1)
In equation (1), D'(xi, y1, t) is the estimated data trace at time t at target
coordinates (xi,
yt). The target coordinates may be defined by a user, and may represent, for
example,
the central coordinates of a bin for which it is desired to estimate a trace.
The bin may
be a vacant bin (for example in the embodiment of Figure 5(b)) or it may
already
contain an acquired data trace (as in the embodiment of Figure 5(a)). The
coordinates
(x~, y~) represent the coordinates of the j'h existing trace data trace D. The
term dips is
the dip associated with the j'h existing trace. W is a weighting function, and
is used to
weight some~traces in the summation more highly than others.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
The weighting function W is centred at the target coordinates. According to
the
invention, the weighting function is dependent on the target co-ordinates. One
suitable
weighting function is a weighting function that is proportional to a sinc
function, where
sinc (x,f) = sin (fx)/(fx), having a locally varying sinc frequency. The
weighting
function may alternatively be proportional to a product of a sinc function in
the primary
(x-) direction with a sinc function in the secondary (y-) direction. The sinc
frequency
of the or each sinc function may depend on the target location.
The summation of equation (1) is carried out over a user-defined area that
surrounds the
target co-ordinates. Conveniently, the user-defined area is rectangular. The
summation
is carried out using only one sample from each input trace. If the reflector
is horizontal
so the dip is zero, the summation is carried out on a flat time slice. Dipping
events may
be tracked by shifting input samples before the interpolation.
As can be seen from equation (1), reflector dip is taken account of by
correcting each
acquired trace used in the summation for reflector dip. Dipping events axe
tracked by
shifting input samples before the interpolation is performed, and this is
equivalent to
mapping steeply dipping events into lower-spatial frequencies before the
interpolation
(essentially by warping the 2D surface on which the interpolation is performed
so that it
follows the dipping energy). Dipping events may also be partially flattened by
performing a normal moveout (NMO) correction before performing the
interpolation; if
this is done, inverse NMO is preferably applied after the interpolation.
The dip may be determined in any suitable way. It may be determined from the
data
traces, or it may come from prior information about the survey location. For
example,
dips in the primary and secondary directions may be scanned separately for
every
sample at every output location. The best dip is picked based on semblance
values
calculated from a user-specified spatial and time window. The dipping events
are
tracked by shifting input samples prior to the interpolation.
In a preferred embodiment of the invention, the weighting function W is given
by:
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
I1
W~.xj -'.xi~,~yj yil,fxi,fyi)-SlnC~.xj -JCi~, fxi)'slnC~yj -yilofyi)
In equation (2), "sinc (x,f)" denotes "sin (xf)/(xf)". The parameters fxi, and
fyi are the
sinc frequencies, and will be further described below.
In this embodiment, the weighting function is the product of a sinc function
in the x-
dimension and a sinc function in the y-dimension, with both functions being
centred on
the target coordinates. The sinc frequency fxi, fyi of each sinc function is
dependent on
the target location,
The general form of a sinc function is illustrated in figure 6, and it will be
seen that the
greatest magnitude of this function occurs for values of x close to zero.
Thus, the
weighting function W of equation (2) gives most weight to traces acquired in
bins in the
near neighbourhood of the vacant bin.
As noted above, a seismic data trace is generally be labelled with one or more
"attributes", which provide information about the acquisition arrangements for
the trace.
Examples of trace attributes are the offset between the source and the
receiver used to
acquire the trace, and the source-receiver azimuth for the trace (that is, the
azimuthal
direction from the source to the receiver). Each acquired data trace 12 shown
in figure
4 will have accompanying trace attributes such as offset and azimuth. In a
case where
the trace location corresponds to the mid-point of the trace, knowledge of the
offset and
azimuth of a trace enables the source and receiver positions to be determined.
In addition to conventional trace attributes mentioned above, an interpolated
trace has
additional attributes that relate to the interpolation process used to obtain
the trace.
Such additional trace attributes include QC (quality control) literals output
to
distinguish the quality of the interpolated traces, for example to identify
those from an
area where data is sparse, and to enable future trace selection. They include,
for
example, (1) the distance from an interpolated trace to the nearest input
trace, (2) the
number of traces used to interpolate a particular trace, (3) the total weight
used in the
interpolation, (4) the gap size Dx in the primary direction, (5) the gap size
Dy in the
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
12
secondary direction, (6) the sine frequency f,~, in the x-direction and (7)
the sinc
frequency fyi in the y-direction.
Trace attributes, or information derived from the trace attributes, may be
required in
subsequent steps of processing the seismic data. It is therefore preferable to
provide an
estimated trace obtained by the method of the invention with one or more trace
attributes such as offset and azimuth. In a preferred embodiment of the
invention,
therefore, the trace attributes for the estimated trace at the target co-
ordinates are also
estimated. This is again done by a weighted interpolation from the attributes
of the
traces used in the interpolation. In a particularly preferred embodiment, the
weighting
function used in the interpolation to determine the attributes of an estimated
trace is the
same weighting function used in the interpolation to estimate the trace. For
example,
the offset of the trace estimated for a target location, may be determined
using:
~O~x;,y~))~W~x~ -xi~,ly~ -yil, fXi, fyi)
O~(x~Y)=
~W x; -x;I~IY; -Yil~fxi~
fyi
.l
In equation (3), O(xi, y1) is the offset of the i'jz seismic data trace (at
coordinates (xi,yi))
and O'(x, y) is the estimated offset of the trace estimated for the target
coordinates (x,y).
The weighting function W is the same weighting function as used in the
determination
of the estimated trace - that is, W is the same weighting function as used in
the
interpolation of equation (1). Other attributes of the trace such as the
azimuth may be
determined using a corresponding equation.
Attributes of the input traces are preferably normalised before the estimated
attribute is
found by interpolation. For example, azimuths of input traces may be adjusted
by ~ 180
degrees, and offset may be adjusted by using absolute values.
The sinc frequencies fxi and fyi are calculated individually for each set of
target
coordinates. The local input density at a target location is assessed in terms
of the gap
size in the primary direction, denoted as DX, and as the gap size in the
secondary
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
13
direction denoted as DY. (The "primary direction" and "secondary direction"
are. the x-
direction and the y-direction respectively. These terms are well-known in
seismic
surveying. For example, in a towed marine survey the "primary direction" is
usually
defined to be parallel to the streamer direction.)
fxl, which is the sine frequency in the primary direction at the target
location (xi, yi), is
then defined as I/DX, and the sine frequency fyi (the sine frequency in the
secondary
direction at the target location (xi, yt)) is defined as 1/DY. The maximum
value for each
sine frequency is the Nyquist frequency which is the reciprocal of the grid
spacing.
DX is the gap size in the primary (x-) direction at the taxget location. This
is the
distance in the primary direction between the nearest trace on one side of the
target
location and the nearest trace on the opposite side of the target location.
Similarly, DY
is the gap size in the secondary (y-) direction at the target location and is
the distance in
the secondary direction between the nearest trace on one side of the target
location and
the nearest trace on the opposite side of the target location. The gap size
DX, DY, and
thus the sine frequency fxl, fy~ will vary from one target location to
another, so that the
weighting function W is dependent on the taxget location. The gap sizes DX and
DY
are indicated in Figure 4 for the central vacant bin.
When it is difficult to calculate a reasonable gap size DX or DY, for example,
at the
edge of the data set, an artificial value is assigned to DX or DY based on the
distance
from the target location to the nearest traces.
Figure 7 is a schematic flow diagram illustrating one embodiment of the
present
invention.
Initially, at step 1, a set of seismic data traces is acquired. At step 2, the
traces are
binned (i.e., sorted into CMP bins). Steps 1 and 2 are conventional steps, and
the result
of step 2 will be generally similar to figure 4.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
14
The invention may alternatively be applied to pre-existing seismic data. In
this case,
step 1 is omitted and is replaced by the step, step 3, of retrieving pre-
existing seismic
data from storage. If necessary, the retrieved seismic data are binned.
However, if the
pre-existing seismic data had been binned before storage, so that the
retrieved seismic
data are already in the form shown in figure 4, it is not necessary to repeat
step 2.
At step 4, the user defines a desired bin array. The bin size and centres of
each bin are
defined at this stage. The bin array defined at step 4 may be the same as, or
may be
different from, the bin array used at step 2 or the bin array of pre-existing,
binned
seismic data. In the "infill" embodiment of Figure 5(b), for example, the bin
array
defined at step 4 would be the same as the existing bin array of the seismic
data traces.
Where the data traces are being regularised to a new bin array, however, the
bin array
defined at step 4 will be the desired new bin array, not the existing bin
array of the
seismic data traces. (Of course, if the desired bin array is identical to the
existing bin
array then step 4 may be omitted.)
At step 5, one of the bins defined at step 4 is selected. The centre of this
bin (or other
desired position in this bin) represent a target location for which it is
desired to estimate
a trace. At step 6, the sinc frequencies fxl, fyi for the target location are
determined, as
the reciprocals of the gap size in the x- and y-directions at the target
location.
At step 7, a seismic data trace is estimated for the target location defined
at step 5. This
is done according to equation (1) above, using a local weighting function W
that
depends on the target location. The weighting function W is preferably the
function
defined in equation (2).
At step ~, one or more attributes, for example such as the source-receiver
offset and/or
azimuth, are estimated for the trace estimated in step 7. This estimation step
is
performed by interpolation from the attributes from the input traces used in
the
interpolation of step 7, according to an equation having the general form of
equation (3)
above. The weighting function W used in step ~ is the same weighting function
as was
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
used in step 7. The attributes of the input traces may be normalised before
the
interpolation is performed.
At step 9, the trace determined in step 7 and the associated trace attributes
determined in
step 8 are output. The trace and attributes may be output in any suitable form
- for
example, they may be displayed on a computer monitor, may be output as a
digitised
trace in a form suitable for further processing, etc.
At step 10 it is determined whether a trace has been estimated for each bin
defined at
step 4 for which a trace is required to be estimated. The program loops over
the output
bin array. It first decides whether to estimate a new trace for a bin. For the
regularisation mode of Figure 5(a) a new trace is usually estimated, unless
there is an
existing trace at the exact output location. For the infill mode of Figure
5(b), a trace is
estimated when the output location is at a vacant bin. If step 10 yields a
"no"
determination, the next bin is selected at step 11, and steps 6-9 are repeated
for the next
bin. If there is still a "no" determination at step 10, steps 11 and 6-9 are
then repeated
as often as necessary until a "yes" determination is obtained at step 10. At
this point, a
trace has been estimated, by interpolation, for each desired bin of the grid
defined at
step 5, as shown in the right hand view of Figure 5(a) or Figure 5(b). These
results may
then be used in any desired way. For example, any conventional processing step
may
be applied to the results obtained by a method of figure 7. The results
obtained by a
method of Figure 7 may alternatively be combined with other traces acquired at
the
same survey location, and the combined data may be subjected to further
processing
steps.
In another embodiment, a method according to figure 7 may be further applied
to a
second set of seismic data traces acquired at the same survey location as the
first set of
seismic data traces. This process regularises both sets of seismic data to a
common bin
array, so that the two data sets may be combined.
In figure 7, the output step 9 is shown as being performed separately for each
trace
estimated for a bin. The invention is not limited to this, and step 9 could be
replaced by
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
16
this step of storing the trace and trace attributes estimated for the selected
bin. If
desired, the results for all estimated traces and trace attributes could be
output in a
single output step, following a "yes" determination at step 10 (although in
practice the
amount of memory required to store the results for all bins is large and,
depending on
the available machine memory, it may be necessary to output the trace and
trace
attributes for a group of bins rather than for every bin).
Figure 7 shows step 8 of determining the estimated attributes(s) for the
selected bin
location as being carried out after the step of determining the trace for the
selected bin
location. In principle, however, step 8 could be carried out before, or in
parallel with,
step 7. Furthermore, in principle, it would be possible to determine the
estimated trace
for each bin location, and subsequently determine the estimated trace
attributes) for
each bin.
The invention has been described with reference to embodiments in which each
trace is
defined by spatial co-ordinates that identify the mid point of the trace. The
invention is
not limited to this however, and may be applied to, for example, common
conversion-
point bins, common offset bins, common shot-point bins, common receiver-point
bins,
etc, provided that a trace for a target location can be estimated from traces
associated
with location surrounding the target location.
The invention has been described with reference to pre-stack seismic data. The
invention may, however, be also applied to post-stack data. Post-stack data
may be
treated as a special case of pre-stack data, with zero offset. (In general, in
the stacking
process two corrections are initially made to eliminate the effects of source-
receiver
offset - normally an NMO (normal moveout) correction is applied, and this is
followed
by a DMO (dip moveout) correction to compensate for mis-positioning owing to
any
dip of the reflecting interface. Application of NMO and DMO corrections
produces
traces which simulate the recording of seismic data with the source and
receiver at the
same location (that is, the traces simulate zero offset traces). It is then
possible to stack
traces that have the same or similar locations ("stacking" trace comprises
summing the
traces and dividing the amplitude of the resultant trace by the number of
traces
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
17
summed). This stacking step improves the signal-to noise ratio of the data,
since it
attenuates random noise, and. it also reduces the number of traces that
require
subsequent processing. In the example of figure 1, for example, a trace
obtained at
receiver 2 following actuation of source 1 could be stacked with a trace
obtained at
receiver 2' after actuation of source 1', since these traces have the same CMP
location.)
The invention may also be used to estimate trace attributes for an estimated
post-stack
trace. It should be noted, however, that a post-stack trace has fewer
independent
attributes than a pre-stack trace. For example, source-receiver offset, source-
receiver
azimuth and source or receiver location have no meaning for a post-stack
trace.
In a particularly preferred embodiment, the neighbourhood over which the
interpolation
is carried out is a rectangle having dimensions that are adaptively determined
based on
the location of the nearest traces on either side of the target location in
the primary (x-)
direction, and on the location of the nearest trace on either side of the
target location in
the secondary (y-) direction. In this embodiment, the neighbourhood over which
the
interpolation is performed is a rectangle having dimension in the primary (x-)
direction
of (2 x XLEN + 8x1 + 8X2), where XLEN is a user specified parameter, 8x1 is
the
distance in the primary (x-) direction between the target location and the
nearest trace to
the target location in the positive x-direction, and 8x2 is the distance in
the primary (x-)
direction between the target location and the nearest trace to the target
location in the
negative x-direction. Similarly, the dimension in the secondary (y-) direction
is (2 x
YLEN + byl + 8y2), where YLEN is a user specified parameter, 8y1 is the
distance in
the secondary (y-) direction between the target location and the nearest trace
to the
target location in the positive y-direction, and bye is the distance in the y-
direction
between the target location and the nearest trace to the target location in
the negative y-
direction. This is shown in Figure 8(c).
It should be noted, however, that the invention is not limited in this way,
and that the
neighbourhood over which the interpolation is performed may be defined in any
convenient way. For example, the neighbourhood over which the interpolation is
carried out may simply be defined as a rectangle having dimensions 2 x XT FN
in the
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
18
primary (x-) direction and 2 x YLEN in the secondary (y-) direction, where
XI,EN and
YLEN are user specified parameters. This is shown in Figure 8(a).
Alternatively, the neighbourhood over which the interpolation is carried out
may be a
rectangle having dimensions that are based on the location of the nearest
trace to the
target location. The dimension in the primary (x-) direction is 2 x (~I EN +
8x), where
~x is the difference in the x-direction between the target location and the
nearest trace to
the target location. The dimension along the secondary (y-) direction is
similarly 2 x
(YLEN + 8y), where 8y is the difference in the y-direction between the target
location
and the nearest trace to the target location. This is shown in Figure 8(b).
In the methods of Figures 8(b) and 8(c), the user preferably specifies upper
limits on the
overall dimensions in the x- and y-direction of the neighbourhood over which
the
interpolation is performed.
In the case of very low density data, it is possible that there will be no
trace, or only a
single trace, in a neighbourhood defined for the interpolation step. If there
is no trace
captured in an interpolation neighbourhood, it is possible either to leave the
vacant bin
empty, or to infill the bin with a zero amplitude trace. If there is only one
trace captured
in an interpolation neighbourhood, it is possible to copy the single captured
trace to the
output bin and flag it as 'copied', to leave the vacant bin empty, or to
infill the vacant
bin with a zero amplitude trace.
In a further preferred embodiment the interpolation step may be modified by
applying a
smoothing scalar s (typically s < 1.0) to the frequency in the sinc function.
In this
embodiment, sinc (x,sf) is used in the weighting function instead of sinc
(x,f), and
similarly sinc (y,sf) is used instead of sinc (y,f). The use of such a scalar
to reduce the
sinc frequency is effective to smooth the data during interpolation.
Figure 9 is a schematic block diagram of an apparatus 14 that is able to
perform a
method according to the present invention.
CA 02546513 2006-05-17
WO 2004/046758 PCT/GB2003/004909
19
The apparatus 14 comprises a programmable data processor 15 with a program
memory
16, for instance in the form of a read only memory (ROM), storing a program
for
controlling the data processor 15 to process seismic data by a method of the
invention.
The apparatus further comprises non-volatile read/write memory 17 for storing,
for
example, any data which must be retained in the absence of a power supply. A
"working" or "scratch pad" memory for the data processor is provided by a
random
access memory RAM 1S. An input device 19 is provided, for instance for
receiving
user commands and data. One or more output devices 20 are provided, for
instance, for
displaying information relating to the progress and result of the processing.
The output
devices) may be, for example, a printer, a visual display unit, or an output
memory.
Sets of seismic data for processing may be supplied via the input device 19 or
may
optionally be provided by a machine-readable data store 21. The results of the
processing may be output via the output device 20 or may be stored.
The program for operating the system and for performing the method described
hereinbefore is stored in the program memory 17, which may be embodied as a
semiconductor memory, for instance of the well known ROM type. However, the
program may well be stored in any other suitable storage medium, such as a
magnetic
data carrier 17a (such as a "floppy disk") or a CD-ROM 17b.
