EP1169655A1

EP1169655A1 - A method for modeling seismic acquisition footprints

Info

Publication number: EP1169655A1
Application number: EP00917892A
Authority: EP
Inventors: Graham A. Winbow; William A. Schneider, Jr.
Original assignee: ExxonMobil Upstream Research Co
Current assignee: ExxonMobil Upstream Research Co
Priority date: 1999-03-12
Filing date: 2000-03-10
Publication date: 2002-01-09
Also published as: NO20014416L; NO20014416D0; WO2000054074A1

Abstract

A model seismic image of a subsurface seismic reflector (216) is constructed, wherein a set of source and receiver pairs is located (202), and a subsurface velocity function (204) is determined. Specular reflection points (206) are determined on the subsurface seismic reflector for each of the source and receiver pairs. A Fresnel zone (208) is determined on the subsurface seismic reflector for each of the specular reflection points using the subsurface velocity function. One or more seismic wavelets are selected (210) and a set of image points is defined containing the subsurface seismic reflector (212). A synthetic seismic amplitude is determined for each of the image points by summing (214) the Fresnel zone synthetic seismic amplitude for all of the Fresnel zones that contain the image point using the seismic wavelets. The model seismic image of the subsurface seismic reflector is constructed (216) using the synthetic seismic amplitudes at the image points.

Description

A METHOD FOR MODELING SEISMIC ACQUISITION FOOTPRINTS

FIELD OF THE INVENTION

This invention relates generally to the field of seismic prospecting. More particularly, the invention relates to a method for constructing a model seismic image of a subsurface seismic reflector.

BACKGROUND OF THE INVENTION

In the oil and gas industry, seismic prospecting techniques are commonly used to aid in the search for and evaluation of hydrocarbon deposits located in subterranean formations. In seismic prospecting, seismic energy sources are used to generate a seismic signal which propagates into the earth and is at least partially reflected by subsurface seismic reflectors. Such seismic reflectors typically are interfaces between subterranean formations having different elastic properties. The reflections are caused by differences in elastic properties, specifically wave velocity and rock density, which lead to differences in impedance at the interfaces. The reflections are recorded by seismic detectors at or near the surface of the earth, in an overlying body of water, or at known depths in boreholes. The resulting seismic data may be processed to yield information relating to the geologic structure and properties of the subterranean formations and potential hydrocarbon content.

The goal of all seismic data processing is to extract from the data as much information as possible regarding the subterranean formations. In order for the processed seismic data to fully represent geologic subsurface properties, the true amplitudes resulting from reflection of the input signal by the geologic target must be accurately represented. This requires that the amplitudes of the seismic data must be processed free from non-geologic seismic effects. Non-geologic amplitude effects include mechanisms that cause the measured seismic amplitudes to deviate from the amplitude caused by the reflection coefficient of the geologic target. These non-geologic amplitude effects can be related to acquisition of the data or to near surface effects. Examples of non-geologic amplitude effects that can be particularly troublesome are source and receiver variation, coherent noises, electrical noise or spikes, and overburden and transmission effects. If uncorrected, these effects can dominate the seismic image and obscure the true geologic picture.

A seismic wave source generates a wave that reflects from or "illuminates" a portion of a reflector. The collection of sources that comprises an entire 3-D survey generally illuminates a large region of the reflector.

Conventional prestack 3-D migration algorithms can produce precise images of the reflector only if illumination is relatively uniform. Lateral velocity variations within the earth and nonuniformly sampled 3-D prestack seismic data, however, generally cause reflectors to be illuminated nonuniformly. Nonuniform illumination is generally due to varying azimuths and source- receiver midpoint locations. Consequently, prestack 3-D migrated images are often contaminated with non-geologic artifacts called an "acquisition footprint". These artifacts can interfere with the ultimate interpretation of seismic images and attribute maps. Understanding and removing the effects of the acquisition footprint has thus become important for seismic acquisition design, processing, and interpretation.

Let S = (S-i, S₂)^τ and G = (G-i, G₂)^τ denote two-dimensional coordinate vectors of a seismic source, commonly called a shot, and a seismic receiver,

typically a geophone, respectively. These vectors S and G are defined with respect to a global Cartesian coordinate system x = (x, y,z) in the plane given

by z = 0. Schleicher et al., "3-D True-Amplitude Finite-Offset Migration", Geophysics, 58, 1112-1126, (1993) showed that for any specified

measurement configuration of sources and receivers, the vector pair (S ,G)

can be described by a single position vector ξ = (ξ-i, ξ₂)^τ according to the

following relations

s(ξ) = s₀ + r_s(ξ - ξ₀)

and

G(ξ) = G₀ + r_G(ξ - ξ₀) . Here S₀ and G₀ are coordinate vectors describing a fixed source-receiver

pair (S₀ , G₀ ) defined by position vector ξ = ξ₀. Configuration matrices Ts

and TG are 2x2 constant matrices, depending only upon the measurement

configuration. The configuration matrices are determined by

and

r_G =

for i = 1 , 2 and j = 1 , 2.

Examples of sets of configuration matrices for particular measurement configurations are:

r_s = I and r_G = I for common offset configuration,

r_s = 0 and r_G = I for common shot configuration, and

Ts = I and TQ = 0 for common receiver configuration.

Here I is the 2x2 identity matrix and 0 is the 2x2 zero matrix. As discussed by Schleicher et al., (1993), supra, common offset 3-D migration can be formulated as a weighted summation along the diffraction

traveltime surface, τ_D , also called the Huygens surface, through the following

Kirchhoff migration equation

U(x) = ∑ Δξ, w(ξ., x) 0(1,, t + τ_D (ξ_{J f} x) ) . (1 )

Here U is the seismic amplitude sum; x is the image point; %_i is the source-

receiver midpoint for trace j; w is a weight chosen to preserve seismic

amplitudes; and τ_D is the diffraction traveltime connecting source S, image

point Q, and receiver G. The dot above U on the right hand side of equation

(1 ) represents time differentiation. The term Δξ. in migration equation (1 ) is

the two-dimensional element of surface area representing trace j, indicating that this equation is a discretized integral formula. Neglect of variations in

Δξ_j , caused by variations in source - receiver midpoint distribution, is in

practice a significant cause of acquisition footprints. Migration implemented by equation (1 ) is referred to in the literature as "True Amplitude Kirchhoff Migration".

Consider an image volume consisting of N points in each of the x, y, and z directions. Then imaging with migration equation (1 ) is an 0(N⁵)

process, because the two-dimensional summations over ξ, must be

performed at each point x of the three-dimensional image space. 0(N⁵) processes such as migration equation (1 ) must usually be implemented on supercomputers or on parallel networks of fast workstations.

Migration equation (1 ) only produces a good image under the assumptions that the source-receiver offset and azimuth are fixed in

magnitude within a narrow range, and the data coverage is dense in ξ -

space. This last assumption means that there are no large data gaps, which

would give large Δξ, values. Failure of the data or the migration

implementation to meet these assumptions generally results in an image that is contaminated with an acquisition footprint. Schleicher et al., (1993), supra, show that synthetic seismic data may be written for one reflector and a collection of sources and receivers as the following data equation

U(ξ., t) = R(ξ.) g(t - τ_R(ξ.)) / Λ(ξ.) . (2)

Here R is the reflection coefficient for the reflector, g(t) represents the seismic

wavelet, TR is the reflection traveltime, and Λ is the geometrical spreading. If

there were more than one reflector, an additional factor to account for the total loss in amplitude due to transmissions across the interfaces along the reflection raypath would be needed. Migration equation (1 ) becomes, upon insertion of synthetic data equation (2), , ,r_{. y}, _{A g} w(ξ_J, x) R(ξ_J) g(t + τ_D(ξ_l, x) - τ_R(ξ_J)) |_t^

U X = > Δξ. =: ^■ (3) i ^J Λ(ξ ⁾

The target-oriented imaging equation (3) can be implemented as an 0(N⁴) process because the image location x may be restricted to only those points in image space that lie on the reflector surface.

In Kirchhoff migration in general, the set of ξ values summed over is

referred to as the "migration aperture". Consider imaging equation (3) for some image point or reflection point x on the reflector. Schleicher et al., "Minimum Apertures and Fresnel Zones in Migration and Demigration", Geophysics, 62, 183-194, (1997), showed that the optimum migration aperture contains only those traces that are reflected from the Fresnel zone area surrounding image point x . This result only applies for migrating synthetic data that are free of diffracting targets, since otherwise a large aperture is needed to adequately focus diffracted energy. The formulas derived by Schleicher et al. (1993) and (1997) apply only to an idealized situation in which offset and azimuth are strictly constant and source - receiver spacing is uniform and dense. In addition, their discussion of Fresnel zones is in relation to theoretical considerations about migration aperture. They did not consider the idea, introduced here for the present invention, of using weighted sums of Fresnel zone amplitudes to construct model seismic images which respect the exact field acquisition geometry in which source - receiver offset and azimuth vary over the survey. The paper "Quantifying Seismic Amplitude Distortions Below Salt" by D.R. Muerdter, R. O. Lindsey, and D. W. Ratcliff: 1997 Offshore Technology conference paper OTC 8339, discusses the effect of irregular reflector illumination. The authors sort the results of raytrace modeling into reflection- point gathers and then prepare and analyze reflection point bin maps around salt structures. However, this method is only an approximation of seismic amplitudes.

Thus, a need exists for a method of removing non-geologic amplitude effects without adversely affecting the quality of the resulting data. There is a need for a method to reveal the acquisition footprint by modeling the migrated image more accurately, more economically, and faster than can be achieved by direct implementation of migration equation (1 ). There is a need for a method that creates maps of the acquisition footprint that are consistent with known wave theory of how seismic images are formed.

SUMMARY OF THE INVENTION

The present invention is a method for constructing a model seismic image of a subsurface seismic reflector. First, a set of source and receiver pairs is located and a subsurface velocity function is determined. Specular reflection points are determined on the subsurface seismic reflector for each of the source and receiver pairs. A Fresnel zone is determined on the subsurface seismic reflector for each of the specular reflection points, using the subsurface velocity function. One or more seismic wavelets are selected. A set of image points is defined containing the subsurface seismic reflector. A

synthetic seismic amplitude is determined at each of the image points by

summing the Fresnel zone synthetic seismic amplitudes for all of the Fresnel

zones that contain the image point, using the seismic wavelets. Finally, the

model seismic image of the subsurface seismic reflector is constructed, using the synthetic seismic amplitudes at the image points.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention and its advantages may be more easily

understood by reference to the following detailed description and the attached

drawings in which:

FIG. 1 is an illustration of the geometry of a dipping planar reflector

relative to source, receiver, and reflection points;

FIG. 2 is a flow chart illustrating the processing steps for one

embodiment of the method of the present invention;

FIG. 3 is a map view of the source - receiver midpoints for the

example;

FIG. 4 is a map view of a portion of the source and receiver locations

for the example;

FIG. 5 is a an illustration of the geometry of the dipping planar

reflections used in the example; FIG. 6 is a plot of 3-D time migrated synthetic data for the planar reflector dipping in the inline direction of the example;

FIG. 7 is a plot of 3-D time migrated synthetic data for the planar reflector dipping 45 degrees between the inline and crossline directions of the example;

FIG. 8 is a plot of 3-D time migrated synthetic data for the planar reflector dipping in the crossline direction of the example.

While the invention will be described in connection with its preferred embodiments, it will be understood that the invention is not limited thereto. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the spirit and scope of the invention, as defined by the appended claims.

DETAILED DESCRIPTION OF THE INVENTION

Traditional seismic image modeling approaches require separate computer runs of modeling software to create synthetic data and of migration software to image the data. These approaches can be cumbersome and time-consuming, often requiring supercomputers. The present invention is developed for modeling 3-D prestack migrated images of single reflectors for realistic acquisition geometries. The present invention combines prestack modeling and imaging into one step, allowing the creation of black box software that is very efficient and easy to use. The present invention, as implemented in the preferred embodiment, is simple and fast, requiring minutes on a workstation. Persons skilled in the art could easily develop computer software for use in implementing the present invention based on the teachings set forth here. In the following, the general theory is described and then acquisition footprint prediction with the present invention for prestack 3-D time migration is demonstrated. The present invention shows how to reveal the acquisition footprint by modeling the migrated image efficiently and accurately. Note, however, that the present invention also applies to general variable velocity media and prestack 3-D depth migration. Thus, with an adequate ray tracing algorithm, wave equation modeling software could be created to efficiently model depth migrated images in real-world situations.

The present invention develops the Fresnel zone theory for variable

source - receiver offset and azimuth. The restriction of the summation over .

for each image point x allows imaging equation (3) to be implemented as an 0(N³) process. This "target-oriented imaging" may be performed efficiently on a single workstation. A particular embodiment of the present invention as a software tool has been produced that specializes the theory to prestack 3-D time migration, for planar reflectors, where the velocity varies with depth only, that is, a v(z) medium. This embodiment models the subsurface illumination or acquisition footprint caused by specific acquisition geometries.

Here the present invention is described for arbitrary velocity models

and planar reflectors. FIG. 1 shows a dipping planar reflector Σ with normal vector 100 making angles α, β, and φ with the x, y, and z axes, respectively.

A reflection path connecting source S and receiver G to specular reflection

point Q on the reflector Σ is illustrated with solid raypath 102. A diffraction

path connecting S and G to nonspecular point Q' on the reflector Σ is

illustrated with dashed raypath 104.

For a wave with pulse length T, the Fresnel zone relative to a surface Σ

is defined as the set of all points Q' on surface Σ for which

|τ(S,Q^,) + τ(Q',G) - τ(S,G)| ≤ T/2 (4)

Here τ(S, Q') is the one-way seismic signal traveltime connecting points S and

Q', τ(Q', G) is the one-way seismic signal traveltime connecting points Q' and

G, and τ(S, G) is the two-way reflection traveltime from points S to Q to G for

specular reflection point Q. Thus, all points Q' on the surface Σ define the

Fresnel zone on that surface, if the sum of traveltimes from S to Q' and from

Q' to G does not differ from the traveltime along specular ray SG by more than

half of pulse length T. Surface Σ may be a planar reflector or the tangent

plane at a specular reflection point Q on the surface of a gently curving

reflector.

The Fresnel zone may be further described in terms of its boundary,

the Fresnel zone ellipse, and using the Fresnel zone matrix H. The Fresnel

zone ellipse is derived by selecting Cartesian coordinates r = (r-t, r₂)^T on the

plane tangent to the reflector Σ at specular reflection point Q so that and r are dip and strike direction, respectively, from Q. The nonspecular reflection point Q', relative to Q, is expressed through nonzero values of r . The

diffraction traveltime TD and reflection traveltime τ_R are defined by

τ_D(S, Q', G) = τ(S, Q') + τ(Q', G), (5)

τ_R(S, G) = τ(S, G) ≡ τ(S, Q) + τ(Q, G). (6)

Expanding diffraction traveltime τo from equation (5) in a second-order Taylor

series in coordinates f = (n, r₂)^T about origin Q yields

τ_D(S,Q',G) « τ_R(S,G) + r^τHr (7)

The first term in this expansion vanishes by Fermat's principle since point Q is the specular reflection point. Inserting this result into the Fresnel zone inequality (4) yields

r^TH? = T (8)

after taking the equality. Taking the equality gives the Fresnel ellipse, the boundary of the Fresnel zone. In the above, H is the 2x2 symmetric Fresnel zone matrix defined by d

H, (9) dr.dτ.

for i = 1,2 and j = 1,2.

The components of H can be explicitly computed as:

5τ_n

- [τ_xxcos²αcos²φ + τ_yy cos² βcos² φ + τ_zz sin φ

5r- sin φ

+ 2τ_xy cos α cos β cos² φ - 2τ_xz cos α cos φ sin² φ (10) - 2τ cos β cos φ sin² φ]

δτ_r D _ _x)cosαcosβcosφ ør..δ 2 A [(τ -τ_x r₂ sin φ

+ τ_xy(cos α-cos β)cosφ (11)

+ τ_xz cos β sin² φ - τ cos α sin² φ]

d τ 1 — - = — — [τ_xxcos²β + τ_wcos²α-2τ cosαcosβ] (12) δr² sin²φ ^{vy xy}

In the Fresnel zone matrix equations (10) - (12), τ_xx, τ_yy, τ_zz, τ_xy, τ_xz, and τ_yz

represent second derivatives of the diffraction traveltime τ_D of equation (5)

with respect to global Cartesian coordinates x = (x,y,z) and evaluated at the

specular reflection point Q on the reflector Σ. Equations (10) - (12) are valid

for general velocity models v(x, y, z). These derivatives are computed for fixed source S and receiver G. The derivatives may be computed precisely, with dynamic ray tracing, or approximately, using finite-differences of the traveltime fields. The three independent components of the symmetric Fresnel zone matrix determine the major and minor elliptical axes and orientation of the Fresnel zone ellipse and thus define the size, shape, and orientation of the Fresnel zone.

Transforming synthetic data equation (2) and imaging equation (3) from

surface midpoint coordinates, ξ_j , for trace j, to Cartesian coordinates, ή. , in

the plane of the reflector, results in

R(ή, ) g(t - τ_R(ή. ))

U( „t) = (13) Λ(ή. ) J(ή . )

and

^J Λ( ,) J( ,)

Here the midpoint coordinate ξ_t has been replaced by the reflection point

for each trace j and J represents the Jacobian of the coordinate transformation

Δή. = J Δξ_r (15)

Transformed synthetic data equation (13) and imaging equation (14) are advantageous because it is particularly easy to restrict the summation over reflection points ή. to the appropriate Fresnel zone. Further, the weight

function w simplifies to

w(ή₎,x) = JΛ(η_J)VdetH (16)

where H = Ht ^x) is the Fresnel zone matrix defined by equations (9) - (12).

Inserting weight equation (16) into imaging equation (14) produces

U(x) = χΔ _] /detHR( _J)g(t + τ_D( _J,x)-τ_R( _J))|_t=0 (17)

Inserting traveltime equation (7) into imaging equation (17) yields

U(x) = ∑Δ _JVdetHR( _J)g(t+r;^THϊ) it=0 (18)

where x_t = ϊji ^x) is the distance vector, in the plane of the reflector, that

connects the reflection point ή_; and the image point x. Imaging equation (18)

represents the best formulation of Fresnel zone imaging theory on the surface

of the imaged reflector. The geometric spreading Λ and the Jacobian J are

eliminated, and the Fresnel zone matrix H controls the amplitudes. Imaging equation (18) precisely and efficiently models the migrated image. Further, imaging equations (17) and (18) are valid for general velocity models. Calculating seismic attributes related to the wavelet in the migrated data requires considering points within an elliptic cylinder which is the extension of the Fresnel zone perpendicular to the reflector surface a distance of approximately a wavelength above and below the reflector. To achieve this with maximum computational efficiency, a first order Taylor expansion is made perpendicular to the reflector and the equation

dτ_n 2cosα

(19)

5n

is employed, where α is the ray incidence angle relative to the normal vector

on the reflector, v is the local seismic velocity and n is the perpendicular

distance from the surface of the reflector. For a small distance δn in the n

direction, that is, for a value of δn less than approximately a wavelength,

equation (18) becomes

U(x) = ∑Δή_J deTH R( _J) (20)

Equation (20) can be used to calculate all seismic attributes on the surface of the reflector given only the rays that reflect directly from the surface of the reflector. Since Fresnel zones vary slowly in a direction perpendicular to the reflector, the Fresnel zone on the reflector surface can be used for such points adjacent the reflector. FIG. 2 is a flow chart that illustrates one embodiment of the method of the present invention for constructing a model seismic image of a subsurface seismic reflector. First, at step 200, a subsurface seismic reflector is selected. Typically, the seismic reflector is selected from among the reflectors of interest to the processing of the seismic data set of interest. Preferably, the seismic reflector is given by specifying a subsurface position and orientation. Second, at step 202, source and receiver locations for a set of source and receiver pairs are determined. Preferably, these source and receiver pairs are selected to provide common-offset gathers in the seismic data set of interest, using standard seismic processing techniques well known to those skilled in the art. These source and receiver pairs will be used in determining specular reflection points below. Next, at step 204, a subsurface velocity function is selected. Preferably, this velocity function is that which is used in the seismic processing of the data set of interest. This subsurface velocity function will be used in determining Fresnel zones below.

Next, at step 206, specular reflection points on the subsurface seismic reflector from step 200 are determined for each of the source-receiver pairs from step 202. Preferably, these specular reflection points are determined by ray tracing. Next, at step 208, a Fresnel zone on the subsurface seismic reflector from step 200 is determined for each of the specular reflection points from step 206, using the subsurface velocity function from step 204. Preferably, the Fresnel zone on the subsurface seismic reflector is determined by size, shape, and orientation of the Fresnel zone ellipse on the subsurface seismic reflector. Preferably, the size, shape, and orientation of the Fresnel zone ellipse are determined from the calculated elements of a Fresnel zone matrix. Preferably, the elements of the Fresnel zone matrix are calculated using equations (10) - (12), above.

Next, at step 210, one or more seismic wavelets are selected. Preferably, the seismic wavelets are selected to match the expected seismic wavelets at depth at the reflector of interest from step 200. Typically, this is determined by examination of the seismic data set of interest. Next, at step 212, a set of image points is defined containing the subsurface seismic reflector from step 200. Preferably, the set of image points is defined depending upon the seismic attributes of interest. Typically, attributes such as amplitude are determined by image points on the reflector, while attributes such as pulse length are determined by including a volume of image points within approximately a spatial wavelength of the reflector.

Next, at step 214, a synthetic seismic amplitude is determined for each of the image points from step 212 by summing Fresnel zone synthetic seismic amplitudes for all of the Fresnel zones from step 208 which contain the image point, using the seismic wavelets from step 210. For a set of image points on the seismic reflector, this summing is preferably done using imaging equation (18) above. However, for a set of image points that include points not on but within approximately a spatial wavelength of the seismic reflector, the summing is preferably done using imaging equations (17) or (20), with equation (20) being computationally more efficient. Use of equations (17), (18), or (20) makes it easy to restrict the summation to reflection points on the appropriate Fresnel zone. Since Fresnel zones vary slowly in a direction perpendicular to the reflector, for such points adjacent the reflector, the

Fresnel zone on the reflector surface is preferably used. Finally, at step 216, the model seismic image of the subsurface seismic reflector from step 200 is constructed using the synthetic seismic amplitudes from step 214 at the image points from step 212. Typically, this construction is accomplished using standard imaging techniques well known to those skilled in the art.

EXAMPLE

Imaging equation (17), (18) and (20) have been implemented for a v(z) medium to test the image modeling technique for time migration applications, which form the majority of cases of practical exploration interest. Imaging equations (17), (18) and (20) can be implemented as code that runs on a single workstation. The code may be specialized for calculating acquisition

footprints or illumination by setting R = 1.0 and holding Δ _j constant. Thus,

any deviation of the output amplitudes from a value of 1.0 represents the footprint. Only source and receiver coordinates, a v(z) velocity function, and the parameters defining a reflecting plane of interest are required to be provided. The program outputs an amplitude horizon slice through the reflector image for plotting, as illustrated in FIGS. 6, 7, and 8, discussed below. Output data sets for seismic attribute calculations may also be generated.

FIG. 3 shows almost 400,000 source-receiver midpoints representing a collection of 2000-meter common-offset traces from a marine seismic survey. Acquisition used a 2-source, 3-cable configuration. FIG. 4 shows the source- receiver positions from just 40,000 of the traces represented in FIG. 3. FIG. 4 illustrates that while source positions 400 are relatively uniform, cable feathering caused significant variations in receiver positions 402. However, missing source lines caused large data gaps 404. FIG. 4 also shows a 180- degree change 406 caused by reversal of the shooting direction.

The 400,000 source and receiver positions were fed into the seismic image modeling program, and using v(z) = 1525 + 0.4 z meters/sec for the velocity function, 3-D migrated images for prestack data from three planar reflectors were obtained. FIG. 5 shows the orientation of these planes to the surface X-Y coordinates 500 of FIGS. 3 and 4. All three planes 501 , 502, and 503 dip at 25 degrees, with the updip directions shown by the arrows 511 , 512, and 513, respectively, in FIG. 5.

FIGS. 6, 7, and 8 are amplitude horizon slices through the migrated images of reflector planes 501 , 502, and 503, respectively, from FIG. 5. These plots were constructed by vertically projecting the migrated image of each reflector onto the surface of the X-Y plane. Each image was computed on a Sun Sparc 10 workstation in a few minutes. The images all show a similar acquisition footprint that appears as vertical stripes in the inline direction. Reflector plane 502 was oblique to the inline direction, which caused the skewed footprint image in FIG. 7. The amplitude variations in these images correlate well with the source-receiver midpoint variations of FIG. 3. Regions of high midpoint density caused amplitude highs, and regions of low midpoint density caused amplitude lows. The images in FIGS. 6, 7, and 8 show highest definition or resolution where the dipping reflectors were shallowest, and this illustrates how Fresnel zones control imaging. These regions of shallowest reflector depth are at the bottom edge of FIG. 6, bottom-left corner of FIG. 7, and the left edge of FIG. 8, recalling the geometry of FIG. 5. Thus, the acquisition footprint is most pronounced on shallow reflectors.

It should be understood that the invention is not to be unduly limited to the foregoing which has been set forth for illustrative purposes. Various modifications and alternatives will be apparent to those skilled in the art without departing from the true scope of the invention, as defined in the following claims.

Claims

What is claimed is:

1. A method for constructing a model seismic image of a subsurface seismic reflector, said method comprising the steps of: locating a set of source and receiver pairs to be used in constructing said model seismic image; determining a subsurface velocity function to be used in constructing said model seismic image; determining specular reflection points on said subsurface seismic reflector for each of said source and receiver pairs; determining a Fresnel zone on said subsurface seismic reflector for each of said specular reflection points, using said subsurface velocity function; selecting one or more seismic wavelets to be used in constructing said model seismic image; defining a set of image points containing said subsurface seismic reflector; determining a synthetic seismic amplitude at each of said image points by summing Fresnel zone synthetic seismic amplitudes for all of said Fresnel zones which contain said image point, using said one or more seismic wavelets; and constructing said model seismic image of said subsurface seismic reflector using said synthetic seismic amplitudes at said image points.

2. The method of claiml , wherein said step of summing said Fresnel zone synthetic seismic amplitude at each said image point is accomplished using the following sum

U(x) = ∑ Δή. Vdet H R(η.) g(t + τ_D(η_J, x) - τ_R(ή.)) |_t=0 ,

where U is said synthetic seismic amplitude sum, x is a image point, . is a

reflection point for trace j, Δ _j is a surface area element for trace j in

Cartesian coordinates in the plane of said subsurface seismic reflector, H is a Fresnel zone matrix, R is a reflection coefficient for said subsurface seismic reflector, g is said seismic wavelet, wherein the dot represents a time

derivative, τ_D is diffraction traveltime, and τ_R is reflection traveltime.

3. The method of claiml , wherein said step of summing said Fresnel zone synthetic seismic amplitude at each said image point is accomplished using the following sum

U(x) = ∑ Δή. Vdet H R(ή.) g(t + ^TH ) |_t=0 ,

where U is said seismic amplitude sum, rj is a distance vector, in the plane of

said subsurface seismic reflector, that connects reflection point and image point x for trace j, Δ . is a surface area element for trace j in Cartesian

coordinates in the plane of said subsurface seismic reflector, H is a Fresnel zone matrix, R is a reflection coefficient for said subsurface seismic reflector, and g is said seismic wavelet, wherein the dot represents a time derivative.

4. The method of claiml , wherein said step of summing said Fresnel zone synthetic seismic amplitude at each said image point is accomplished using the following sum

2cosα „

U(x) = ∑Δή. Λ/drtH R(ή.) g δn + ~ r.Tι H_ι r.

V

where U is said seismic amplitude sum, r| is a distance vector, in the plane of

said subsurface seismic reflector, that connects reflection point ή. and image

point x for trace j, Δfj_j is a surface area element for trace j in Cartesian

coordinates in the plane of said subsurface seismic reflector, H is a Fresnel zone matrix, R is a reflection coefficient for said subsurface seismic reflector,

α is a reflection angle, v is a local seismic velocity at said reflection point . ,

and g is said seismic wavelet, wherein the dot represents a time derivative.

5. The method of claim 1 , wherein said step of determining said Fresnel zone comprises the steps of: calculating a Fresnel zone matrix for each of said specular reflection points on said subsurface seismic reflector; and determining size, shape, and orientation of said Fresnel zone from said Fresnel zone matrix.

6. The method of claim 1 , wherein said subsurface velocity function varies with depth only.