CA2546513A1 - Processing seismic data - Google Patents
Processing seismic data Download PDFInfo
- Publication number
- CA2546513A1 CA2546513A1 CA002546513A CA2546513A CA2546513A1 CA 2546513 A1 CA2546513 A1 CA 2546513A1 CA 002546513 A CA002546513 A CA 002546513A CA 2546513 A CA2546513 A CA 2546513A CA 2546513 A1 CA2546513 A1 CA 2546513A1
- Authority
- CA
- Canada
- Prior art keywords
- trace
- estimated
- traces
- target location
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Abandoned
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/57—Trace interpolation or extrapolation, e.g. for virtual receiver; Anti-aliasing for missing receivers
Abstract
A method of processing seismic data comprises the step of: obtaining an estimated data trace for a target location by interpolation from a first set of spatially irregularly sampled seismic data traces. The interpolation step comprises applying a weighting function dependent on the target location. The weighting function may be proportional to a sinc function, where sinc (x,f) =
sin (xf)/xf, with a sinc frequency that varies locally dependent on the target location.
sin (xf)/xf, with a sinc frequency that varies locally dependent on the target location.
Description
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.
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.
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 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.
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.
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.
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 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 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.
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.
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:
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 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 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.
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 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 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 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 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.
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.
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 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.
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.
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.
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 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 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.
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.
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:
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 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 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.
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 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 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 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 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.
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.
Claims (30)
1. A method of processing seismic data comprising the step of:
obtaining an estimated data trace for a 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.
obtaining an estimated data trace for a 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.
2. A method as claimed in claim 1 wherein the weighting function is dependent on the value at the target location of the local input density in a co-ordinate direction.
3. A method as claimed in claim 1 or 2 wherein the weighting function is proportional to a sinc function having a locally varying sinc frequency, where sinc (x,f) = sin (fx)/(fx).
4. A method as claimed in claim 3 wherein the weighting function is proportional to a sinc function of the form sinc (x,sf), where s (s < 1) is a scalar.
5. A method as claimed in claim 3 or 4 wherein the sinc frequency is determined from the reciprocal of the gap size at the target location.
6. A method as claimed in any preceding claim wherein the neighbourhood over which the interpolation step is performed is determined adaptively.
7. A method as claimed in any preceding claim wherein a dimension of the neighbourhood over which the interpolation step is performed is subject to a user-defined limit.
8. A method as claimed in any preceding claim and comprising obtaining a plurality of estimated seismic data traces, each for a respective one of a plurality of target locations, from the set of spatially irregularly sampled seismic data traces.
9. A method as claimed in any preceding claim and further comprising the step of estimating at least one attribute for the or each estimated data trace.
10. A method as claimed in claim 9 wherein the at least one attribute for a selected one of the estimated data traces is estimated using the same weighting function as used in the step of estimating the selected estimated data trace.
11. A method as claimed in claim 9 and comprising estimating, for the or each estimated data trace, at least one attribute indicative of the respective interpolation step.
12. A method as claimed in claim 9 or 10 wherein the at least one attribute includes the source-receiver offset.
13. A method as claimed in claim 9, 10 or 11 wherein the at least one attribute includes the source-receiver azimuth.
14. A method as claimed in claim 8, or in any of claims 9 to 13 when dependent directly or indirectly from claim 8, and further comprising obtaining a plurality of estimated seismic data traces, each for one of the plurality of target locations, by interpolation from another set of spatially irregularly sampled seismic data traces.
15. A method as claimed in claim 8, or in any of claims 9 to 14 when dependent directly or indirectly from claim 8, wherein the plurality of target locations define a regular array.
16. A method as claimed in claim 15 wherein the plurality of target locations define a regular two-dimensional array.
17. A method as claimed in any preceding claim wherein the or each target location is a common midpoint co-ordinate.
18. A method as claimed in any preceding claim wherein the or each interpolation step comprises applying:
where D'(x t, y i, t) is the data trace at time t for a target location having spatial co-ordinate (x i,y i), W is a weighting function dependent on the target location, x j, y j are the spatial co-ordinates of the j th data trace, and dip j is the local dip associated with the j th trace.
where D'(x t, y i, t) is the data trace at time t for a target location having spatial co-ordinate (x i,y i), W is a weighting function dependent on the target location, x j, y j are the spatial co-ordinates of the j th data trace, and dip j is the local dip associated with the j th trace.
19. A method as claimed in claim 18 wherein the weighting function is:
20. 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 in any of claims 1 to 19.
21. An apparatus for processing seismic data comprising:
means for obtaining an estimated data trace for a target location by interpolation from a set of spatially irregularly sampled seismic data traces using a weighting function dependent on the target location.
means for obtaining an estimated data trace for a target location by interpolation from a set of spatially irregularly sampled seismic data traces using a weighting function dependent on the target location.
22. An apparatus as claimed in claim 21 and adapted to obtain a plurality of estimated seismic data traces, each for a respective one of a plurality of pre-determined spatial co-ordinates.
23. An apparatus as claimed in claim 21 or 22 and adapted to obtain the or each estimated data trace using a weighting function proportional to a sinc function having a locally varying sinc frequency, where sinc (x,f) = sin (fx)/(fx).
24. An apparatus as claimed in claim 21, 22 or 23 and further comprising means for estimating at least one attribute for the or each estimated data trace.
25. An apparatus as claimed in claim 22 and comprising means for obtaining a plurality of estimated seismic data traces, each for one of the plurality of pre-determined spatial co-ordinates, from a second set of spatially irregularly sampled seismic data traces.
26. An apparatus claimed in any of claims 21 to 25 and adapted to perform the or each interpolation step using:
where D'(x i, y i, t) is the data trace at time t for a target location having spatial co-ordinate (x i,y i), W is a weighting function dependent on the target location, x j, y j are the spatial co-ordinates of the j th data trace, and dip j is the local dip associated with the j th trace.
where D'(x i, y i, t) is the data trace at time t for a target location having spatial co-ordinate (x i,y i), W is a weighting function dependent on the target location, x j, y j are the spatial co-ordinates of the j th data trace, and dip j is the local dip associated with the j th trace.
27. An apparatus claimed in claim 26 and adapted to perform the or each interpolation step using the weighting function:
28. An apparatus as claimed in any of claims 21 to 27 and comprising a programmable data processor.
29. A storage medium comprising a program for a data processor of an apparatus as defined in claim 28.
30. A storage medium containing a program for controlling a data processor to perform a method as defined in any of claims 1 to 20.
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| GB0226899A GB2395559B (en) | 2002-11-19 | 2002-11-19 | Processing seismic data |
| GB0226899.3 | 2002-11-19 | ||
| PCT/GB2003/004909 WO2004046758A1 (en) | 2002-11-19 | 2003-11-13 | Processing seismic data |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CA2546513A1 true CA2546513A1 (en) | 2004-06-03 |
Family
ID=9948065
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CA002546513A Abandoned CA2546513A1 (en) | 2002-11-19 | 2003-11-13 | Processing seismic data |
Country Status (4)
| Country | Link |
|---|---|
| AU (1) | AU2003302127A1 (en) |
| CA (1) | CA2546513A1 (en) |
| GB (1) | GB2395559B (en) |
| WO (1) | WO2004046758A1 (en) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2444953B (en) * | 2006-12-19 | 2009-07-22 | Westerngeco Seismic Holdings | Method for obtaining an image of a subsurface by regularizing irregularly sampled seismic data |
| GB2450122B (en) | 2007-06-13 | 2009-08-05 | Westerngeco Seismic Holdings | Method of representing signals |
| US8339898B2 (en) * | 2008-05-25 | 2012-12-25 | Westerngeco L.L.C. | Processing seismic data using combined regularization and 4D binning |
| GB2501934B (en) * | 2012-05-11 | 2014-07-09 | Statoil Petroleum As | Method and apparatus for processing seismic data |
| CN105277974A (en) * | 2014-07-23 | 2016-01-27 | 中国石油化工股份有限公司 | Stratum data interpolation method |
| US10408958B2 (en) * | 2015-07-22 | 2019-09-10 | Fairfield Industries Incorporated | Systems and methods for detecting subsurface features using 3D angle gathers |
| US10393901B2 (en) | 2015-08-31 | 2019-08-27 | Pgs Geophysical As | Wavefield interpolation and regularization in imaging of multiple reflection energy |
| CN108919352B (en) * | 2018-05-17 | 2019-11-12 | 中国海洋石油集团有限公司 | A kind of fidelity imaging method based on scalar wave wave field extrapolation |
Family Cites Families (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4887244A (en) * | 1988-06-28 | 1989-12-12 | Mobil Oil Corporation | Method for seismic trace interpolation using a forward and backward application of wave equation datuming |
| US4964098A (en) * | 1990-03-15 | 1990-10-16 | Exxon Production Research Company | Method for seismic trace interpolation |
| FR2675590B1 (en) * | 1991-04-19 | 1993-07-16 | Elf Aquitaine | PROCESSING PROCESS FOR OBTAINING A SUM SECTION WITH ZERO OFFSETS. |
| GB2256048B (en) * | 1991-05-23 | 1994-08-31 | Geco As | Method of seismic processing |
| US5319554A (en) * | 1991-08-30 | 1994-06-07 | Shell Oil Company | Method for data interpolation and signal enhancement |
| US5235556A (en) * | 1992-01-10 | 1993-08-10 | Halliburton Geophysical Services Inc. | Interpolation of aliased seismic traces |
| GB9320540D0 (en) * | 1993-10-06 | 1993-11-24 | Ensign Geophysics Ltd | Seismic data acquisition |
| US5677892A (en) * | 1996-08-14 | 1997-10-14 | Western Atlas International, Inc. | Unaliased spatial trace interpolation in the f-k domain |
| US6021379A (en) * | 1997-07-29 | 2000-02-01 | Exxon Production Research Company | Method for reconstructing seismic wavefields |
-
2002
- 2002-11-19 GB GB0226899A patent/GB2395559B/en not_active Expired - Fee Related
-
2003
- 2003-11-13 AU AU2003302127A patent/AU2003302127A1/en not_active Abandoned
- 2003-11-13 WO PCT/GB2003/004909 patent/WO2004046758A1/en not_active Application Discontinuation
- 2003-11-13 CA CA002546513A patent/CA2546513A1/en not_active Abandoned
Also Published As
| Publication number | Publication date |
|---|---|
| GB2395559B (en) | 2006-01-18 |
| GB2395559A (en) | 2004-05-26 |
| WO2004046758A1 (en) | 2004-06-03 |
| AU2003302127A1 (en) | 2004-06-15 |
| GB0226899D0 (en) | 2002-12-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP1779144B1 (en) | Offeset-azimuth binning for migration and velocity analysis | |
| AU2009215000B2 (en) | Method for determining geological information related to a subsurface volume of interest | |
| US6735527B1 (en) | 3-D prestack/poststack multiple prediction | |
| CA2712618C (en) | Method for indexing a subsurface volume for the purpose of inferring geologic information | |
| CA2558000C (en) | Method, media, and signals for processing seismic data to obtain a velocity field | |
| US8082106B2 (en) | 3D surface related multiple elimination for wide azimuth seismic data | |
| US6487502B1 (en) | System for estimating the locations of shaley subsurface formations | |
| US20050114032A1 (en) | Seismic data interpolation system | |
| EP2946232B1 (en) | Wavefield modelling and 4d-binning for seismic surveys with different acquisition datums | |
| US6493634B1 (en) | Method for determining stacking velocity parameters or other reflection geometry information from seismic gather data using multiple attributes and 3-D visualization | |
| GB2420408A (en) | Method for processing at least two sets of seismic data | |
| CA2546513A1 (en) | Processing seismic data | |
| US6625543B1 (en) | Output based azimuth moveout re-gridding of seismic data | |
| EP1348138B1 (en) | A seismic method and apparatus for generating a semblance panel and computing the reflector dip | |
| US5629905A (en) | Method of datuming seismic data and method of processing seismic data | |
| Schonewille et al. | A general reconstruction scheme for dominant azimuth 3D seismic data | |
| EP2743737A2 (en) | Methods and systems for quality control of seismic illumination maps | |
| AU2002216298A1 (en) | A seismic method and apparatus for generating a semblance panel and computing the reflector dip | |
| WO1998015849A1 (en) | Method of generating a fold distribution and of evaluating a seismic survey | |
| US20220413175A1 (en) | Enhanced projection on convex sets for interpolation and deblending | |
| Williams et al. | Moving towards full-sampling in land 3D acquisition | |
| AU2015201019A1 (en) | Method for indexing a subsurface volume for the purpose of inferring geologic information |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| EEER | Examination request | ||
| FZDC | Discontinued application reinstated | ||
| FZDE | Discontinued |
Effective date: 20121113 |