CA2533012C - Method of estimating the illumination fold in the migrated domain - Google Patents
Method of estimating the illumination fold in the migrated domain Download PDFInfo
- Publication number
- CA2533012C CA2533012C CA2533012A CA2533012A CA2533012C CA 2533012 C CA2533012 C CA 2533012C CA 2533012 A CA2533012 A CA 2533012A CA 2533012 A CA2533012 A CA 2533012A CA 2533012 C CA2533012 C CA 2533012C
- Authority
- CA
- Canada
- Prior art keywords
- seismic
- travel time
- illumination
- source
- fold
- 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.)
- Expired - Fee Related
Links
- 238000005286 illumination Methods 0.000 title claims abstract description 77
- 238000000034 method Methods 0.000 title claims description 52
- 239000013598 vector Substances 0.000 claims abstract description 22
- 108050005509 3D domains Proteins 0.000 claims abstract description 10
- 238000013508 migration Methods 0.000 claims description 49
- 230000005012 migration Effects 0.000 claims description 48
- 239000002689 soil Substances 0.000 claims description 22
- 238000011161 development Methods 0.000 claims description 13
- 238000012937 correction Methods 0.000 claims description 9
- 230000014509 gene expression Effects 0.000 claims description 7
- 230000000694 effects Effects 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 5
- 230000018109 developmental process Effects 0.000 description 12
- 238000002310 reflectometry Methods 0.000 description 10
- 230000005526 G1 to G0 transition Effects 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 238000013507 mapping Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000005303 weighing Methods 0.000 description 1
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. for interpretation or for event detection
-
- 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. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Electroluminescent Light Sources (AREA)
Abstract
The invention relates to an estimate of the seismic illumination fold (x,p) in the migrated 3D domain at an image point x, for a dip of vector p characterised in that the illumination fold I (z, p; s, r) is estimated for each (source s, receiver r) pair in the seismic survey, by applying the following steps: - determination of the reflection travel time tr (.xr (p) ; s , r)from the source s to the specular reflection point z, on the plane reflector passing through the image point x and perpendicular to the dip vector p, and then return to the reflector r; starting from the diffraction travel time td(.z;s, r ) from the source to the said image point x and then return to the reflector r; - incrementing the said illumination fold I (X, p;s, r) related to the said (source s, receiver r) pair as a function of the difference between the diffraction travel time td ( x; s, r) and the reflection travel time tr(xr (p)isrr).
Description
METHOD OF ESTIMATING THE ILLUMINATION FOLD IN THE
MIGRATED DOMAIN
The field of the invention is seismic prospecting.
The invention relates to seismic imagery and more particularly to estimating the illumination fold in the migrated 3D domain necessary to obtain ~ high resolution image of the underground structure at true amplitude.
Seismic prospecting generally consists of emitting seismic waves into the sub-soil using one or several seismic sources, and recording seismic data corresponding to seismic waves reflected on geological interfaces in the sub-soil (also called reflectors) on the ground surface,_ as a function of time, using receivers (also called geophones or hydrophones depending on whether the survey is being made on land or at sea) and then processing these data to extract useful information about the geology of the sub-soil.
Data processing (or seismic traces) recorded by receivers are useful particularly for designing images with information about the geology of the sub-soil.
The objective of high resolution seismic imagery with true amplitude is to provide an image of the sub-soil that is not disturbed, particularly by wave field propagation effects in the sub-soil and by acquisition effects (we then talk about acquisition patterns).
Acoustic reflectivity of the sub-soil is represented as an image using a so-called migration method to restore correct forms of geological interfaces. To achieve this, seismic data are converted from the domain in which they are acquired to the migrated domain in which geological interfaces are represented in their genuine locations.
In the migrated domain, the vertical scale may be time (we then talk about time migration) or depth (we then talk about depth migration).
Migration is usually done before summation (we then talk about pre-stack migration), with direct mapping of the seismic signal from seismic traces towards the migrated domain.
The precise identification of the sub-soil geology at a particular point is directly related to the capability of a seismic emission to illuminate this point. If a point in the sub-soil is highly illuminated, in all likelihood the seismic survey can evaluate the geology of the sub-soil at this point. On the other hand, if a point of the sub-soil is not well illuminated, it is quite possible that the seismic survey will provide an imprecise or even .incorrect evaluation of the sub-soil geology at this point.
In order to obtain a precise representation of the sub-soil, the true amplitudes of the reflection of the seismic wave on a reflector must be represented perfectly.
Therefore, it is necessary to know precisely if a change in the amplitude at a specific point in the sub-soil is the result of a change in the geology of the sub-soil at this point or if it is the result of a.particular illumination at this point.
But reflectors are usually illuminated non-uniformly. In particular, the illumination depends on:
- the configuration of seismic data acquisition;
- lateral and vertical variations of the propagation velocity of seismic waves in the sub-soil.
In order to obtain a high resolution image of the reflectivity of the sub-soil, it is necessary to make an amplitude correction to take account of non-uniformity of illumination of reflectors when migrating seismic reflections.
As has already been described above, this correction (done using a weighting factor) is necessary because:
- the reflectivity of the sub-soil cannot be precisely imaged without correcting amplitude distortions related to propagation of the wave field in the sub-soil (spherical divergence, absorption, directivity, etc.);
- to improve the resolution of the reflectivity image, it is important to take account of the variability of the migration illumination fold for each inclination of a reflector used during the migration.
For a given image point in the migrated domain, the illumination fold corresponds to the number of source receiver pairs for which contributions for building up the amplitude of the reflectivity image interfere constructively (stationary phase condition).
The document "Bleistein, N. 1987, On the Imaging of Reflectors in the Earth: Geophysics Soc. of Expl.
Geophysics, 931-942" proposes a method used to obtain an amplitude weighting factor for a regular acquisition configuration, capable of making a high resolution representation of the reflectivity of the sub-soil in true amplitude.
This is an "inverse" method composed of a modelling step (modelling of seismic traces from a reflectivity model) followed by a migration step (from modelled seismic traces to imaged reflectivity), the weighting factor then being chosen so as to obtain perfect agreement between the imaged reflectivity and the reflectivity model.
By following this "inverse" method proposed by Bleistein and as shown in Figure 1, it can be seen that:
- for a given acquisition configuration (di.stribution of (source S, receiver R) pairs), - for a given image point x in the migrated domain, - for a given reflector inclination (indicated by a dip vector p), the illumination fold ~(x,p) is determined by the number of source-receiver pairs such that their migration contributions constructively interfere with the image of the dip reflector p at image point x.
Expressions developed by Bleistein for the illumination fold ~ ( .z,p ) include the effect of the entire geometric configuration of the acquisition in a single step.
Consequently, the method proposed by Bleistein can only be used for a regular and parametric distribution of all source-receiver pairs.
Consequently, this restrictive application condition limits its use to a small number of acquisition configurations.
Therefore for acquisition configurations with irregular geometries, it is necessary to define individual illumination folds ~ ( x,p,S,Y ) , for each (source s, receiver s°) pair involved in the migration.
MIGRATED DOMAIN
The field of the invention is seismic prospecting.
The invention relates to seismic imagery and more particularly to estimating the illumination fold in the migrated 3D domain necessary to obtain ~ high resolution image of the underground structure at true amplitude.
Seismic prospecting generally consists of emitting seismic waves into the sub-soil using one or several seismic sources, and recording seismic data corresponding to seismic waves reflected on geological interfaces in the sub-soil (also called reflectors) on the ground surface,_ as a function of time, using receivers (also called geophones or hydrophones depending on whether the survey is being made on land or at sea) and then processing these data to extract useful information about the geology of the sub-soil.
Data processing (or seismic traces) recorded by receivers are useful particularly for designing images with information about the geology of the sub-soil.
The objective of high resolution seismic imagery with true amplitude is to provide an image of the sub-soil that is not disturbed, particularly by wave field propagation effects in the sub-soil and by acquisition effects (we then talk about acquisition patterns).
Acoustic reflectivity of the sub-soil is represented as an image using a so-called migration method to restore correct forms of geological interfaces. To achieve this, seismic data are converted from the domain in which they are acquired to the migrated domain in which geological interfaces are represented in their genuine locations.
In the migrated domain, the vertical scale may be time (we then talk about time migration) or depth (we then talk about depth migration).
Migration is usually done before summation (we then talk about pre-stack migration), with direct mapping of the seismic signal from seismic traces towards the migrated domain.
The precise identification of the sub-soil geology at a particular point is directly related to the capability of a seismic emission to illuminate this point. If a point in the sub-soil is highly illuminated, in all likelihood the seismic survey can evaluate the geology of the sub-soil at this point. On the other hand, if a point of the sub-soil is not well illuminated, it is quite possible that the seismic survey will provide an imprecise or even .incorrect evaluation of the sub-soil geology at this point.
In order to obtain a precise representation of the sub-soil, the true amplitudes of the reflection of the seismic wave on a reflector must be represented perfectly.
Therefore, it is necessary to know precisely if a change in the amplitude at a specific point in the sub-soil is the result of a change in the geology of the sub-soil at this point or if it is the result of a.particular illumination at this point.
But reflectors are usually illuminated non-uniformly. In particular, the illumination depends on:
- the configuration of seismic data acquisition;
- lateral and vertical variations of the propagation velocity of seismic waves in the sub-soil.
In order to obtain a high resolution image of the reflectivity of the sub-soil, it is necessary to make an amplitude correction to take account of non-uniformity of illumination of reflectors when migrating seismic reflections.
As has already been described above, this correction (done using a weighting factor) is necessary because:
- the reflectivity of the sub-soil cannot be precisely imaged without correcting amplitude distortions related to propagation of the wave field in the sub-soil (spherical divergence, absorption, directivity, etc.);
- to improve the resolution of the reflectivity image, it is important to take account of the variability of the migration illumination fold for each inclination of a reflector used during the migration.
For a given image point in the migrated domain, the illumination fold corresponds to the number of source receiver pairs for which contributions for building up the amplitude of the reflectivity image interfere constructively (stationary phase condition).
The document "Bleistein, N. 1987, On the Imaging of Reflectors in the Earth: Geophysics Soc. of Expl.
Geophysics, 931-942" proposes a method used to obtain an amplitude weighting factor for a regular acquisition configuration, capable of making a high resolution representation of the reflectivity of the sub-soil in true amplitude.
This is an "inverse" method composed of a modelling step (modelling of seismic traces from a reflectivity model) followed by a migration step (from modelled seismic traces to imaged reflectivity), the weighting factor then being chosen so as to obtain perfect agreement between the imaged reflectivity and the reflectivity model.
By following this "inverse" method proposed by Bleistein and as shown in Figure 1, it can be seen that:
- for a given acquisition configuration (di.stribution of (source S, receiver R) pairs), - for a given image point x in the migrated domain, - for a given reflector inclination (indicated by a dip vector p), the illumination fold ~(x,p) is determined by the number of source-receiver pairs such that their migration contributions constructively interfere with the image of the dip reflector p at image point x.
Expressions developed by Bleistein for the illumination fold ~ ( .z,p ) include the effect of the entire geometric configuration of the acquisition in a single step.
Consequently, the method proposed by Bleistein can only be used for a regular and parametric distribution of all source-receiver pairs.
Consequently, this restrictive application condition limits its use to a small number of acquisition configurations.
Therefore for acquisition configurations with irregular geometries, it is necessary to define individual illumination folds ~ ( x,p,S,Y ) , for each (source s, receiver s°) pair involved in the migration.
5 Obviously, the said individual ratios contribute to the total illumination fold ~(x,p) according to I(x,p) -~I ( x~h~s~~" ) ) .
s,a.
However, the study carried out by Bleistein is based on a number of developments that were intended to estimate the illumination fold in the migrated domain for any geometric acquisition configuration, and for any distribution of the source-receiver pairs.
For years, various approximations have been proposed for the seismic wave that is the key element for estimating the stationary phase condition (constructive interference condition).
The first developments consisted of modelling specular reflection points along a horizon (~D map) selected in the migrated domain (3D) by ray tracing for 30 each of the (source s~ receiver i°) pairs.
This reflection modelling step determines the following quantities:
- location of reflection points .~,.(s,y~) along the selected horizon;
- the specular reflection travel time tr (~,.as,f~);
The diffraction travel time td along the selected horizon can then be estimated by a development of the reflection travel time tr, around the specular reflection point .z,.(s,f~) in the Taylor series (determination of first and second derivatives ~~a,.(a;,.) and ~z,~t,.(x,.) ) .
s,a.
However, the study carried out by Bleistein is based on a number of developments that were intended to estimate the illumination fold in the migrated domain for any geometric acquisition configuration, and for any distribution of the source-receiver pairs.
For years, various approximations have been proposed for the seismic wave that is the key element for estimating the stationary phase condition (constructive interference condition).
The first developments consisted of modelling specular reflection points along a horizon (~D map) selected in the migrated domain (3D) by ray tracing for 30 each of the (source s~ receiver i°) pairs.
This reflection modelling step determines the following quantities:
- location of reflection points .~,.(s,y~) along the selected horizon;
- the specular reflection travel time tr (~,.as,f~);
The diffraction travel time td along the selected horizon can then be estimated by a development of the reflection travel time tr, around the specular reflection point .z,.(s,f~) in the Taylor series (determination of first and second derivatives ~~a,.(a;,.) and ~z,~t,.(x,.) ) .
Figure 2 illustrates this procedure during which the diffraction travel time td(x) at image point x is estimated using the reflection travel time tr(:x,.) and more particularly by its development around the specular reflection point x,, in the Taylor series.
Summation of all wavelet mappings obtained for each source-receiver pair then provides a means of estimating the illumination fold in the migrated domain along the ' selected horizon.
, By modelling the seismic wavelet according to a pulse peak, with infinite frequency band, document "Cain, G., Cambois, G., Gehin, M and Hall, R., 1998 Reducing risk in seismic acquisition and interpretation of complex targets using a Gocad-based 3-D modelling tool, 68th Ann.
Internet. Mtg, Soc. of Expl. Geophys., 2072-2075" has shown that an impact map can be obtained, and in this case only the specular contribution on a selected horizon is actually taken into account.
Considering a realistic approximation of a seismic wavelet with limited band, the following documents:
- "Schneider, W., Jr and Winbow, G., 1999, Efficient and accurate modelling of 3D seismic illumination, 69th Ann. Internet. Mtg, Soc. Of Expl. Geophys, Expanded Abstract, 633-636"a - "Laurain, R. and Vinje, V. 2001, Pre-stack migration and illumination maps, 7lth Ann. Internet., Mtg: Soc. Of Expl. Geophys, Expanded Abstract, 929-932"~
have shown that it is possible to improve estimating the illumination fold along selected horizons.
These approaches, in particular based on scanning of selected horizons, have been found to be reliable and precise. However, their efficiency is limited particularly due to the fact that:
- compensation for the illumination fold is limited to a finite number of selected horizons and therefore, unlike the method developed by Bleistein, this method is not completely 3D;
- for each (source s, receiver f°) pair involved in the migration, a dedicated modelling step must be carried out after the migration. However, this modelling step to determine the (x,.,~xt,.,~x,xt,.) quantities mentioned previously by ray tracing along each of the selected horizons is particularly expensive in terms of time.
In order to overcome these limitations, and in order to. approximate a 3D approach, a two-step solution has been proposed to estimate and then to compensate.for the illumination ratio, and particularly in documents:
- "Albertin, U., Bloor, R., Beasley, C., Chang, W., Jaramillo, H., Mobley, E. and Yingst, D., 1999, Aspects of true amplitude migration, 69th Ann. Internat. Mtg:
Soc. Of Expl. Geophys, Expanded Abstract, 1247-1250";
- "Audebert, F., Froidevaux, P. Huard, I., Nicoletis, L. and Svay-Lucas, J., 2000, A multi angle toolbox for restored amplitude images and AVA gathers, 70th Ann. Internat, Mtg: Soc. Of Expl. Geophys, Expanded Abstract, 1020-1023".
The first step consists of making an analysis of the multi-dip illumination fold using predetermined tables for Kirchhoff migration for diffraction travel times td ( ~- ; s , ~ ) , in the migrated domain .
In the second step, the inverse values of the estimated illumination folds are used to make the weighted summation of seismic data.
Summation of all wavelet mappings obtained for each source-receiver pair then provides a means of estimating the illumination fold in the migrated domain along the ' selected horizon.
, By modelling the seismic wavelet according to a pulse peak, with infinite frequency band, document "Cain, G., Cambois, G., Gehin, M and Hall, R., 1998 Reducing risk in seismic acquisition and interpretation of complex targets using a Gocad-based 3-D modelling tool, 68th Ann.
Internet. Mtg, Soc. of Expl. Geophys., 2072-2075" has shown that an impact map can be obtained, and in this case only the specular contribution on a selected horizon is actually taken into account.
Considering a realistic approximation of a seismic wavelet with limited band, the following documents:
- "Schneider, W., Jr and Winbow, G., 1999, Efficient and accurate modelling of 3D seismic illumination, 69th Ann. Internet. Mtg, Soc. Of Expl. Geophys, Expanded Abstract, 633-636"a - "Laurain, R. and Vinje, V. 2001, Pre-stack migration and illumination maps, 7lth Ann. Internet., Mtg: Soc. Of Expl. Geophys, Expanded Abstract, 929-932"~
have shown that it is possible to improve estimating the illumination fold along selected horizons.
These approaches, in particular based on scanning of selected horizons, have been found to be reliable and precise. However, their efficiency is limited particularly due to the fact that:
- compensation for the illumination fold is limited to a finite number of selected horizons and therefore, unlike the method developed by Bleistein, this method is not completely 3D;
- for each (source s, receiver f°) pair involved in the migration, a dedicated modelling step must be carried out after the migration. However, this modelling step to determine the (x,.,~xt,.,~x,xt,.) quantities mentioned previously by ray tracing along each of the selected horizons is particularly expensive in terms of time.
In order to overcome these limitations, and in order to. approximate a 3D approach, a two-step solution has been proposed to estimate and then to compensate.for the illumination ratio, and particularly in documents:
- "Albertin, U., Bloor, R., Beasley, C., Chang, W., Jaramillo, H., Mobley, E. and Yingst, D., 1999, Aspects of true amplitude migration, 69th Ann. Internat. Mtg:
Soc. Of Expl. Geophys, Expanded Abstract, 1247-1250";
- "Audebert, F., Froidevaux, P. Huard, I., Nicoletis, L. and Svay-Lucas, J., 2000, A multi angle toolbox for restored amplitude images and AVA gathers, 70th Ann. Internat, Mtg: Soc. Of Expl. Geophys, Expanded Abstract, 1020-1023".
The first step consists of making an analysis of the multi-dip illumination fold using predetermined tables for Kirchhoff migration for diffraction travel times td ( ~- ; s , ~ ) , in the migrated domain .
In the second step, the inverse values of the estimated illumination folds are used to make the weighted summation of seismic data.
For each image point x and for an entire reflector dip range p, it has been proposed to increment an illumination occurrence counter I (.z,p ) every time that a migration contribution reaches the image point x with the dip component p.
In a manner similar to that proposed by Cain et al (see above), the criterion used to increase the illumination occurrence counter takes account of a pulse model of the seismic wavelet, with infinite frequency band.
This high frequency approximation is such that for a given dip component, only specular migration contributions corresponding to the said dip component are used for estimating the illumination.
And there are a number of limitations to this high frequency approximation:
- the estimated illumination is biased since the limited band nature of the seismic signal is ignored;
- in order to avoid division by zeros, the illumination folds are stored by dip classes. However, the size of each class is a critical element. If the class size is too small, risks of numerical instability are very high; if the class size is too large, the precision of the estimated illumination coverage is low.
The invention is designed to provide a method for estimating the seismic illumination fold in the migrated 3D domain that is not subject to these limitations.
More precisely, the purpose of the invention is to propose a method that takes account of the limited band nature of the seismic signal in order to estimate more realistic illumination folds, and which at the same time solves the problem of the class size.
In a manner similar to that proposed by Cain et al (see above), the criterion used to increase the illumination occurrence counter takes account of a pulse model of the seismic wavelet, with infinite frequency band.
This high frequency approximation is such that for a given dip component, only specular migration contributions corresponding to the said dip component are used for estimating the illumination.
And there are a number of limitations to this high frequency approximation:
- the estimated illumination is biased since the limited band nature of the seismic signal is ignored;
- in order to avoid division by zeros, the illumination folds are stored by dip classes. However, the size of each class is a critical element. If the class size is too small, risks of numerical instability are very high; if the class size is too large, the precision of the estimated illumination coverage is low.
The invention is designed to provide a method for estimating the seismic illumination fold in the migrated 3D domain that is not subject to these limitations.
More precisely, the purpose of the invention is to propose a method that takes account of the limited band nature of the seismic signal in order to estimate more realistic illumination folds, and which at the same time solves the problem of the class size.
Consequently, according to a first aspect of the invention, it divulges a method for estimating the seismic illumination fold ~(x,p) in the migrated 3D
domain at an image point x for a dip of vector p characterised in that it estimates the illumination fold I (x, p,s, y-) for each ( source s , receiver s° ) pair in the seismic survey, by applying the following steps:
- determination of the reflection travel time tr ( x,. ( p ) ; s f" ) from the source s to the specular reflection point :x,. on the plane reflector passing through the image point x and perpendicular to the dip vector p and then return to the reflector starting from the diffraction travel time td(x;s,~) from the source s to the said image point .~ and then return to the reflector ~';
- incrementing the said illumination fold I (x, p,s,f"") related to the said (source s, receiver f") pair as a function of the difference between the diffraction travel time td ( x;s,~ ) and the reflection travel time tr(.x,. (p)isj") Preferred but non-limitative aspects of the method according to the first aspect of the invention are as follows:
- the method also comprises a step to summate each illumination fold I(.~,p,s,i°) related to a (source s, receiver i"") pair so as to determine the total illumination fold I (x, p) - ~I(x, p, s, s") , s,p - during the increment step, the illumination fold I (x, p, s,i-) may be incremented using an increment function 10, i ( td, tr i s, j' ) using the equation I (x, p, s, i°°) - I (x, p, s, ~) +
1 (tdr tri s~j~ ) r - the increment function i may be a function of the seismic wavelet s (t) , - the increment function i may be expressed as a function of the derivative of the seismic wavelet with respect to time, - it is envisaged that an a priori correction w ( x,s,~ ) of the illumination fold could be taken into account, which usually has to be considered during migration, - the determination step includes the second order Taylor series development of the diffraction travel time td ( x;s,f~ ) around the image point x , - the specular reflection point zr(p) is determined along the length of the said reflector such that the diffraction travel time at the said specular reflection point x,. ( p ) is stationary, , - the determination step uses isochronic migration maps td ( z;s,~ ) specified for each (source s , receiver pair involved in the migration and each image point x in the migrated 3D domain, - the seismic illumination fold I(.~,p) in the migrated 3D domain is estimated during the Kirchhoff summation migration of seismic data recorded during the 3D seismic prospecting.
According to another aspect, the invention proposes a method for correction of seismic data amplitudes recorded during a 3D seismic survey in order to compensate for the effect of non-uniform illumination of sub-soil reflectors, characterised in that it comprises steps consisting of:
- estimating the illumination fold I(x,p) using the method according to the first aspect of the invention, - using the inverse I-1 ( x,p ) of the said ratio as a weighting factor to be applied to each of the said seismic data amplitudes.
According to another aspect, the invention proposes a method for selection of an acquisition geometry among a plurality of acquisition geometries as a function of the target of 3D seismic prospecting, characterised in that it consists of the following steps:
- determining the illumination fold I(x,p) by the method according to the first aspect of the invention, for each of the acquisition geometries considered, - selecting the acquisition geometry providing the optimum illumination fold as a function of the target.
Other aspects, purposes and advantages of the invention will become clearer after reading the following detailed description of a preferred embodiment of it, given as a non-limitative example with reference to the attached drawings on which:
- Figure 1, already commented, illustrates the constructive or non-constructive contribution as a function of s(t~-tr), of the seismic trace migrated to the image of the reflector with. a dip along .~~
- Figure 2, already commented, illustrates the estimate made according to the state of the art for the diffraction travel time at the image point td(x) starting from a Taylor series development of the reflection travel time at the specular reflection point tr(xr);
domain at an image point x for a dip of vector p characterised in that it estimates the illumination fold I (x, p,s, y-) for each ( source s , receiver s° ) pair in the seismic survey, by applying the following steps:
- determination of the reflection travel time tr ( x,. ( p ) ; s f" ) from the source s to the specular reflection point :x,. on the plane reflector passing through the image point x and perpendicular to the dip vector p and then return to the reflector starting from the diffraction travel time td(x;s,~) from the source s to the said image point .~ and then return to the reflector ~';
- incrementing the said illumination fold I (x, p,s,f"") related to the said (source s, receiver f") pair as a function of the difference between the diffraction travel time td ( x;s,~ ) and the reflection travel time tr(.x,. (p)isj") Preferred but non-limitative aspects of the method according to the first aspect of the invention are as follows:
- the method also comprises a step to summate each illumination fold I(.~,p,s,i°) related to a (source s, receiver i"") pair so as to determine the total illumination fold I (x, p) - ~I(x, p, s, s") , s,p - during the increment step, the illumination fold I (x, p, s,i-) may be incremented using an increment function 10, i ( td, tr i s, j' ) using the equation I (x, p, s, i°°) - I (x, p, s, ~) +
1 (tdr tri s~j~ ) r - the increment function i may be a function of the seismic wavelet s (t) , - the increment function i may be expressed as a function of the derivative of the seismic wavelet with respect to time, - it is envisaged that an a priori correction w ( x,s,~ ) of the illumination fold could be taken into account, which usually has to be considered during migration, - the determination step includes the second order Taylor series development of the diffraction travel time td ( x;s,f~ ) around the image point x , - the specular reflection point zr(p) is determined along the length of the said reflector such that the diffraction travel time at the said specular reflection point x,. ( p ) is stationary, , - the determination step uses isochronic migration maps td ( z;s,~ ) specified for each (source s , receiver pair involved in the migration and each image point x in the migrated 3D domain, - the seismic illumination fold I(.~,p) in the migrated 3D domain is estimated during the Kirchhoff summation migration of seismic data recorded during the 3D seismic prospecting.
According to another aspect, the invention proposes a method for correction of seismic data amplitudes recorded during a 3D seismic survey in order to compensate for the effect of non-uniform illumination of sub-soil reflectors, characterised in that it comprises steps consisting of:
- estimating the illumination fold I(x,p) using the method according to the first aspect of the invention, - using the inverse I-1 ( x,p ) of the said ratio as a weighting factor to be applied to each of the said seismic data amplitudes.
According to another aspect, the invention proposes a method for selection of an acquisition geometry among a plurality of acquisition geometries as a function of the target of 3D seismic prospecting, characterised in that it consists of the following steps:
- determining the illumination fold I(x,p) by the method according to the first aspect of the invention, for each of the acquisition geometries considered, - selecting the acquisition geometry providing the optimum illumination fold as a function of the target.
Other aspects, purposes and advantages of the invention will become clearer after reading the following detailed description of a preferred embodiment of it, given as a non-limitative example with reference to the attached drawings on which:
- Figure 1, already commented, illustrates the constructive or non-constructive contribution as a function of s(t~-tr), of the seismic trace migrated to the image of the reflector with. a dip along .~~
- Figure 2, already commented, illustrates the estimate made according to the state of the art for the diffraction travel time at the image point td(x) starting from a Taylor series development of the reflection travel time at the specular reflection point tr(xr);
- Figure 3 illustrates the method according to the invention according to which the specular reflection point x is estimated, together with the reflection travel time tr(xr) starting from a Taylor series development of the diffraction travel time, at the image point td(x).
- Figure 4 is a flowchart diagrammatically showing the main steps in the method according to the invention.
The invention proposes a method of estimating the illumination fold in the migrated domain I ( x , p ; s , F ) for any migrated position x, regardless of its dip component p and regardless of the distribution of seismic traces, using only the most basic quantities involved in migration by a Kirchhoff summation, namely:
- the coordinates (s,:r) of source-receiver pairs, where s represents the position of a source S and z' represents the position of a receiver R;*
- migration isochronic maps td(x;s,r) for each of the (source s, receiver r) pairs involved in the migration.
The invention is equally applicable to imagery in time and imagery in depth.
The invention divulges a method of the type including a modelling step followed by a migration step.
Concerning the modelling step of the reflection point ~ and the reflection travel time tr(.~r), the method according to the invention applies the stationary phase condition to a Taylor series development of diffraction travel times td(x).
The said development is a second order development.
Concerning the migration step, the method according to the invention uses isochronic migration maps available for the migration conventionally done by a Kirchhoff summation .
For each of the (source s, receiver r) pairs, the illumination fold I ( x, p ; s , r ) is updated with reference to Figure 3 using the following equation:
I(x,pis,r) - I(x,pis,r) + I(td,tr; s,r) Equation (1) where is the image point in the migrated domain x -(xlr x2r x3)Ti - p is the dip vector representing the selected dip component, p - (pi r pa r p3 ) T
- s defines the coordinates of the source, s - (s1, s2r S3)Ti - r defines the coordinates of the receiver, r -T, 1 5 (rlr r2r r3) r - td(x; s, r) is the diffraction travel time from the source s, as far as the migrated position x, and return as far as the receiver.r (path indicated by td on Figure 3 and shown in solid lines);
- tr(xr(p)i s, r) is the reflection travel time from the source s as far as point x of specular reflection on the plane reflector passing through image point ~: and perpendicular to the dip vector Via, and return as far as the receiver ~' (path indicated by tr on Figure 3 and shown in dashed lines);
- I(td,tri s, ~') is an increment function for which the expression will be described in detail below.
As already mentioned above, the method according to the invention advantageously uses isochronic migration maps td(x; s, r) involved in the migration by Kirchhoff summation, for each of the ( source s , receiver ~ ) pairs involved in the migration. These maps are "input data" to the method according to the invention.
It also uses first and second derivatives O Xtd ( x ; s , r ) , ~X,Xta ( x ; S . r ) of migration isochrones when they are available. If not, the said derivatives are not input data, but may be calculated as part of the method according to the invention.
Figure 4 is a flow chart showing the different steps in the method according to a preferred embodiment of the invention.
The said method consists of carrying out a step 1 used for each image point x in the migrated 3D domain, for each of the (source s, receiver r) pairs involved in the migration.
The said step 1 includes two operations 2, 3, for each dip component p selected for the illumination calculation I ( x ; p ; s , r ) .
The first operation 2 consists of determining the reflection travel time tr ( xr ( p ) ; s , r ) to the specular reflection point xr(p) on -the plane reflector passing through the image point x and perpendicular to the dip vector p .
~peration 2 is determined from the diffraction travel time td(x; s, r) known from the isochronic migration maps.
The second operation 3 consists of incrementing the illumination fold in the migrated domain I(x,p;s,r) of an increment function I(td, tr; s, r), the said increment function taking account of:
- the difference {td ( x ; s , r ) - tr ( xr ( p ) ; s , r ) ~
between the travel times mentioned above-;
- the seismic wavelet s(t);
- and possibly the migration amplitude term w(x;s.r).
At the end of operations 2 and 3, the illumination 5 amplitude ratio I ( x, p ; s , r ) is known for a given image point x and a dip component p. .
By repeating the operations mentioned above for each of the (source s, receiver r) pairs, the final step is to determine the illumination ratio in the migrated 10 domain I ( x , p ) .
This is done by adding therelative illumination to each pair ( source s , receiver r ) according to I ( x , p ) -~I (x.h%s.-r').
s,;.
The reflection point xr and the reflection travel l5 time tr are determined (during operation 2 mentioned above) starting from a stationary phase analysis, by a second order Taylor series development of the diffraction travel time td around the image point x .
The following expressions for xr and tr (equations (3) and (4) ) are obtained from this analysis (see equation (2) below):
td(~;~,.~'>=tr((~;~s.r>+(~xtr>r(.~-.~r>+2 (.~-%r)T(~~..utr)(x-.~r)+...
Equation (2) .~r( p ) - ~-M.F-1. b Equation (3) tr(xr(~);s.r) - td(x;S.r) y.j7r. F1. j7 Equation (4) where, in equations (3) and (4).
- M is a (3 x 2) matrix described by two vectors extending along the length of the reflection plane.
Consequently, these two vectors are perpendicular to the dip vector p.
- b is a (2 x 1) vector of the second order derivatives of the diffraction travel time along the reflection plane, b - MT . ( 0 Xtd ) Equation ( 5 ) - F is a (2 x 3) matrix of first order derivatives of the diffraction travel time along the reflection plane, F =_ MT- (0 X~xtd) ~M Equation (6) The following describes how equations (3) and (4) are obtained.
Note that the dip vector p represents the apparent dip of a reflector in the migrated domain, this reflector passing through the image point x.
The dip vector p is expressed according to:
la=(ax3/axl, a~3/ax2, 1) Equation (7a) where ax3/ax~ represents the apparent dip along the x1 axis and ax3/7x2 represents the apparent dip along the x~ axi s .
A and ~r represent the inclination and azimuth respectively of the dipped reflector passing through the image point ~ at which the illumination is evaluated.
Thus, the following expressions are obtained:
cosh = (1+pl+p2)-1~2 Equation (7b) tanyr = p2/pl Equation (7c) As we have seen above, starting from the Taylor series development of the diffraction travel time (knoian from the isochronic migration map), the following equation is applicable:
td(xr) - td(x) + (oxtd)T(xr -x)+2 (xr -x)T~Ox,xtd~ (xr-x ) + ... Equation ( 2 ) The reflection point xr is firstly estimated along the reflection plane for which the diffraction travel time is stationary (condition for specular reflection).
In order to consider positions only belonging to the dipped plane, a new vector xp is introduced such that:
~:r = x + M. xp = x+ 0 1 xPl Equation (7d) xp2 - h~ -PZ
Note that in this case, as required, the two vectors of the matrix M are actually perpendicular to the dip vector p .
Rewriting equation (2) using equation (7d) gives the following expression: .
td(xp) - td(x) + (~xtd)T M. xp + ~ ..xP .MT.~x,xtd. Mxp +... Equation (7e) namely Td ( xp) - Td ( x ) +b T. xp + ~ . xP . Fxp + ...
Equation (7f) In this equation (7f):
b T = ( ~ xtd ) T M, where b - MT . ( ~ xtd ) , and F - MT
~x~xtd~M.
Consequently, the position xrp for which the diffraction travel time td(xp)~is stationary is deduced:
0xptd(xp) ~ x,p=oT _~ xrp=-F 1~ b =~ xr=x+M. xrp Equation (7g) Equations (7g) and (2) are combined to obtain an expression for the reflection travel time in equation (4) : tr = td(.xr) ,., td(x) - 1 bT.F 1. b .
As already mentioned, the illumination fold is incremented during operation 3 as a function of the difference between the diffraction travel time td and the reflection travel time tr.
Even if the reflection travel time tr(xr(p); s, r) is different from the diffraction travel time td(x;s,r), due to the limited band nature of seismic data, the migration of the reflection that took place at x can constructively contribute to the image of the reflector with a dip along x.
The selection between illuminating migration contributions (constructive interference) and non illuminating migration contributions (destructive interference) is made using the increment function 1 (td, tr; S , ~ ) ~
According to one embodiment of the invention, the increment function is expressed according to i(td,tr;
s,r) - s(td(x;s,r) - tr(xr(p); s, r)), where:
- s(t) is a seismic wavelet which may be monochromatic, with limited band or infinite band and it can be chosen to be dependent on time and space, - s(t) is the derivative with respect to time of the seismic wavelet s (t) - d[s (t) ]/dt.
A derivative with respect to time is usually applied to each seismic trace before the Kirchhoff migration (in time or in depth). Thus, particularly in order. to make the studied quantities uniform, the derivative s(t) of the seismic wavelet with respect to time is chosen instead of using the seismic wavelets s(t) itself.
According to another embodiment of the invention, particularly when the derivatives of seismic traces with respect to time are not determined, the increment function takes account of the seismic wavelet s(t) instead of the derivative s(t) of the wavelet with respect to time; the result is then i (td, tr; s , r ) -S(td(xi~r-~') - tr(xr(p)iSrz'))~
A migration amplitude correction term w(x;s,r) is conventionally calculated during the migration, and in particular takes account of a regular acquisition geometry.
Therefore, this term w(x;s,r) may advantageously be used in the method according to the invention, and act as "preconditioning" to the correction made by determination of the previously described illumination fold.
The increment function may then be weighted by this preconditioning term, such that the illumination fold I ( x, p ; s , ~r ) ) related to a (source s , receiver r ) pair is incremented by:
I(tdrtri ~'rr)~w(~aSr-r').
The evaluation of the illumination fold in the migrated domain according to the invention is optimum since it is applicable to any acquisition geometry.
Based on an inverse approach, the illumination fold I(x,p) estimated by the invention agrees with the theoretical results established by Bleistein. This is not the case for other methods proposed to solve a similar problem.
As already mentioned earlier, the evaluation of the illumination fold proposed by the invention uses only basic quantities available for migration conventionally done by a Kirchhoff summation.
5 Advantageously, the illumination fold according to the invention may be estimated while seismic traces are being migrated.
This is not the case for some other proposed methods (Schneider et al., Laurain and Vinje) for which the said 10 evaluation is done post-migration.
According to another aspect, the invention also proposes a method for correction of the amplitudes of seismic data recorded during 3D seismic prospecting in order to compensate for the effect of non-uniform 15 illumination of subsoil reflectors, characterized in that it comprises steps consisting of:
- estimating the illumination fold I(a:,la) by the method according to the first aspect of the invention described in detail above, 20 - using the inverse I-'~ ( x, p ) of the said ratio as the weighing factor to be applied to each of the said seismic data amplitudes.
According to yet another aspect, the invention also proposes a method for selection of one of a plurality of acquisition geometries as a function of the target of 3D
seismic prospecting, characterized in that it comprises steps consisting of:
- for each acquisition geometry considered, determining the illumination fold I(x,p) estimated according to the first aspect of the invention given in detail above, - selecting the acquisition geometry providing the optimum illumination fold as a function of the target.
Finally, it will be noted that the illumination fold I(x,p) estimated according to the first aspect of the invention described in detail above, may also be used after the migration to eliminate acquisition patterns, and to deconvolute the migrated image.
- Figure 4 is a flowchart diagrammatically showing the main steps in the method according to the invention.
The invention proposes a method of estimating the illumination fold in the migrated domain I ( x , p ; s , F ) for any migrated position x, regardless of its dip component p and regardless of the distribution of seismic traces, using only the most basic quantities involved in migration by a Kirchhoff summation, namely:
- the coordinates (s,:r) of source-receiver pairs, where s represents the position of a source S and z' represents the position of a receiver R;*
- migration isochronic maps td(x;s,r) for each of the (source s, receiver r) pairs involved in the migration.
The invention is equally applicable to imagery in time and imagery in depth.
The invention divulges a method of the type including a modelling step followed by a migration step.
Concerning the modelling step of the reflection point ~ and the reflection travel time tr(.~r), the method according to the invention applies the stationary phase condition to a Taylor series development of diffraction travel times td(x).
The said development is a second order development.
Concerning the migration step, the method according to the invention uses isochronic migration maps available for the migration conventionally done by a Kirchhoff summation .
For each of the (source s, receiver r) pairs, the illumination fold I ( x, p ; s , r ) is updated with reference to Figure 3 using the following equation:
I(x,pis,r) - I(x,pis,r) + I(td,tr; s,r) Equation (1) where is the image point in the migrated domain x -(xlr x2r x3)Ti - p is the dip vector representing the selected dip component, p - (pi r pa r p3 ) T
- s defines the coordinates of the source, s - (s1, s2r S3)Ti - r defines the coordinates of the receiver, r -T, 1 5 (rlr r2r r3) r - td(x; s, r) is the diffraction travel time from the source s, as far as the migrated position x, and return as far as the receiver.r (path indicated by td on Figure 3 and shown in solid lines);
- tr(xr(p)i s, r) is the reflection travel time from the source s as far as point x of specular reflection on the plane reflector passing through image point ~: and perpendicular to the dip vector Via, and return as far as the receiver ~' (path indicated by tr on Figure 3 and shown in dashed lines);
- I(td,tri s, ~') is an increment function for which the expression will be described in detail below.
As already mentioned above, the method according to the invention advantageously uses isochronic migration maps td(x; s, r) involved in the migration by Kirchhoff summation, for each of the ( source s , receiver ~ ) pairs involved in the migration. These maps are "input data" to the method according to the invention.
It also uses first and second derivatives O Xtd ( x ; s , r ) , ~X,Xta ( x ; S . r ) of migration isochrones when they are available. If not, the said derivatives are not input data, but may be calculated as part of the method according to the invention.
Figure 4 is a flow chart showing the different steps in the method according to a preferred embodiment of the invention.
The said method consists of carrying out a step 1 used for each image point x in the migrated 3D domain, for each of the (source s, receiver r) pairs involved in the migration.
The said step 1 includes two operations 2, 3, for each dip component p selected for the illumination calculation I ( x ; p ; s , r ) .
The first operation 2 consists of determining the reflection travel time tr ( xr ( p ) ; s , r ) to the specular reflection point xr(p) on -the plane reflector passing through the image point x and perpendicular to the dip vector p .
~peration 2 is determined from the diffraction travel time td(x; s, r) known from the isochronic migration maps.
The second operation 3 consists of incrementing the illumination fold in the migrated domain I(x,p;s,r) of an increment function I(td, tr; s, r), the said increment function taking account of:
- the difference {td ( x ; s , r ) - tr ( xr ( p ) ; s , r ) ~
between the travel times mentioned above-;
- the seismic wavelet s(t);
- and possibly the migration amplitude term w(x;s.r).
At the end of operations 2 and 3, the illumination 5 amplitude ratio I ( x, p ; s , r ) is known for a given image point x and a dip component p. .
By repeating the operations mentioned above for each of the (source s, receiver r) pairs, the final step is to determine the illumination ratio in the migrated 10 domain I ( x , p ) .
This is done by adding therelative illumination to each pair ( source s , receiver r ) according to I ( x , p ) -~I (x.h%s.-r').
s,;.
The reflection point xr and the reflection travel l5 time tr are determined (during operation 2 mentioned above) starting from a stationary phase analysis, by a second order Taylor series development of the diffraction travel time td around the image point x .
The following expressions for xr and tr (equations (3) and (4) ) are obtained from this analysis (see equation (2) below):
td(~;~,.~'>=tr((~;~s.r>+(~xtr>r(.~-.~r>+2 (.~-%r)T(~~..utr)(x-.~r)+...
Equation (2) .~r( p ) - ~-M.F-1. b Equation (3) tr(xr(~);s.r) - td(x;S.r) y.j7r. F1. j7 Equation (4) where, in equations (3) and (4).
- M is a (3 x 2) matrix described by two vectors extending along the length of the reflection plane.
Consequently, these two vectors are perpendicular to the dip vector p.
- b is a (2 x 1) vector of the second order derivatives of the diffraction travel time along the reflection plane, b - MT . ( 0 Xtd ) Equation ( 5 ) - F is a (2 x 3) matrix of first order derivatives of the diffraction travel time along the reflection plane, F =_ MT- (0 X~xtd) ~M Equation (6) The following describes how equations (3) and (4) are obtained.
Note that the dip vector p represents the apparent dip of a reflector in the migrated domain, this reflector passing through the image point x.
The dip vector p is expressed according to:
la=(ax3/axl, a~3/ax2, 1) Equation (7a) where ax3/ax~ represents the apparent dip along the x1 axis and ax3/7x2 represents the apparent dip along the x~ axi s .
A and ~r represent the inclination and azimuth respectively of the dipped reflector passing through the image point ~ at which the illumination is evaluated.
Thus, the following expressions are obtained:
cosh = (1+pl+p2)-1~2 Equation (7b) tanyr = p2/pl Equation (7c) As we have seen above, starting from the Taylor series development of the diffraction travel time (knoian from the isochronic migration map), the following equation is applicable:
td(xr) - td(x) + (oxtd)T(xr -x)+2 (xr -x)T~Ox,xtd~ (xr-x ) + ... Equation ( 2 ) The reflection point xr is firstly estimated along the reflection plane for which the diffraction travel time is stationary (condition for specular reflection).
In order to consider positions only belonging to the dipped plane, a new vector xp is introduced such that:
~:r = x + M. xp = x+ 0 1 xPl Equation (7d) xp2 - h~ -PZ
Note that in this case, as required, the two vectors of the matrix M are actually perpendicular to the dip vector p .
Rewriting equation (2) using equation (7d) gives the following expression: .
td(xp) - td(x) + (~xtd)T M. xp + ~ ..xP .MT.~x,xtd. Mxp +... Equation (7e) namely Td ( xp) - Td ( x ) +b T. xp + ~ . xP . Fxp + ...
Equation (7f) In this equation (7f):
b T = ( ~ xtd ) T M, where b - MT . ( ~ xtd ) , and F - MT
~x~xtd~M.
Consequently, the position xrp for which the diffraction travel time td(xp)~is stationary is deduced:
0xptd(xp) ~ x,p=oT _~ xrp=-F 1~ b =~ xr=x+M. xrp Equation (7g) Equations (7g) and (2) are combined to obtain an expression for the reflection travel time in equation (4) : tr = td(.xr) ,., td(x) - 1 bT.F 1. b .
As already mentioned, the illumination fold is incremented during operation 3 as a function of the difference between the diffraction travel time td and the reflection travel time tr.
Even if the reflection travel time tr(xr(p); s, r) is different from the diffraction travel time td(x;s,r), due to the limited band nature of seismic data, the migration of the reflection that took place at x can constructively contribute to the image of the reflector with a dip along x.
The selection between illuminating migration contributions (constructive interference) and non illuminating migration contributions (destructive interference) is made using the increment function 1 (td, tr; S , ~ ) ~
According to one embodiment of the invention, the increment function is expressed according to i(td,tr;
s,r) - s(td(x;s,r) - tr(xr(p); s, r)), where:
- s(t) is a seismic wavelet which may be monochromatic, with limited band or infinite band and it can be chosen to be dependent on time and space, - s(t) is the derivative with respect to time of the seismic wavelet s (t) - d[s (t) ]/dt.
A derivative with respect to time is usually applied to each seismic trace before the Kirchhoff migration (in time or in depth). Thus, particularly in order. to make the studied quantities uniform, the derivative s(t) of the seismic wavelet with respect to time is chosen instead of using the seismic wavelets s(t) itself.
According to another embodiment of the invention, particularly when the derivatives of seismic traces with respect to time are not determined, the increment function takes account of the seismic wavelet s(t) instead of the derivative s(t) of the wavelet with respect to time; the result is then i (td, tr; s , r ) -S(td(xi~r-~') - tr(xr(p)iSrz'))~
A migration amplitude correction term w(x;s,r) is conventionally calculated during the migration, and in particular takes account of a regular acquisition geometry.
Therefore, this term w(x;s,r) may advantageously be used in the method according to the invention, and act as "preconditioning" to the correction made by determination of the previously described illumination fold.
The increment function may then be weighted by this preconditioning term, such that the illumination fold I ( x, p ; s , ~r ) ) related to a (source s , receiver r ) pair is incremented by:
I(tdrtri ~'rr)~w(~aSr-r').
The evaluation of the illumination fold in the migrated domain according to the invention is optimum since it is applicable to any acquisition geometry.
Based on an inverse approach, the illumination fold I(x,p) estimated by the invention agrees with the theoretical results established by Bleistein. This is not the case for other methods proposed to solve a similar problem.
As already mentioned earlier, the evaluation of the illumination fold proposed by the invention uses only basic quantities available for migration conventionally done by a Kirchhoff summation.
5 Advantageously, the illumination fold according to the invention may be estimated while seismic traces are being migrated.
This is not the case for some other proposed methods (Schneider et al., Laurain and Vinje) for which the said 10 evaluation is done post-migration.
According to another aspect, the invention also proposes a method for correction of the amplitudes of seismic data recorded during 3D seismic prospecting in order to compensate for the effect of non-uniform 15 illumination of subsoil reflectors, characterized in that it comprises steps consisting of:
- estimating the illumination fold I(a:,la) by the method according to the first aspect of the invention described in detail above, 20 - using the inverse I-'~ ( x, p ) of the said ratio as the weighing factor to be applied to each of the said seismic data amplitudes.
According to yet another aspect, the invention also proposes a method for selection of one of a plurality of acquisition geometries as a function of the target of 3D
seismic prospecting, characterized in that it comprises steps consisting of:
- for each acquisition geometry considered, determining the illumination fold I(x,p) estimated according to the first aspect of the invention given in detail above, - selecting the acquisition geometry providing the optimum illumination fold as a function of the target.
Finally, it will be noted that the illumination fold I(x,p) estimated according to the first aspect of the invention described in detail above, may also be used after the migration to eliminate acquisition patterns, and to deconvolute the migrated image.
Claims (14)
1. Method for estimating the seismic illumination fold in the migrated 3D domain at at least one image point ~, for at least one dip of vector ~, wherein the illumination fold for each (source ~, receiver ~) pair in the seismic survey is estimated, by applying the following steps:
- determination of the reflection travel time from the source ~ to the specular reflection point ~ on the plane reflector passing through the image point ~ and perpendicular to the dip vector ~ and then returning to the reflector ~;
starting from the diffraction travel time t d(~;~,~) from the source ~ to the said image point ~ and then returning to the reflector ~;
- incrementing the said illumination fold related to the (source ~ , receiver ~) pair as a function of the difference between the diffraction travel time and the reflection travel time
- determination of the reflection travel time from the source ~ to the specular reflection point ~ on the plane reflector passing through the image point ~ and perpendicular to the dip vector ~ and then returning to the reflector ~;
starting from the diffraction travel time t d(~;~,~) from the source ~ to the said image point ~ and then returning to the reflector ~;
- incrementing the said illumination fold related to the (source ~ , receiver ~) pair as a function of the difference between the diffraction travel time and the reflection travel time
2. Method according to claim 1, comprising the step of summating each of the said illumination folds related to a (source ~, receiver ~) pair so as to determine the total illumination fold
3. Method according to claim 1 or 2, wherein, during the incrementing step, the illumination fold is incremented using an increment function according to the increment function taking account of the difference between the diffraction travel time and the reflection travel time .
4 . Method according to claim 3, wherein the increment function i is a function of the seismic wavelet s (t) .
5. Method according to claim 4, wherein the increment function i is expressed as a function of the derivative of the seismic wavelet s (t) according to:
6. Method according to claim 4, wherein the increment function i is expressed as a function of the derivative ~ (t) of the seismic wavelet s (t) with respect to time according to:
7 . Method according to any one of claims 3 to 6, in which an a priori correction of the illumination fold is taken into account by migration, comprising the step of incrementing the illumination fold related to a (source ~ , receiver ~ ) pair by .
8 . Method according to any one of claims 1 to 7, wherein the determination step includes the second order Taylor series development of the diffraction travel time around the image point ~ :
9. Method according to claim 8, wherein the specular reflection point is determined along the length of the reflector such that the diffraction travel time at the said specular reflection point is stationary, according to the equation:
10. Method according to claim 8 or 9, wherein the specular reflection point ~, and the reflection travel time are determined according to the following expressions:
where:
- M is a (3 x 2) matrix described by two vectors extending along the length of the reflector, and therefore perpendicular to the dip vector ~;
- ~ is a (2 x 1) vector of first order derivatives of the diffraction travel time along the reflection plane:
- F is a (2 x 2) matrix of second order derivatives of the diffraction travel time along the reflection plane:
where:
- M is a (3 x 2) matrix described by two vectors extending along the length of the reflector, and therefore perpendicular to the dip vector ~;
- ~ is a (2 x 1) vector of first order derivatives of the diffraction travel time along the reflection plane:
- F is a (2 x 2) matrix of second order derivatives of the diffraction travel time along the reflection plane:
11. Method according to any one of claims 1 to 10, wherein the determination step uses isochronic migration maps specified for each (source ~, receiver ~) pair involved in the migration at each image point ~ in the migrated 3D domain.
12. Method according to any one of claims 1 to 11, wherein the seismic illumination fold in the migrated 3D domain is estimated during the Kirchhoff summation migration of seismic data recorded during the 3D seismic prospecting.
13. Method for correction of seismic data amplitudes recorded during 3D seismic prospecting in order to compensate for the effect of non-uniform illumination of sub-soil reflectors, comprising the steps of:
- estimating the illumination fold using the method according to any one of claims 1 to 12, - using the inverse of the ratio as a weighting factor to be applied to each of the said seismic data amplitudes.
- estimating the illumination fold using the method according to any one of claims 1 to 12, - using the inverse of the ratio as a weighting factor to be applied to each of the said seismic data amplitudes.
14. Method for selection of an acquisition geometry among a plurality of acquisition geometries as a function of the target of 3D seismic prospecting, comprising the steps of:
- determining the illumination fold by the method according to any one of claims 1 to 12, for each of the acquisition geometries considered, - selecting the acquisition geometry providing the optimum illumination fold as a function of the target.
- determining the illumination fold by the method according to any one of claims 1 to 12, for each of the acquisition geometries considered, - selecting the acquisition geometry providing the optimum illumination fold as a function of the target.
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| FR0308862 | 2003-07-21 | ||
| FR0308862A FR2858063B1 (en) | 2003-07-21 | 2003-07-21 | METHOD FOR ESTIMATING THE ILLUMINATION COVERAGE RATE IN THE MIGRE DOMAIN |
| PCT/IB2004/002618 WO2005008292A1 (en) | 2003-07-21 | 2004-07-20 | Method of estimating the illumination fold in the migrated domain |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CA2533012A1 CA2533012A1 (en) | 2005-01-27 |
| CA2533012C true CA2533012C (en) | 2014-08-26 |
Family
ID=33560960
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CA2533012A Expired - Fee Related CA2533012C (en) | 2003-07-21 | 2004-07-20 | Method of estimating the illumination fold in the migrated domain |
Country Status (6)
| Country | Link |
|---|---|
| US (1) | US7406383B2 (en) |
| CA (1) | CA2533012C (en) |
| FR (1) | FR2858063B1 (en) |
| GB (1) | GB2419412B (en) |
| NO (1) | NO337503B1 (en) |
| WO (1) | WO2005008292A1 (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104407383A (en) * | 2014-12-09 | 2015-03-11 | 中国石油天然气集团公司 | Parallel three-dimensional seismic illumination and analysis implementation method |
Families Citing this family (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2007143355A2 (en) | 2006-05-17 | 2007-12-13 | Nexus Geosciences Inc. | Diplet-based seismic processing |
| US8194498B2 (en) * | 2008-01-18 | 2012-06-05 | Westerngeco L.L.C. | Using a wave propagator for transversely isotropic media |
| US8699298B1 (en) | 2008-06-26 | 2014-04-15 | Westerngeco L.L.C. | 3D multiple prediction and removal using diplets |
| WO2010082938A1 (en) * | 2009-01-19 | 2010-07-22 | Landmark Graphics Corporation | Data acquisition and prestack migration based on seismic visibility analysis |
| MX2012003959A (en) | 2009-10-02 | 2012-05-08 | Bp Corp North America Inc | Migration-based illumination determination for ava risk assessment. |
| CA2797434C (en) * | 2010-05-12 | 2017-09-19 | Shell Internationale Research Maatschappij B.V. | Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis |
| IT1404170B1 (en) * | 2011-02-10 | 2013-11-15 | Eni Spa | METHOD OF SEISMIC SURVEY OF THE SUBSULE |
| US9116255B2 (en) | 2011-05-27 | 2015-08-25 | Conocophillips Company | Two-way wave equation targeted data selection for improved imaging of prospects among complex geologic structures |
| US9279896B2 (en) | 2011-05-27 | 2016-03-08 | Conocophillips Company | Reciprocal method two-way wave equation targeted data selection for improved imaging of complex geologic structures |
| CA2837649A1 (en) | 2011-09-28 | 2013-04-04 | Conocophillips Company | Reciprocal method two-way wave equation targeted data selection for seismic acquisition of complex geologic structures |
| US9733374B2 (en) | 2012-12-17 | 2017-08-15 | Cgg Services Sas | Repeatability indicator based on shot illumination for seismic acquisition |
| EP2743737A2 (en) | 2012-12-17 | 2014-06-18 | CGG Services SA | Methods and systems for quality control of seismic illumination maps |
| CN104237939B (en) * | 2014-09-12 | 2016-10-26 | 中国科学院地质与地球物理研究所 | A kind of 3 D seismic observation system multi-frequency migration noise acquisition methods and device |
| US10353092B2 (en) * | 2015-12-10 | 2019-07-16 | Pgs Geophysical As | Velocity model update with an inversion gradient |
| CN105467445B (en) * | 2015-12-15 | 2018-02-02 | 中国石油天然气集团公司 | The method for building up and device of a kind of 3 D seismic observation system |
| CN113866821B (en) * | 2021-09-26 | 2022-08-02 | 吉林大学 | Passive source interference offset imaging method and system based on illumination direction constraint |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5537319A (en) * | 1993-12-01 | 1996-07-16 | Schlumberger Technology Corporation | Method for load balancing seismic migration processing on a multiproccessor computer |
| FR2821677B1 (en) * | 2001-03-05 | 2004-04-30 | Geophysique Cie Gle | IMPROVEMENTS TO TOMOGRAPHIC INVERSION PROCESSES OF POINTED EVENTS ON MIGREE SEISMIC DATA |
| US6665618B1 (en) * | 2002-08-14 | 2003-12-16 | Conocophillips Company | Seismic survey design technique |
-
2003
- 2003-07-21 FR FR0308862A patent/FR2858063B1/en not_active Expired - Fee Related
-
2004
- 2004-07-20 US US10/565,403 patent/US7406383B2/en active Active
- 2004-07-20 WO PCT/IB2004/002618 patent/WO2005008292A1/en active Application Filing
- 2004-07-20 CA CA2533012A patent/CA2533012C/en not_active Expired - Fee Related
- 2004-07-20 GB GB0602786A patent/GB2419412B/en not_active Expired - Fee Related
-
2006
- 2006-02-20 NO NO20060818A patent/NO337503B1/en not_active IP Right Cessation
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104407383A (en) * | 2014-12-09 | 2015-03-11 | 中国石油天然气集团公司 | Parallel three-dimensional seismic illumination and analysis implementation method |
| CN104407383B (en) * | 2014-12-09 | 2017-01-18 | 中国石油天然气集团公司 | Parallel three-dimensional seismic illumination and analysis implementation method |
Also Published As
| Publication number | Publication date |
|---|---|
| US7406383B2 (en) | 2008-07-29 |
| WO2005008292A1 (en) | 2005-01-27 |
| FR2858063A1 (en) | 2005-01-28 |
| FR2858063B1 (en) | 2005-09-16 |
| GB0602786D0 (en) | 2006-03-22 |
| CA2533012A1 (en) | 2005-01-27 |
| GB2419412B (en) | 2006-11-01 |
| US20060190179A1 (en) | 2006-08-24 |
| NO337503B1 (en) | 2016-04-25 |
| GB2419412A (en) | 2006-04-26 |
| NO20060818L (en) | 2006-04-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CA2533012C (en) | Method of estimating the illumination fold in the migrated domain | |
| EP2920615B1 (en) | Process for creating image gathers | |
| US6574564B2 (en) | 3D prestack seismic data migration method | |
| CA2521919C (en) | Method for simulating local prestack depth migrated seismic images | |
| US8457899B2 (en) | Method of processing data obtained from seismic prospecting | |
| AU710968B2 (en) | Method and apparatus for seismic signal processing and exploration | |
| CA2775561C (en) | Migration-based illumination determination for ava risk assessment | |
| AU604630B2 (en) | Method of processing seismic data | |
| GB2342165A (en) | 3D prestack seismic data migration method | |
| WO2001096905A1 (en) | Method of frequency domain seismic attribute generation | |
| Thiel et al. | 2D reservoir imaging using deep directional resistivity measurements | |
| Bleistein et al. | Migration/inversion: think image point coordinates, process in acquisition surface coordinates | |
| US6954403B2 (en) | Concurrent phase angle graphic analysis | |
| WO2018087647A1 (en) | Simultaneous source acquisition and separation method | |
| Sun et al. | Wavepath migration versus Kirchhoff migration | |
| US7069200B2 (en) | Method intended to obtain reflection travel times from an interpretation of migrated cylindrical wave seismic data | |
| Brandsberg-Dahl et al. | Beam-wave imaging | |
| Woelz et al. | Shear waves in near surface 3D media–SH–wavefield separation, refraction time migration and tomography | |
| WO2000054074A1 (en) | A method for modeling seismic acquisition footprints | |
| Rickett | Model-space vs data-space normalization for finite-frequency depth migration | |
| Koyan et al. | 3D ground-penetrating radar data analysis and interpretation using attributes based on the gradient structure tensor | |
| Calderón et al. | Efficient computation of 3D DMO responses for 3D acquisition footprint analysis | |
| EA044564B1 (en) | BUILDING A SPEED MODEL | |
| Sigernes et al. | Seismic fault imaging | |
| Rickett | Model-space vs. data-space normalization for recursive depth migration |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| EEER | Examination request | ||
| MKLA | Lapsed |
Effective date: 20180720 |