U.S. patent application number 14/828030 was filed with the patent office on 2015-12-10 for process for creating image gathers.
The applicant listed for this patent is TOTAL E&P USA, INC.. Invention is credited to Bertrand Duquet, Huimin Guan, Paul Williamson.
Application Number | 20150355355 14/828030 |
Document ID | / |
Family ID | 50681593 |
Filed Date | 2015-12-10 |
United States Patent
Application |
20150355355 |
Kind Code |
A1 |
Guan; Huimin ; et
al. |
December 10, 2015 |
Process for Creating Image Gathers
Abstract
The process of obtaining seismic data includes deploying a
seismic energy source and seismic receivers, actuating the source,
and detecting seismic energy resulting therefrom at the receivers.
The process further includes digitally sampling seismic energy
detected at the receivers indexed with respect to time to form a
plurality of traces and sorting the traces to form a plurality of
shot gathers. In addition, the process includes applying a depth
migration technique to the shot gathers to generate two images
according to the cross-correlation imaging condition (I.sub.0) and
the gradient-based imaging condition (I.sub.1) for each shot
gather. A reflection angle (.theta.) or a general domain parameter
(.alpha.) is computed at each subsurface position and the images
are mapped according to the corresponding reflection angle or the
general domain parameter to form common image gathers. Amplitude
correction is applied to common image gathers and the images and
common image gathers are stored on non-transitory computer-readable
media.
Inventors: |
Guan; Huimin; (Houston,
TX) ; Williamson; Paul; (Pau, FR) ; Duquet;
Bertrand; (Aberdeen, GB) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
TOTAL E&P USA, INC. |
Houston |
TX |
US |
|
|
Family ID: |
50681593 |
Appl. No.: |
14/828030 |
Filed: |
August 17, 2015 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
13675699 |
Nov 13, 2012 |
9128205 |
|
|
14828030 |
|
|
|
|
Current U.S.
Class: |
367/53 ;
367/38 |
Current CPC
Class: |
G01V 2210/512 20130101;
G01V 1/282 20130101; G01V 2210/51 20130101; G01V 1/28 20130101;
G01V 1/301 20130101; G01V 1/325 20130101; G01V 1/30 20130101; G01V
1/303 20130101; G01V 1/362 20130101 |
International
Class: |
G01V 1/30 20060101
G01V001/30; G01V 1/36 20060101 G01V001/36 |
Claims
1. A process comprising: obtaining seismic data by deploying a
seismic energy source and seismic receivers, actuating the source,
and detecting seismic energy reflected from a plurality of
subsurface positions resulting therefrom at the receivers; forming
a plurality of traces by digitally sampling the seismic energy
detected at the receivers and indexing the seismic energy with
respect to time; sorting the traces to form a plurality of shot
gathers; applying a depth migration technique to the shot gathers
to generate two images according to the cross-correlation imaging
condition (I.sub.0) and the gradient-based imaging condition
(I.sub.1) for each shot gather; computing a reflection angle
(.theta.) or a general domain parameter a at each subsurface
position; mapping the images according to the corresponding
reflection angle or general domain parameter to form common image
gathers; applying amplitude correction to common image gathers; and
storing the images and common image gathers on non-transitory
computer-readable media.
2. The process of claim 1, wherein the step of applying a depth
migration technique is performed using reverse-time migration.
3. The process of claim 1, wherein the step of computing a
reflection angle (.theta.) or a general domain parameter (.alpha.)
at each subsurface position comprises computing the reflection
angle at at least one subsurface position using the following
formula: .theta.({right arrow over
(x)})=cos.sup.-1(I.sub.1/I.sub.0) wherein {right arrow over (x)} is
an image point corresponding to the at least one subsurface
position.
4. The process of claim 1, wherein the step of computing the
reflection angle (.theta.) or a general domain parameter (.alpha.)
comprises computing the general domain parameter using the
following formula: .alpha. = c 0 2 c s c r cos ( .theta. ) = I 1 /
I 0 ##EQU00012## wherein c.sub.0 is the phase velocity along a
symmetry axis, and c.sub.s and c.sub.r are the phase velocities of
the source and receiver wave fields respectively.
5. The process of claim 1, further comprising, prior to the step of
computing a reflection angle (.theta.) or a general domain
parameter (.alpha.), the step of: converting the images to
reflector dip angle domain.
6. The process of claim 5, further comprising, after the step of
mapping the images according to the corresponding reflection angle
or the general domain parameter, the step of: transforming the
images to space domain.
7. The process of claim 1, further comprising, after the step of
computing a reflection angle (.theta.) or a general domain
parameter (.alpha.), the steps of: mapping the source illumination
to the reflection angle domain; and applying illumination
compensation to common image gathers.
8. The process of claim 1, further comprising generating a model of
a subsurface formation including the subsurface positions based at
least in part on the images or common image gathers of the images
and common image gathers.
Description
BACKGROUND
[0001] 1. Field of the Disclosure
[0002] The disclosure relates generally to the field of seismic
data processing. In particular, methods of the disclosure relate to
the extraction of common image gathers (CIGs) in the angle
domain.
[0003] 2. Description of the Related Art
[0004] Seismic surveying may be used to determine structures,
compositions and fluid content of subsurface Earth formations. For
instance, seismic surveying may be used to infer the presence of
useful materials, such as oil and gas, in subsurface Earth
formations. In seismic exploration for oil and gas, the Earth's
subsurface may be illuminated by a seismic source at or near the
surface of the Earth. As used herein, the term "illumination" means
at least that seismic energy from the source is incident on a
subsurface point. Scattered or reflected energy from the
illuminated subsurface point may be recorded by one or more sensors
or receivers deployed for detecting seismic energy originating from
the source. Impedance boundaries are frequently located at
boundaries between Earth formations having different composition.
Waves propagate into the Earth and may be reflected from the
impedance boundaries and travel upwardly until being detected by
seismic sensors. Structure and composition of the Earth's
subsurface may be inferred from the travel time of the reflected
seismic energy, from the geographic position of the source and each
of the sensors, and from the amplitude and phase of the various
frequency components of the reflected seismic energy. Reflected
waves are typically recorded by a receiver array. The receivers can
be positioned on the Earth's surface, on the ocean bottom, towed
near a water surface, or in a well and can be arranged in any
geometrical pattern in two or three dimensions. In a seismic
survey, the source and receiver array are often relocated to a
number of overlapping areas in order to uniformly illuminate the
subsurface in a region.
[0005] FIG. 1 illustrates diagrammatically an example of a survey
of seismic data with a source S of seismic waves and an array of
receivers G. It also shows a point C of the subsurface which is
assumed to contribute to the signal sensed by one of the receivers
G. The horizontal coordinates of point C of the subsurface are
denoted by x, y (or only one spatial coordinate if 2D imaging
instead of 3D imaging is considered), while its depth is denoted by
z. FIG. 1 also provides a simplified representation (dashed lines)
of the propagation of seismic waves from the source S to the point
C and from the point C to the receiver G. The waves may be
refracted at discontinuities of the geological layers where the
impedance changes and reflected or diffracted at different
positions including that of point C.
[0006] The data recorded in a seismic survey include, for each shot
from a source S and for each receiver G, a seismic trace which is a
time series of the signal sensed by the receiver G. The traces for
a number of shots must be transformed to provide an image of the
subsurface which will be the result of stacking or integrating a
large amount of information. An important step of the
transformation is the migration which consists in processing the
data with respect to a model such that the data received at the
surface are mapped into subsurface to represent the reflectivity of
the discontinuities in the model. The model is usually a map of the
propagation velocity of the acoustic waves in the subsurface. This
model is not known a priori and it is a main challenge of many
seismic imaging technologies to determine a model that will
properly account for the field data.
[0007] In pre-stack depth migration (PSDM) methods, migrated data
are computed for each shot using the velocity model and arranged in
an output cube containing migrated values associated with positions
in the subsurface. The cubes obtained for different shots may then
be stacked and/or analyzed to form Common Image Gathers (CIGs) and
check consistency of the model. The model may be corrected and the
process iterated until a satisfactory image is obtained.
[0008] CIGs are popular tools for evaluating the migration velocity
field, for subsurface attribute analysis, and for imaging
enhancement. CIGs are created either during the migration process
or from data extracted from the output cubes of migration, sorted
in a convenient way for analysis so as to check the velocity model.
A CIG is a bi-dimensional data structure defined for a given
horizontal position (x, y), with a first axis representing the
depth z and a second axis representing a domain parameter a
referred to for sorting the multiple images of the migration
process. It contains reflectivity values obtained from the output
cubes resulting from the migration, forming an image which can be
analyzed to check and/or correct the velocity model. Examples of
commonly used domain parameters .alpha. include the surface offset,
subsurface offset or the scattering or reflection angle at the
subsurface position (x, y, z), etc. Because CIGs in the scattering
angle domain closely represent the angle-dependent reflectivity of
subsurface reflectors, they are the most useful for velocity model
analysis, subsurface attribute analysis and image improvement in
complex media.
[0009] The computation of CIGs is not simple in all wave-field
extrapolation methods and it can be very expensive for reverse-time
migration (RTM), which by itself is already a process requiring a
large amount of computation time and computer memory. RTM is a
two-way migration solution which can accurately describe wave
propagation in complex media. It is increasingly used in seismic
exploration by virtue of advances in computer power and
programming.
[0010] In "Offset and angle-domain common image-point gathers for
shot-profile migration", Geophysics, Vol. 67, No. 3, 2002, pp.
883-889, J. Rickett and P. Sava established the notion of
subsurface offset CIGs which requires the extension of the imaging
condition through the computation of the correlation function along
the spatial horizontal dimension. This type of gather is the most
common way to output wave-equation-based migration images. It
corresponds to velocity updating techniques based on focusing
analysis, which look for the highest correlation at zero-time lag
and/or zero-offset, and only small values elsewhere. A method has
been proposed to derive scattering angle CIGs from subsurface
offset CIGs by applying local slant stack. A similar method has
also been proposed to extend the imaging condition in the time
domain and to convert the time-shifted image gathers to CIGs into
the angle domain.
[0011] The local offset and time lag gathers are formed during the
migration process, and significant computer memory is required to
store the intermediate CIGs. For 3-D cases, a 5-D array is
necessary where local offsets in both x and y directions are used;
if all space and time lags are included, a 7-D array is necessary.
While applying an imaging condition may normally require only a
small proportion of the memory needed in the whole computation of
RTM, the memory requirements may increase dramatically if a large
number of lags are computed in more than one dimension. Another
major part of the computation cost is the use of slant stacking to
extract angle gathers from local offset (or time lag) gathers.
While local slant stack may be an efficient procedure in 2-D, the
gather conversion is very complex in 3-D. While potentially less
expensive to derive angle gathers from CIGs with only non-zero time
lag the resolution of angle gathers from this approach is not as
good as those provided by local offset gathers. Furthermore, this
methodology does not provide azimuth information for 3-D cases.
[0012] In another example, angle domain CIGs can be formed during
the migration process by wave field decomposition. Source and
receiver wave fields in the frequency domain are transformed into
the wave-number domain, cross-correlation is applied between each
component of the source wave field and that of the receiver wave
field, and then the resulting partial images are mapped and binned
according to the corresponding azimuth and reflection angle to form
angle-domain CIGs. Wave field decomposition can be performed in the
time and space domain by local slant stack as well. A 4-D
spatial/temporal Fourier transform is applied to both the source
and receiver wave fields to convert them into frequency wave-number
domain. The procedure to generate angle dependent partial images
typically involves an expensive multi-dimensional convolution of
seven loops. A 3-D inverse spatial Fourier transform is applied to
the angle dependent image formed initially in the wave-number
domain. Spatial windowing, ALFT (anti-aliasing Fourier transform)
and the facts that the norm of slowness in a small window in a
homogeneous media varies only in a small range and the seismic
events are closer to linear help to reduce the cost to make the
approach feasible in practice. The cost of this migration has been
estimated to be 5-10 times the cost of the RTM itself. Besides
intensive computation, this method requires significant disk space
to store the wave fields at all time steps and significant memory
to store the 5-D angle gathers.
[0013] In still another example, instead of wave-field
decomposition, the direction of wave propagation is computed at
each time step, partial images are computed together with the
corresponding angles, and then mapped accordingly into angle
gathers. This approach determines the dominant direction of wave
propagation during the migration process. The additional
computation cost is caused by computing the Poynting vectors of the
wave fields, which is similar to that of a RTM. However when the
wave field is complex it is difficult to detect the direction of
wave propagation accurately.
SUMMARY
[0014] A process is disclosed. The process includes obtaining
seismic data by deploying a seismic energy source and seismic
receivers, actuating the source, and detecting seismic energy
resulting therefrom at the receivers. The process further includes
digitally sampling seismic energy detected at the receivers indexed
with respect to time to form a plurality of traces and sorting the
traces to form a plurality of shot gathers. In addition, the
process includes applying a depth migration technique to the shot
gathers to generate two images according to the cross-correlation
imaging condition (I.sub.0) and the gradient-based imaging
condition (I.sub.1) for each shot gather. A reflection angle
(.theta.) or a general domain parameter (.alpha.) is computed at
each subsurface position and the images are mapped according to the
corresponding reflection angle or the general domain parameter to
form common image gathers. Amplitude correction is applied to
common image gathers, and the images and common image gathers are
stored on non-transitory computer-readable media.
BRIEF DESCRIPTION OF THE DRAWINGS
[0015] The present disclosure is best understood from the following
detailed description when read with the accompanying figures. It is
emphasized that, in accordance with the standard practice in the
industry, various features are not drawn to scale. In fact, the
dimensions of the various features may be arbitrarily reduced for
clarity of discussion.
[0016] FIG. 1 is a schematic diagram illustrating the acquisition
of seismic data.
[0017] FIG. 2 is a graphical depiction of a layered model with one
reflector in accordance with one embodiment of the present
disclosure.
[0018] FIG. 3A is a common image gather at point A in FIG. 2 at
correct velocity.
[0019] FIG. 3B is a common image gather at point A in FIG. 2 at
velocity 10% smaller than correct velocity.
[0020] FIG. 3C is a common image gather at point A in FIG. 2 at
velocity 10% larger than correct velocity.
DETAILED DESCRIPTION
[0021] It is to be understood that the following disclosure
provides many different embodiments, or examples, for implementing
different features of various embodiments. Specific examples of
components and arrangements are described below to simplify the
present disclosure. These are, of course, merely examples and are
not intended to be limiting. In addition, the present disclosure
may repeat reference numerals and/or letters in the various
examples. This repetition is for the purpose of simplicity and
clarity and does not in itself dictate a relationship between the
various embodiments and/or configurations discussed. Moreover, the
formation of a first feature over or on a second feature in the
description that follows may include embodiments in which the first
and second features are formed in direct contact, and may also
include embodiments in which additional features may be formed
interposing the first and second features, such that the first and
second features may not be in direct contact.
[0022] In certain embodiments of the present disclosure, the method
includes obtaining seismic data by deploying a seismic energy
source and seismic receivers, actuating the source, and detecting
seismic energy resulting therefrom at the receivers. The seismic
energy is detected and digitally sampled at the receivers, and
indexed with respect to time to form traces. The traces may then be
sorted to form a plurality of shot gathers. A depth migration
technique may then be applied to the shot gathers to generate two
images according to the cross-correlation imaging condition
(I.sub.0) and the gradient-based imaging condition (I.sub.1) for
each shot gather, i.e., one with and one without the reflection or
scattering angle information are generated for each shot gather in
the migration process.
[0023] Depth migration may be applied to seismic data in order to
convert it from traces in time, additionally indexed by the
coordinates of source and receiver locations to a subsurface image
in space coordinates, including depth. This method may require a
velocity model. Any suitable depth migration technique may be used
for the present disclosure. Non-limiting depth migration techniques
include reverse time migration (RTM) and wave-equation
migration.
[0024] Reflection angles or general domain parameters, which are
functions of reflection angles, may then be computed via a division
between these two images. Pre-stack images may then be mapped onto
angle domain or the domain of general domain parameters to form
CIGs. Amplitude correction may be applied afterwards. Amplitude
correction is a process applied to migration images or image
gathers to make the amplitude of the images or image gathers
proportional to the reflectivity of the subsurface reflector, where
possible.
[0025] In certain embodiments, the present disclosure may be
applied to isotropic media. In certain other embodiments, the
present disclosure may be applied to anisotropic media. In certain
embodiments, the angle gathers keep the high resolution of the
original images.
[0026] In certain embodiments of the present disclosure, the
following algorithms may be used. As disclosed in Stolk, C. C., M.
V. De Hoop and T. J. P. M. OP'T Root, 2009, Inverse scattering of
seismic data in the reverse time migration (RTM) approach,
Proceedings of the Project Review, Geo-Mathematical Imaging Group,
Purdue University, Vold, 91-108, which is herein incorporated by
reference in its entirety, based on an inverse scattering
formulation, an imaging condition for RTM is
I ( x .fwdarw. ) = .pi. .intg. R 1 .omega. 3 p s ( x .fwdarw. ,
.omega. ) 2 j = 0 n ( k js k jr p _ s ( x .fwdarw. , .omega. ) p r
( x .fwdarw. , .omega. ) ) .omega. ( 1 ) ##EQU00001##
in which p.sub.s({right arrow over (x)}, .omega.) is the complex
conjugate of the source wave field; p.sub.r({right arrow over (x)},
.omega.) is the receiver wave field; k.sub.0s=k.sub.0r=.omega./c,
.omega. is the circular frequency, c is the phase velocity of wave
propagation; and k.sub.js and k.sub.jr (j=1, 2, 3) are the
components of wave number in j direction for the source and
receiver wave fields respectively, n is the number of spatial
dimension, z denotes an image point.
[0027] Equation (1) is equivalent to
I ( x .fwdarw. ) = .pi. .intg. R 1 .omega. 3 p s ( x .fwdarw. ,
.omega. ) 2 ( .omega. 2 c 2 ( x ) p _ s ( x .fwdarw. , .omega. ) p
r ( x .fwdarw. , .omega. ) - .gradient. p _ s ( x .fwdarw. ,
.omega. ) .gradient. p _ r ( x .fwdarw. , .omega. ) ) .omega. ( 2 )
##EQU00002##
[0028] Ignoring source-side illumination compensation, the imaging
condition has a simpler form as follows.
I ( x .fwdarw. ) = .intg. R ( .omega. 2 c 2 ( x .fwdarw. ) p _ s (
x .fwdarw. , .omega. ) p r ( x .fwdarw. , .omega. ) - .gradient. p
_ s ( x .fwdarw. , .omega. ) .gradient. p r ( x .fwdarw. , .omega.
) ) .omega. ( 3 ) ##EQU00003##
[0029] I({right arrow over (x)}) can be separated into two
parts:
I ( x .fwdarw. ) = I 0 ( x .fwdarw. ) + I 1 ( x .fwdarw. ) ( 4 ) I
0 ( x .fwdarw. ) = .intg. R .omega. 2 c 2 ( x .fwdarw. ) p _ s ( x
.fwdarw. , .omega. ) p r ( x .fwdarw. , .omega. ) .omega. ( 5 ) I 1
( x .fwdarw. ) = .intg. R - .gradient. p _ s ( x .fwdarw. , .omega.
) .gradient. p r ( x .fwdarw. , .omega. ) ) .omega. ( 6 )
##EQU00004##
[0030] The imaging condition presented in Equation (5) is referred
to as the cross-correlation imaging condition, while the imaging
condition presented in Equation (6) is referred to as the
gradient-based imaging condition. Defining .theta.({right arrow
over (x)}) as the angle between the source and receiver wave-number
vectors, for isotropic case, equations (5) and (6) may be rewritten
as
I 0 ( x .fwdarw. ) = .intg. R k 2 p _ s ( x .fwdarw. , .omega. ) p
r ( x .fwdarw. , .omega. ) .omega. ( 7 ) I 1 ( x .fwdarw. ) =
.intg. R - .gradient. p _ s ( x .fwdarw. , .omega. ) .gradient. p r
( x .fwdarw. , .omega. ) ) .omega. = .intg. R k 2 cos .theta. p _ s
( x .fwdarw. , .omega. ) p r ( x .fwdarw. , .omega. ) .omega. ( 8 )
##EQU00005##
[0031] In which, k=.omega./c, is the wave number. Equation (7) is
similar to that of the conventional cross-correlation imaging
condition, and it is proportional to the cross-correlation of the
first-order derivatives with respect to time of the source and
receiver wave fields. The local reflection angle can be derived
as
.theta.({right arrow over
(x)})=cos.sup.-1(I/I.sub.0-1)=cos.sup.-1(I.sub.1/I.sub.0) (9)
[0032] Equations (3)-(6) are extended to anisotropic media, such as
the transverse isotropic media with a vertical symmetry axis (VTI),
the transverse isotropic media with a tilted symmetry axis (TTI),
the orthorhombic anisotropic media etc. by replacing c with
c.sub.0, the phase velocity along the symmetry axis. It can be
derived that
I 1 / I 0 = c 0 2 c s c r cos ( .theta. ) ( 10 ) ##EQU00006##
[0033] In which, c.sub.s and c.sub.r are the phase velocities of
the source and receiver wave fields, respectively. The right side
of the equation,
c 0 2 c r c s cos ( .theta. ) , ##EQU00007##
is referred to as the general domain parameter a for common image
gathers. The division process comprises minimizing a cost function
defined by a reflection angle and local values of images generated
by using the imaging conditions I.sub.1 and I.sub.0 in a
neighbourhood of each subsurface position.
[0034] For isotropic media, image I can be derived by applying the
Laplacian operator to the image obtained by cross correlating the
source and receiver wave field:
I.sub.2({right arrow over (x)})=.intg..sub.0.sup.t.sup.maxS({right
arrow over (x)},t)R(x,t.sub.max-t)dt (11)
In which S and R are the source and receiver wave field in the time
and space domain, respectively; t.sub.max is the maximum length of
time of the shot gather being migrated. A Laplacian filter is
applied to the image:
.DELTA.I.sub.2({right arrow over (x)})={right arrow over
(.gradient.)}{right arrow over (.gradient.)}I.sub.2 (12)
[0035] The derivative operators are permuted with the integration
in time and derived under the integral sign:
.DELTA.I.sub.2({right arrow over
(x)})=.intg..sub.0.sup.t.sup.max[.DELTA.S({right arrow over
(x)},t)R({right arrow over (x)},t.sub.max-t)+2{right arrow over
(.gradient.)}S({right arrow over (x)},t){right arrow over
(.gradient.)}R({right arrow over (x)},t.sub.max-t)+S({right arrow
over (x)},t).gradient.R({right arrow over (x)},t.sub.max-t)]dt
(13)
[0036] Both fields obey the wave equation:
{ .DELTA. S ( x .fwdarw. , t ) = 1 c ( x .fwdarw. ) 2
.differential. 2 S ( x .fwdarw. , t ) .differential. t 2 .DELTA. R
( x .fwdarw. , t ma x - t ) = 1 c ( x .fwdarw. ) 2 .differential. 2
R ( x .fwdarw. , t ma x - t ) .differential. t 2 ( 14 )
##EQU00008##
[0037] Considering:
.differential. 2 .differential. t 2 S ( x .fwdarw. , t ) R ( x
.fwdarw. , t ma x - t ) + S ( x .fwdarw. , t ) .differential. 2
.differential. t 2 R ( x .fwdarw. , t ma x - t ) = .differential. 2
.differential. t 2 ( S ( x .fwdarw. , t ) R ( x .fwdarw. , t ma x -
t ) ) - 2 .differential. .differential. t S ( x .fwdarw. , t )
.differential. .differential. t R ( x .fwdarw. , t ma x - t ) ( 15
) ##EQU00009##
an integration by parts is performed of the time integral:
.intg. 0 t ma x 1 c ( x .fwdarw. ) 2 ( .differential. 2
.differential. t 2 S ( x .fwdarw. , t ) R ( x .fwdarw. , t ma x - t
) + S ( x .fwdarw. , t ) .differential. 2 .differential. t 2 R ( x
.fwdarw. , t ma x - t ) ) t = 1 c ( x .fwdarw. ) 2 .differential.
.differential. t ( S ( x .fwdarw. , t ) R ( x .fwdarw. , t ma x - t
) ) | 0 t ma x + .intg. 0 t ma x 1 c ( x .rarw. ) 2 2 S . ( x
.fwdarw. , t ) R . ( x .fwdarw. , t ma x - t ) t ( 16 )
##EQU00010##
[0038] The first term on the right side of Equation 14 is zero in
the embodiments of the present disclosure. Applying Equation (12)
and (14) to Equation (13):
- 1 2 .DELTA. I 2 ( x .fwdarw. ) = - .intg. 0 t ma x [ 1 c ( x
.fwdarw. ) 2 S . ( x .fwdarw. , t ) R . ( x .fwdarw. , t ma x - t )
+ .gradient. .fwdarw. S ( x .fwdarw. , t ) .gradient. .fwdarw. R (
x .fwdarw. , t ma x - t ) ] t ( 17 ) ##EQU00011##
The right-hand side of equation (17) is equivalent to the imaging
condition presented in Equation (3). The equivalence allows the
computation of I by applying the Laplacian operator to I.sub.2.
[0039] In certain embodiments of the present disclosure, the
computation cost added to the migration process for isotropic media
is minimal. For anisotropic media, the computation cost increases
by less than one time of the migration cost due to the computation
of the gradients of both the source and receiver wave fields.
[0040] Equation (9) or (10) may only be accurate when an image
point is illuminated at most once by a given shot. In certain
embodiments, such as for a complicated model, this condition may
not be fulfilled; if acoustic/elastic waves arrive at a point in
the subsurface from a given shot at two or more distinct times
and/or from two or more directions, Equation (9) or (10) may not be
accurate. To alleviate the interference of such multiple arrivals,
every image point may be decomposed into images at different dip
angles and the reflection angle may be computed for each dip
component. Then the individual components are mapped to different
reflection angle bins. Such a process separates multiple arrivals
contributing to an image point with different dip angles. Those
multiple arrivals contributing to an image point with the same dip
angle are not separated, so, in certain embodiments, the computed
reflection angle can be a weighted average, close to that of the
most energetic arrival among them.
[0041] In certain embodiments, Equation (9) or (10) may be
implemented by: [0042] (1) For a shot gather, applying RTM and
generating two images according to the imaging conditions, I.sub.1
and I.sub.0; The images according to I.sub.1 and I.sub.0 are in the
(x, y, z) domain. [0043] (2) At each image point, decomposing both
images from the (x, y, z) domain to the (x, y, z, p) domain, in
which p represents the reflector dip angle. [0044] (3) Computing
the reflection angle with Eq. (9) or the general domain parameter a
with Eq. (10) for each p-component at an image point followed by
mapping the images according to the corresponding reflection angle
or general domain parameter, and stacking every component in each
reflection angle or general domain parameter bin. In certain
embodiments, equation (9) and (10) are implemented by least-squares
division. In other embodiments, simple division may be used. [0045]
(4) Repeating procedure (1)-(3) for every shot, stacking
contributions from all shots. [0046] (5) At each image point,
mapping the source illumination similarly according to the computed
reflection angle to the angle domain or according to the general
domain parameter to the a domain for every shot, and stacking
contributions from all shots to form gathers of source
illumination. [0047] (6) Applying illumination compensation to the
common image gathers by, for each angle or .alpha. bin at each
position, dividing the stacked image amplitude from step (4) by the
corresponding stacked illumination from step (5). In certain other
embodiments, Equation (9) or (10) may be implemented by: [0048] (1)
For a shot gather, applying RTM and generating two images according
to the imaging conditions I.sub.1 and I.sub.0. [0049] (2) At each
image point, decomposing the two images in (x, y, z) domain to (x,
y, z, p) domain. [0050] (3) Applying the Hilbert transform to each
trace from both images in the (x, y, z, p) domain, computing the
reflection angle with Equation (9) or the general domain parameter
a with Equation (10) for each p-component, followed by mapping the
images according to the corresponding reflection angle or domain
parameter a, and stacking all contributions to each reflection
angle bin or .alpha. bin. In certain embodiments, least-squares
division may be used to implement equation (9) and (10). In other
embodiments, simple division may be used. [0051] (4) Repeating
procedure (1)-(3) for every shot, stacking contributions from all
shots. [0052] (5) At each image point, mapping the source
illumination similarly according to the computed reflection angle
to the angle domain or according to the general domain parameter to
the a domain for every shot, and stacking contributions from all
shots to form gathers of source illumination. [0053] (6) Applying
illumination compensation to the common image gathers by dividing
the value of the stacked image amplitude by the corresponding value
of stacked illumination.
[0054] In certain other embodiments, Equation (9) or (10) may be
implemented by computing the reflection angle or the general domain
parameter with a weighted average of the values produced in the
previous two implementations.
[0055] In certain embodiments, I.sub.1 and I.sub.0 are in the (x,
y, z) domain. As one of ordinary skill in the art with the benefit
of this disclosure will realize, I.sub.1 and I.sub.0 could also be
represented in other domains.
[0056] As indicated above, in certain embodiments, illumination
compensation may be applied to the common image gathers. In
illumination compensation, the migration image is normalized by
source-side illumination, which is defined as the auto-correlation
of the source wave field and measures the strength of energy from
the seismic source at subsurface.
[0057] In certain embodiments, before the step of computing the
reflection angle or the general domain parameter, the images may be
decomposed into reflector dip angle domain by local slant stack or
.tau.-p transform, herein .tau. represents depth. In certain
embodiments, after computing the reflection angle or the general
domain parameter, the source illumination is decomposed by mapping
to the reflection angle domain. Illumination compensation may then
be applied to the common image gathers by any suitable methods.
[0058] The images, source illumination, and common image gathers
may then be stored on non-transitory computer-readable media or
printed.
[0059] FIGS. 2 and 3A-C are indicative of the process used in at
least one embodiment of the present disclosure. FIG. 2 is a
preliminary graphical depiction of a layered model with one
reflector. In the model, the reflection point can be point A, which
is located on a portion of the reflector without a dip angle. The
axes for the model are x[.DELTA.x] and z[.DELTA.z]. Density of the
upper media is 2000 kg/m.sup.3 with an assumed velocity of 2000
m/s. Density of the lower media is 2000 kg/m.sup.3 with an assumed
velocity of 2400 m/s.
[0060] FIGS. 3A-C show example CIGs developed in accordance with
one embodiment of the present disclosure. FIG. 3A depicts the CIGs
when the model has the correct velocity. FIG. 3B depicts the CIGs
when the velocity is 10% smaller than the correct velocity with
FIG. 3C depicting the CIGS when the velocity is 10% larger than the
correct velocity. The amount of curvature of the line shown in
FIGS. 3B and 3C allows adjustment of the model velocity to attain
the correct velocity.
[0061] The foregoing outlines features of several embodiments so
that a person of ordinary skill in the art may better understand
the aspects of the present disclosure. Such features may be
replaced by any one of numerous equivalent alternatives, only some
of which are disclosed herein. One of ordinary skill in the art
should appreciate that they may readily use the present disclosure
as a basis for designing or modifying other processes and
structures for carrying out the same purposes and/or achieving the
same advantages of the embodiments introduced herein. One of
ordinary skill in the art should also realize that such equivalent
constructions do not depart from the spirit and scope of the
present disclosure, and that they may make various changes,
substitutions and alterations herein without departing from the
spirit and scope of the present disclosure.
[0062] The Abstract at the end of this disclosure is provided to
comply with 37 C.F.R. .sctn.1.72(b) to allow the reader to quickly
ascertain the nature of the technical disclosure. It is submitted
with the understanding that it will not be used to interpret or
limit the scope or meaning of the claims.
[0063] Moreover, it is the express intention of the applicant not
to invoke 35 U.S.C. .sctn.112, paragraph 6 for any limitations of
any of the claims herein, except for those in which the claim
expressly uses the word "means" together with an associated
function.
* * * * *