GEOPHYSICS, VOL. 74, NO. 5 共SEPTEMBER-OCTOBER 2009兲; P. R85–R97, 16 FIGS.

10.1190/1.3190131

Ray-based stochastic inversion of prestack seismic

data for improved reservoir characterization

Dennis van der Burg1, Arie Verdel2, and Kees Wapenaar3

With inversion, detailed information about a reservoir that is avail-

ABSTRACT able at well locations is extrapolated away from the wells to all loca-
tions in the reservoir on the basis of seismic reflections. Examples of
Trace inversion for reservoir parameters is affected by an- reservoir properties that can be inverted for are thickness, wave
gle averaging of seismic data and wavelet distortion on the propagation velocity, porosity, and pore-fluid type of individual lay-
migration image. In an alternative approach to stochastic
ers. The initial layered reservoir model is derived from well-log data,
trace inversion, the data are inverted prestack before migra-
seismic reflection picks, and geologic information. The wavelet re-
tion using 3D dynamic ray tracing. This choice makes it pos-
sible to interweave trace inversion with Kirchhoff migration. quired for inversion is derived from a seismic-to-well match 共White
The new method, called ray-based stochastic inversion, is a and Simm, 2003兲. For a recent overview of inversion techniques, see
generalization of current amplitude versus offset/amplitude Veeken and Da Silva 共2004兲.
versus angle 共AVO/AVA兲 inversion techniques. The new Stochastic trace inversion is a generic term for a specific subset of
method outperforms standard stochastic inversion tech- trace inversion techniques that quantifies uncertainties in reservoir-
niques in cases of reservoir parameter estimation in a struc- parameter estimates. This is achieved by applying Bayes’ rule
turally complex subsurface with substantial lateral velocity 共Duijndam, 1987; Leguijt, 2001; Tarantola, 2005兲. Consequently,
variations and significant reflector dips. A simplification of the method also requires initial uncertainty estimates of all model
the method inverts the normal-incidence response from res- parameters based on prior knowledge. In Bayesian inversion, for-
ervoirs with approximately planar layering at the subsurface
ward modeling a trace from the migration image and the subsequent
target locations selected for inversion. It operates along ray-
comparison of modeled and actual trace is done many times to evalu-
paths perpendicular to the reflectors, the direction that offers
optimal resolution to discern layering in a reservoir. In a test ate the uncertainties in the reservoir parameters. This procedure is
on field data from the Gulf of Mexico, reservoir parameter es- repeated for every trace in the inversion range. A Markov-chain
timates obtained with the simplified method, the estimates Monte Carlo algorithm 共Sambridge and Mosegaard, 2002兲 can be
found by conventional stochastic inversion, and the actual used to generate parameter updates to efficiently sample the posteri-
values at a well drilled after the inversion are compared. Al- or probability density function 共PDF兲 in a statistically valid way. The
though the new method uses only 2% of the prestack data, the method is a global optimization technique. Local optimization tech-
result indicates it improves accuracy on the dipping part of niques such as conjugate gradient often have difficulties because
the reservoir, where conventional stochastic inversion suffers they cannot escape local extremes that might be present in the corre-
from wavelet stretch caused by migration. sponding posterior PDF.
The forward-modeling stage in conventional trace inversion is de-
picted in greater detail in the process flow at the bottom of Figure 1.
Initial estimates of rock and pore-fluid properties acquired from the
INTRODUCTION
layered reservoir model at the current trace position are inserted into
Techniques that invert seismic data can be used to estimate rock 1D rock/fluid property models appropriate for each identified layer
and pore-fluid properties of oil and gas reservoirs in the subsurface. in the inversion target. Using these rock models, the elastic-layer

Manuscript received by the Editor 24 November 2008; revised manuscript received 27 February 2009; published online 18 September 2009.
1
Formerly Delft University of Technology, Delft, The Netherlands; presently Petroleum Geo-Services B. V., Leiden, The Netherlands. E-mail: [Link]
.berg@[Link].
2
Shell International Exploration and Production B.V., Rijswijk, The Netherlands. E-mail: [Link]@[Link].
3
Delft University of Technology, Department of Geotechnology, Delft, The Netherlands. E-mail: [Link]@[Link].
© 2009 Society of Exploration Geophysicists. All rights reserved.

R85

Downloaded from [Link]

by CNOOC user
 R86 van der Burg et al.

properties — P-wave velocity VP, S-wave velocity VS, and density ␳ move some of the limitations of 1D convolution. For example, 1D
are calculated. Subsequently, using the Zoeppritz equations 共given convolution does not take into account that the image has a lateral
in, e.g., Young and Braile, 1976兲, the reflection coefficients R共 ␪ 兲 at resolution dependent on the migration aperture, depth of observa-
each layer interface are calculated as a function of angle of incidence tion, and dominant wavelength 共Chen and Schuster, 1999兲. A fact of-
␪ , locally assuming a 1D layered earth. With thicknesses from the ten ignored is that the migration aperture limits the maximum dip of
initial layered reservoir model, a spiky reflectivity trace r共t兲 is built
structures that can be seen on the migration image; the steeply dip-
and convolved with a wavelet from a seismic-to-well tie to create the
ping parts are not imaged 共see, e.g., Hertweck et al., 2003 and Tox-
modeled trace.
opeus et al., 2008兲. Note that illumination and resolution on the
migration image can be analyzed efficiently using ray tracing
PROBLEMS WITH TRADITIONAL 共Lecomte, 2006兲. Additional limitations worth mentioning are the
STOCHASTIC INVERSION neglect of converted waves and anisotropy. Locally converted shear
waves influence the reflection response as a function of angle of a
For the scheme displayed in Figure 1, 1D convolution is used to
stack of thin layers 共Simmons and Backus, 1994兲. Anisotropy in the
model traces in the migrated domain 共Oldenburg et al., 1983; van
overburden can also influence this response 共Wright, 1987兲.
Riel and Berkhout, 1985兲, i.e.,
In this section, we give special attention to two further problems
s共t兲 ⳱ w共t兲 ⴱ r共t兲, 共1兲 of standard inversion: the effects of averaging reflection coefficients
over specular reflection angle ␪ and of pulse distortion on the migra-
in which the asterisk denotes temporal convolution, s共t兲 is the mod-
tion image varying with reflector dip angle ␤ . These angles are indi-
eled seismic trace, w共t兲 is the wavelet, and r共t兲 is the subsurface
spiky reflectivity — the impulse response of a 1D layered earth for cated in the upper-left portion of Figure 1 and in Figure 2, respective-
primary unconverted reflected waves. ly. In our examples, we restrict ourselves to primary unconverted
Because inversion is an iterative procedure, the forward modeling P-wave reflections in a 2.5D setting. In such a setting, point sources
must be computationally fast, justifying the choice for 1D convolu- and 3D spreading are used but the medium does not vary in the direc-
tion. However, with the increased computing power available today, tion perpendicular to the plane containing the sources and receivers
a more sophisticated forward modeler should be considered to re- 共Bleistein et al., 2001兲.

s r

Porosity,
shale
content Reflection
saturation, coeffi-
grain cients
density, (per angle)
etc.

Figure 1. Overview of conventional trace inversion with detailed forward-modeling stage.

Downloaded from [Link]

by CNOOC user
 Ray-based stochastic inversion R87

Averaging over reflection angle ␪ m0共VP共xB兲,␪ ⳱ 0,␤ ⳱ 0兲 VP共xA兲 1

n⬘0共␤ 0兲 ⳱ ⳱ .
A shortcoming of standard inversion that we suspect to be harmful m0共VP共xA兲,␪ ⳱ 0,␤ ⳱ ␤ 0兲 VP共xB兲 cos ␤ 0
for reservoir parameter estimation in a laterally strongly variable 共3兲
subsurface is that the migration image in practice is a band-limited
image of angle-averaged reflection coefficients R̄共 ␪¯ 兲, yet usually Inversion is performed usually after 1D depth-to-vertical two-way
R共 ␪ ⳱ 0兲 is assumed when constructing r共t兲 for 1D convolution 共or time conversion of the depth-migrated data. In this procedure, a scal-
when ␪ ⫽ 0, a horizontally layered model is used to compute ␪ 兲. We ing of the migration image along the vertical with local velocity oc-
do not investigate this effect in our paper 共see van der Burg 关2007兴 curs 共see equation 13兲, effectively removing the velocity dependen-
for a brief analysis兲, but our new inversion scheme properly takes cy in the stretch equation above. This yields the expression for the
into account ␪ by means of ray tracing. dip-dependent migration-induced wavelet stretch n0共 ␤ 兲 on the
depth-to-vertical time-converted zero-offset migrated image:

1
Wavelet stretch induced by reflector dip ␤ n 0共 ␤ 兲 ⳱ . 共4兲
cos ␤
On the migration image, wavelet distortion inevitably occurs
共Black et al., 1993; Brown, 1994; Tygel et al., 1994兲, whereas in the In practice, this means the wavelet representing the position of a
inversion the wavelet, w共t兲 is assumed stationary. This has an impact reflector on the migration image in vertical two-way time is
on parameter estimation in a reservoir with strong structural dip vari- stretched more with increasing reflector dip. As pointed out by Levin
ations. 共1998兲, the display of the migration image with vertical traces causes
this stretch. In our new ray-based inversion method, we ensure the
inversion occurs in a frame oriented perpendicular to the layering in
Theory the reservoir zone by incorporating reflector-oriented ray tracing.
Wavelet distortion is a consequence of varying reflection angle ␪ ,
reflector dip ␤ , and/or P-wave propagation velocity VP. The govern- Numerical test: Effect on conventional inversion
ing expression that measures wavelet stretch in a 2.5D setting along To analyze the effect of migration wavelet stretch n0共 ␤ 兲 on the in-
the vertical direction is 共Tygel et al., 1994兲 version of traces from the migration image, consider the following
simple but illustrative example: A series of synthetic data tests was
2 performed in a 2.5D normal-incidence 共 ␪ ⳱ 0兲 setting using the
m0共VP,␪ ,␤ 兲 ⳱ cos ␪ cos ␤ , 共2兲
VP model depicted in Figure 2. The experiments involve standard inver-
sion for layer properties VP and h 共⬍ ␭d, the dominant wavelength兲
with m0 denoting the ratio between the wavelet length in the two- using the wavelet derived at ␤ ⳱ 0 on a single trace from a set of ide-
way recording time domain and the length in the depth domain, al migration images shown in Figure 3, each corresponding to a flank
VP共x兲 the local P-wave velocity, ␪ 共x兲 the angle of ray incidence, and with a different dip. 共In ideal migration images, all migration arti-
␤ 共x兲 the local reflector dip. For blocky velocity models, the stretch- facts besides wavelet stretch are suppressed, described in detail lat-
evaluation point x on the depth-migrated image is chosen just above er.兲
or below the velocity discontinuity. The normal-incidence data set was characterized by a Hanning-
For a reflector on a zero-offset depth-migrated image, the stretch tapered zero-phase band-pass wavelet with corner frequencies of
induced by a reflector dip ␤ 0 at position xA, relative to the stretch at 4 – 12– 50– 75 Hz. In the experiments, exact mean values and the
point xB with zero dip, is position of ⌺ 1 共the upper boundary of the layer兲 were supplied to the
inversion as prior knowledge. Figure 4 gives the layer VP and h esti-
mates for the various dip angles. For higher dip angles, inversion re-
sults clearly deviate considerably more than two standard deviations
from the desired values.

⫺
two-way time
O

⫺
t
T

Figure 3. Migration-induced dip-dependent wavelet stretch n0共 ␤ 兲

Figure 2. Model of a structure with flank dipping at angle ␤ and a on a detail of ideal migration results for the subsurface model of Fig-
well at zero dip. The single contrasting thin layer has VP⫽2500 m/s ure 2 with dip angle as indicated. Trace separation is 10 m. Dashed
and h ⳱ 10 m. For our tests, ␪ ⳱ 0. lines denote contrasts ⌺ i.

Downloaded from [Link]

by CNOOC user
 R88 van der Burg et al.

In the next section, we present our unique inversion method, assumed to satisfy the standard ray-theoretical validity conditions
which inverts in the prestack unmigrated domain and therefore is un- 共Červený, 2001兲. The reservoir is identified by a clearly distinguish-
affected by migration stretch. able reference reflector. Instead of 1D convolution, the new method
uses 3D dynamic ray tracing in the forward-modeling kernel.
RAY-BASED STOCHASTIC INVERSION Concentrating on unconverted, primary P-wave reflections, the
To compute uncertainties in a laterally variable subsurface, specu- key vehicle for ray-based inversion is formed by a single pair of
lar reflection-angle and raypath information is needed to determine P-wave rays leaving the specular reflection point xR on the reference
the local reflection coefficient. We propose incorporating this infor- reflector at angles ␪ to the normal vector n̂共xR兲 and arriving at the
mation inside the inversion by replacing the 1D convolutional model source and receiver positions xs and xr, respectively, at the surface
with 3D dynamic ray tracing in a stochastic inversion scheme 共Fig- 共see Figure 1兲.
ure 5兲. To emphasize the evaluation along raypaths, the new inver- Layer parameters in the inversion target are updated iteratively
sion scheme is called ray-based stochastic inversion 共van der Burg et using a guided Monte Carlo algorithm. The new method minimizes
al., 2004兲, detailed in van der Burg 共2007兲. In the following, we refer the mismatch between the target reflection response—forward mod-
to it as the ray-based or new method as opposed to standard stochas- eled by 3D dynamic ray tracing to the layer interfaces in the target—
tic trace inversion, which we refer to as the conventional or old meth- and the real recordings from the prestack unmigrated data. The mod-
od. eled response is defined uniquely by the source-receiver pair 共xs,xr兲,
Earlier efforts have been made in prestack inversion for reservoir the initial directions 共 ␪ ,␾ 兲 共measured from n̂共xR兲 in the plane of
characterization using ray theory. Smith and Gidlow 共1987兲 use ray propagation at angle ␾ with the azimuth兲, the velocity model VP共x兲,
tracing through 1D models to compute reflection angles in common- and the source wavelet.
midpoint 共CMP兲 gathers, whereas in our method ray tracing can be
performed in a 3D laterally strongly varying subsurface. Also, there Simplification: Inversion along normal-incidence rays
has been a renewed effort to retrieve elastic parameters of the subsur-
If we restrict ourselves to R共 ␪ ⳱ 0兲, the inversion is performed on
face from prestack data using full waveform inversion 共Mora, 1987,
the normal-incidence gather. For the practical advantage of reducing
1989; Pratt, 1999; Pratt and Shipp, 1999兲. Unlike these methods,
data overhead, we also restrict ourselves to a 2.5D setting. In the
which concentrate on retrieving the elastic parameters on large to
overburden, 3D dynamic normal-incidence ray tracing is performed
medium scales, our approach inverts for a large range of reservoir
to the reference reflector to find those source/receiver locations of
rock and pore-fluid parameters on the reservoir scale at target level.
the survey that contain reflection information from the reservoir and
Another important difference is that the approaches mentioned
to calculate overburden amplitude effects. In the target, 1D convolu-
above are deterministic inversions that do not give the uncertainty
tion is used to model the target reflection response on the normal-in-
estimates in the results obtained with our method.
cidence gather after preprocessing to compensate for the overburden
Because we are using dynamic ray tracing, we neglect 共as in con-
amplitude effects. The flow chart is depicted in Figure 6.
ventional stochastic inversion兲 the effect of the locally converted
Transmission and spreading effects in the target are not modeled
shear wave or, alternatively, we assume its effects can be removed in
by 1D convolution. However, as long as the target interval behaves
preprocessing. Also, it is good to note that the result of the stochastic
locally as a stack of plane-parallel layers, contains a moderate num-
inversion depends on the quality of the prior information.
ber of layers with impedance contrasts that are not too large, and has
a total interval thickness of a few dominant wavelengths at most, the
Principle error will be negligible 共van der Burg, 2007兲. The 1D convolutional
For ray-based inversion, the subsurface is parameterized as an
overburden macromodel overlying a layered target reservoir and is Prestack Macro--
unmigrated velocity PSDM data
data model

Reflector
picking

Source
Normal vector PSDM data in
wavelet,
field along vertical two-
density
reflector way time
model

3D elastodynamic
ray tracing

Reservoir-
layer
parameters

Ray-based inversion

Figure 4. Estimates of VP and h from standard inversion for the dips Figure 5. Flow chart for the new ray-based inversion 共left兲 and the
indicated on the horizontal axis. Error bars denote posterior standard standard method 共right兲. The new scheme uses 3D ray-based model-
deviations. Dashed horizontal lines denote desired values for VP ing and is applied to prestack unmigrated data. PSDM⳱ prestack
共gray兲 and h 共black兲. depth migrated.

Downloaded from [Link]

by CNOOC user
 Ray-based stochastic inversion R89

共0兲
kernel offers a practical advantage because it is readily available in U3,prep 共xs ⳱ xr,xR兲 ⬇ R共xR,␪ ⳱ 0兲, 共8兲
common inversion software. The simplified method, called 1D con-
volutional ray-based stochastic inversion, is introduced in van der a necessary condition for the application of a 1D convolutional in-
Burg et al. 共2005兲. version kernel on the data set.
The depth-migrated image is used in the new method to pick the
Theory reference reflector indicating the target. Standard inversion operates
on traces from the true-amplitude prestack depth-migrated 共PSDM兲
Consider a 2.5D setting with a caustic-free isotropic-elastic sub-
image. It is constructed from the measured vertical component of the
surface consisting of an inhomogeneous overburden overlying a set
of N ⳮ n plane-parallel layers 共Figure 1兲. The isotropic point source particle velocity u̇3 ⳱ ⳵ u3 / ⳵ t as follows 共modified from Schleicher
is S共x,t兲 ⳱ w共t兲␦ 共x ⳮ xs兲, with w共t兲 denoting the real-valued wave- et al., 1993兲:

冕冕
let, assuming normalized source strength. The dip angle of the pack-
age of layers in the target is ␤ , and vertical thicknesses are hi ⬍ ␭d, 2 1 ⳵
具R共x,␪ ⳱ 0兲典 ⯝ ⳮ
with ␭d the dominant wavelength. Real thicknesses are h⬘i ⳱ hi cos ␲ xr苸A T共xr兲C0共xr兲 ⳵ t
␤ . The overburden P-velocity and density macromodels are as-
sumed to be known. ⫻ u̇3共xr;t兲t⳱tddxr,1dxr,2, 共9兲
For the chosen configuration, the vertical component of the parti-
cle-displacement vector measured from an unconverted primary in which td represents the two-way traveltime between x and xs ⳱ xr,
P-wave normal-incidence reflection at a reflection point xR on the dependent on the P-wave velocity model VP共x兲. The time derivative
lowest interface in the target ⌺ N using ray theory can be written as ⳵ / ⳵ t operating on u̇3 in the migration equation ensures the phase of
the wavelet is preserved on the migration output, i.e., it compensates
u3共xs ⳱ xr,xR;t兲 ⳱ U共0兲
3 共xs ⳱ xr,xR兲w共t ⳮ ␶ 共xs ⳱ xr,xR兲兲,
for the phase-shift effect caused by the double integration 共Newman,
共5兲 1990兲. The variable T takes into account target and overburden trans-
mission losses, C0 accounts for the free surface, and 具 典 denotes that
with U共0兲
3 denoting the real-valued amplitude function of the vertical
component of the vectorial particle displacement. The wavelet is R is spatially band limited.
placed at ␶ , the normal-incidence two-way traveltime to reflection
point xR. The expression on the right-hand side of equation 5 repre-
sents the leading term of the formal asymptotic ray-series expansion SYNTHETIC DATA TEST
solution of the general elastodynamic wave equation.
We compare the new method with conventional inversion for de-
The amplitude function U共0兲 3 can be calculated using the expres-
sion termining VP and h of the layers from a thin-layered structure that re-
sembles a real data case discussed later. With this synthetic test, we
U共0兲
3 共xs ⳱ xr,xR兲 ⳱ CB共xs,xT兲R共xR,␪ ⳱ 0兲
intend to focus on the effect of wavelet stretch; hence, we model only
noise-free primary P-wave reflections in the synthetic data.
⫻ TN共xs ⳱ xr兲关LB共xs,xT兲
Ⳮ LT共xT,xR兲兴ⳮ1, 共6兲 -
PSDM data
with xT being the intersection point of the normal-incidence ray with
interface ⌺ n, dividing the raypath in parts through the overburden
and the target. The overburden amplitude effects are denoted by CB,
which are caused by transmission through interfaces and the pres-
ence of a free surface. The expression R共 ␪ ⳱ 0兲 is the normal-inci-
dence Zoeppritz reflection coefficient measured from the ray-inci-
i⳱n T i 共 ␪ ⳱ 0兲T i 共 ␪ ⳱ 0兲 is the product
ⳮ Ⳮ
dence direction, and TN ⳱ 兿Nⳮ1
of transmission losses in the target 共while the ray pair crosses N ⳮ n
interfaces through the target: N ⬎ n and n ⱖ 1; for n ⳱ 1, the over-
burden does not contain interfaces兲. Finally, L is the relative geo-
metric spreading as defined in Červený 共2001兲. It includes the effects
of reflector curvature.
The zero-offset particle-displacement data set u3共xs ⳱ xr,xR;t兲 is
preprocessed so that all overburden amplitude effects caused by CB
and LB are removed, i.e.,

共0兲
U3,prep 共xs ⳱ xr,xR兲 ⳱ 冋 册
LB
CB
U共0兲
3 共xs ⳱ xr,xR兲

⳱ RTNLB共LB Ⳮ LT兲ⳮ1 . 共7兲

If amplitude losses within the target resulting from TN and LT are ne- Figure 6. Flow chart for 1D convolutional ray-based inversion.
glected 共i.e., TN ⬇ 1 and LT / LB Ⰶ 1兲, it follows that PSDM⳱ prestack depth migrated.

Downloaded from [Link]

by CNOOC user
 R90 van der Burg et al.

Model geometry lel layer assumption for the new method is not satisfied fully because
dips are not exactly the same along reflector normals.
The model is depicted in Figure 7. It consists of a target of five thin
layers below the sixth interface ⌺ 6 that serves as the reference reflec-
tor and the boundary between overburden and target. Layer proper- Normal-incidence data set
ties are invariant in the x2-direction, so the variations are restricted to The normal-incidence data set u̇3共xs ⳱ xr;t兲 on which 1D convo-
the 共x1,x3兲-plane and data acquisition 共using point sources and 3D lutional ray-based inversion will be applied is generated using dy-
spreading兲 is performed along the x1-direction to obtain a 2.5D con- namic ray tracing 共equations 5 and 6兲. It is computed for 601 source-
figuration. Everywhere, the layer S-wave velocity is VP / 1.7. receiver positions with 10-m spacing at the free surface. The data set
The position of the layer interfaces in the overburden including must be preprocessed before it can be inverted by the new method.
the target reference reflector ⌺ 6 is given by the following equation: The overburden amplitude losses, calculated by dynamic ray tracing
to the reference reflector, are removed using equation 7 to compen-
ⳮ共x1 ⳮ ␮兲2 sate for the effect of the laterally varying overburden on target-re-
x3,i共x1兲 ⳱ xmax
3,i ⳮ ⌬x3,i exp 共10兲
2␴ 2 flection amplitude.
For the wavelet in equation 5, the zero-phase Gabor wavelet is
for i 苸 兵1,2, . . . ,6其, with ␴ ⳱ 1000 m, ␮ ⳱ 3000 m, ⌬x3,i ⳱ x3,i max
used 共Hubral and Tygel, 1989; their equation 1兲:
ⳮ x3,i , where x3,i and x3,i can be read from Figure 7.
min max min
2
The positions of the five target interfaces below the reference re- w共t兲 ⳱ cos共2␲ f dt兲eⳮ共2␲ f dt/␥ 兲 , 共12兲
flector are calculated by applying a simple translation in depth to ⌺ 6:
with t denoting the two-way traveltime, f d ⳱ 35 Hz the dominant
frequency, and parameter ␥ ⳱ 4. Note that at zero dip in the chosen
x3,i共x1兲 ⳱ x3,iⳮ1共x1兲 Ⳮ hi ∀ i 苸 兵7,8, . . . ,11其, 共11兲 model, the two-way traveltime coincides with vertical two-way trav-
eltime. Thus, equation 12 also applies to the zero-dip wavelet
with hi denoting the 共laterally constant兲 vertical layer thicknesses 共in
present on the depth-to-time converted migrated image.
meters兲 in the target, given by hi ⳱ 兵50,8,5,10,7其 for i
苸 兵7,8,9,10,11其. This target configuration is favorable for conven-
tional inversion in the sense that dips ␤ for all target interfaces are Ideal migrated data set
the same for a given x1, which means that the associated dip-depen- The conventional method inverts traces from the PSDM 1D verti-
dent migration wavelet stretch n0共 ␤ 兲 in the target will be constant cal depth-to-time converted image of the normal-incidence data for
per trace from the migrated image. However, the locally plane-paral- layer thickness and P-wave velocity in the target. In this case, it in-
verts the ideal migrated image to isolate the effect of dip-dependent
a) migration wavelet stretch on conventional inversion; the effects of
illumination and limited lateral resolution are absent.
This ideal migration image is generated by first applying a 1D ver-
tical depth-to-time conversion to the model of the elastic properties
in depth, using the exact P-wave velocity model:
x3

tv共x1,x3兲 ⳱ 冕 2
VP共x1,x⬘3兲
dx⬘3, 共13兲
0

with tv共x1,x3兲 denoting the vertical two-way traveltime correspond-

ing to some depth location x and VP共x1,x3兲, the laterally and vertical-
ly varying P-velocity. In this way, the exact position of all interfaces
b) in vertical two-way traveltime and the size of the impedance con-
trasts associated with the interfaces are known.
Next, for every 10 m in the x1-direction, an acoustic impedance
trace is synthesized from the density and P-wave velocity models in
vertical two-way traveltime, with which the normal-incidence re-
flection coefficients are calculated. Subsequently, a spiky reflectivi-
ty trace r共tv兲 can be computed — the impulse response of a 1D lay-
ered earth, considering the earth as a linear system — with the ex-
pression 共van Riel and Berkhout, 1985; their equation 2兲
N
r共tv兲 ⳱ 兺 R j␦ 共tv ⳮ ␶ j兲, 共14兲
j⳱1

with N denoting the number of reflectors, R j the normal-incidence

Figure 7. 共a兲 Density distribution. Density in the target is constant at
2000 kg/ m3. 共b兲 Enlargement of the dashed area, showing the reflection coefficients, ␦ Dirac’s delta function, and ␶ j the lag time of
P-wave velocity distribution in the target. Overburden velocity is the jth reflector. Finally, r共tv兲 is convolved 共see equation 1兲 with a
constant at 2000 m / s. varying dip-dependent wavelet, which has the proper migration

Downloaded from [Link]

by CNOOC user
 Ray-based stochastic inversion R91

wavelet stretch n0共 ␤ 兲, depending on the local dips ␤ of the reflectors played after flattening along the reference reflector and showing
encountered. The stretch is calculated using equation 4 and applied only a portion of the bottom of the first target layer. Also, for im-
to the original Gabor wavelet at zero dip. proved display, a six-point moving-average filter is applied on the h
and VP estimates for each layer in the lateral direction. The filter
Inversion: Conventional versus ray based width is equal to the lateral resolution at target level on the migration
image.
The inversion goal is to estimate h and VP for each layer in the tar- Layer thicknesses clearly are overestimated by the old method in
get. Ray-based inversion inverts for true thickness h⬘, which after- the parts of the reservoir that have strong dips. In contrast, the thick-
ward is converted to vertical thickness h using the local dip at the re- ness estimates from the new method do not suffer from this overesti-
flection points on the top interfaces of the target layers known from mation. The VP estimates for the simplified method are closer to the
ray tracing 共van der Burg, 2007兲. Furthermore, the h and VP esti- desired values as well, although errors increase for the lower layers
mates from ray-based inversion found at irregularly distributed re- because of the neglect of spreading and transmission losses in the
flection points on each target interface are interpolated to a regular target. Posterior standard deviations 共not shown兲 do not differ much
lateral interval of 10 m, coinciding with the migration output grid, so with this noise-free data.
a direct comparison with the estimates obtained by conventional in-
version can be made. FIELD DATA TEST

Prior knowledge We tested 1D convolutional ray-based inversion and compared it

with conventional inversion using a real data set from the Gulf of
Knowledge of the reservoir before inversion for both types of in- Mexico 共courtesy Shell Offshore Inc., New Orleans, U.S.A.兲. The
versions includes the correct number of target layers; layer-␳ ; VS deepwater Gulf of Mexico field in which the inversion tests are done
⳱ VP / 1.7; and ␳ , VP, and VS outside the target. The correct wavelet is is a hydrocarbon-bearing reservoir consisting of sheet sands and
assumed to be derived from the normal-incidence section 共ray-based shales. The reservoir contains a horizontal and dipping part with dips
inversion兲 and from the horizontal part of the zero-offset migrated up to 31° 共see Figure 9兲. The test strategy is depicted in Figure 9. For
section 共conventional inversion兲. a fair comparison, both methods perform a seismic-to-well tie at the
For all layers, the prior mean ␮共VP兲 coincides with the true VP, same calibration well on a horizontal part of the target, and both use
whereas the uncertainty is described by a standard deviation of the same prior information.
␴ 共VP兲 ⳱ 250 m / s. For all h 共and h⬘兲, prior mean ␮共h兲 coincides Based on results from synthetic tests, we suspect the artifacts
with true h 共and h⬘兲 with an uncertainty described by a standard devi- caused by migration deteriorate the standard inversion results in the
ation of ␴ 共h兲 ⳱ 5 m for the first 共thick兲 target layer and ␴ 共h兲 ⳱ 1 m dipping part of the reservoir. This is verified by checking the inver-
for the other target layers. The normal distributions were bounded at sion results with the logs of an evaluation well drilled right through
a minimum value of zero; during the stochastic inversion, samples the slope after the initial inversion was complete, representing a
drawn outside this bound were rejected. The two-way traveltime to blind test. The comparison in this field data test is complicated by re-
the reference reflector was known up to a standard deviation of 2 ms. maining multiple energy, locally converted shear waves, and aniso-
tropy. None of these effects is taken into account in old or new meth-
Inversion results ods. The presence of noise and uncertainties in the prior information
The posterior reservoir models obtained from standard and ray- of the wavelet are additional complications in the comparison. For
based inversions are depicted in Figure 8. To facilitate the observa- initial test descriptions, see van der Burg et al. 共2007兲.
tion of the thickness estimates, the estimated reservoir model is dis-
Seismic data description
A 3D high-resolution seismic survey was conducted over the area.
The seismic data were recorded in an acquisition campaign that was

Figure 8. Reservoir model as found by conventional inversion 共top兲 Figure 9. The capabilities in lateral prediction of target reservoir pa-
and 1D convolutional ray-based inversion 共bottom兲 flattened along rameters away from well I 共the exploration well兲 to the dipping part
the reference reflector. Colors indicate VP misestimates. Notice the h of the target at well II are tested for conventional versus ray-based in-
overestimates for the conventional method at the dipping parts. version.

Downloaded from [Link]

by CNOOC user
 R92 van der Burg et al.

a follow-up of a previous 3D seismic survey. It was deployed to sup- depopulation was done from 80 to 48 offsets, with output offsets
port well placement, reduce uncertainty, and allow prestack inter- ranging from 450 to 6325 m and with a 125-m offset increment.
pretation. To succeed in these goals, a very large, usable signal band- The left side of Figure 11 depicts a portion of the 575-m common-
width of up to 80 Hz was needed. The desired high-frequency pres- offset gather for the dip-line selected for inversion. The reflectors
ervation and structural detail required a fine spatial sampling at the around 4000 ms on the right side of this gather mark the target for in-
surface of 6.25 m inline ⫻12.5 m crossline 共later interpolated to version. Also visible are some diffractions, mainly from the overbur-
6.25 m兲 and a temporal sampling of 2 ms. den, and a reflection-triplication caused by the convex shape of the
Because the geometry of the target was known from the previous reflectors. On these high-quality data, prestack interpretation is at-
survey, the acquisition configuration could be optimized by sailing
tainable, and prestack inversion should be feasible.
the marine acquisition vessel approximately in dip lines over the tar-
get 共see Figure 10兲.

True-amplitude PSDM
Processing before migration
To prepare the data for true-amplitude PSDM, Shell applied A common-offset true-amplitude prestack P-wave Kirchhoff
preprocessing in which the relative amplitude behavior of target in- depth migration 共Schleicher et al., 1993兲 was applied on the offset
terface reflections was preserved as much as possible. The most im- gathers from the preprocessed prestack data using the P-wave veloc-
portant process was a 3D normal moveout/dip moveout 共NMO/ ity model obtained from velocity analysis during a preceding
DMO兲—inverse NMO/DMO sequence. There were two reasons for prestack time migration and traveltime inversion. This velocity
this. First, we wanted to obtain an early structural image of the target model was defined on a grid with an inline/crossline/depth spacing
via prestack time migration before applying the computationally ex- of 100 m. The migration operator grid was sampled twice as dense-
pensive true-amplitude PSDM. Second, it suppresses acquisition ly, with a spacing of 50 m.
footprints and regularizes data. During the inverse DMO, an offset The migrated data were available in vertical two-way traveltime,
directly suitable for conventional inversion. The prestack unmi-
grated data had to be downsampled from 2 to 4 ms two-way travel-
time to obtain a manageable data size for migration; 4 ms is also the
output sampling in vertical two-way time after migration. The spa-
tial output sampling in the inline and crossline directions was 12.5 m
— twice the size of the inline midpoint distance.
After prestack migration, a stack was made for the 16 nearest off-
sets from 450 to 2325 m to increase the signal-to-noise ratio and fa-
cilitate structural interpretation. Before inversion, a phase rotation
of 90° was applied to the image, an operation sometimes applied to
Figure 10. Dip line in a 3D seismic data cube. Around this sail line, improve the interpretation of inversion results in the target. In Figure
the subsurface is assumed laterally invariant in the crossline direc- 11 共middle兲, a portion of the near-stack migrated section along the
tion. A specular ray pair is shown on the top interface of the target
with reflection angle ␪ . dip line is displayed.

Figure 11. Portion of the 575-m common-offset gather 共left兲 and the near-stack migrated image 共middle兲 along the dip line. The rightmost two
panels indicate well paths projected on the migrated image in depth, inline 共top兲, and crossline.

Downloaded from [Link]

by CNOOC user
 Ray-based stochastic inversion R93

Seismic data selection for inversion cating an inversion convergence. Notice that a five-point moving-
average filter was applied on the results obtained for the separate
The selection 共from the total data volume兲 of the migrated and
traces. In choosing the filter width, care was taken not to smooth
prestack unmigrated data to be inverted is performed on the basis of
more than the lateral resolution at target level on the migration im-
three criteria: conformity to a 2.5D setting, proximity of wells, and
age.
data quality. Although no fundamental limitations exist that prevent
applying 1D convolutional ray-based inversion in a 3D setting, the
comparative tests against conventional inversion are performed in a 1D convolutional ray-based inversion
2.5D setting to reduce data overhead by confining the analysis to the In the second part of the comparative test, the new simplified
migrated and prestack unmigrated data along a single sail line. method is used to invert the prestack near-offset section for the un-
On the right side of Figure 11, two paths of wells traversing the known reservoir-layer parameters.
target are displayed in the inline and crossline directions. The
crossline is taken at the point where well II intersects the reference Prestack seismic data selection
reflector. It shows that the inline is not placed ideally but is close
enough to wells I and II, taking into account that the target seems to To select from the 575-m common-offset gather those source-re-
remain approximately laterally constant 共2.5D兲 in the crossline di- ceiver pairs that contain reflection information on the specified in-
rection between the sail line and well II. version target, ray tracing is done to the reference reflector in the mi-
For the prestack data, the gather for the first available offset gration P-wave velocity model. In the process, the two-way travel-
共450 m兲 suffered significantly more from high-frequency noise times to the reference reflector are calculated; these are required for
around the target zone at 4000 ms than the next offset gather tying the inversion window to the traces.
共575 m兲, so the latter was preferred for ray-based inversion. The ray tracing is performed in a 2.5D setting with the source-re-
ceiver pairs on the dip line and no crossline subsurface property vari-
ations. Source-receiver distance is set to 575 m, and separation be-
Model geometry tween midpoints is 6.25 m to mimic the acquisition configuration of
The reservoir model consists of a sequence of seven layers situat- the real data. Figure 13 shows every tenth ray. The angles of inci-
ed directly above the reference horizon. It is a subset from a larger se- dence are up to ␪ ⳱ 6°.
quence of layers originally used for conventional inversion; the sub- When forward modeling the portion of the offset gather contain-
set meets the application requirements for 1D convolutional ray- ing the target reflection response, the small additional traveltimes in
based inversion. The inversion interval was chosen to extend lateral- the target caused by ␪ ⫽ 0° will be neglected, whereas the actual ray-
ly from the horizontal part at 15,100 m to the steepest part at incidence angles are taken into account when computing the reflec-
13,800 m horizontal distance, where dip increases to a maximum of tion coefficients. The reference reflector is a mildly smoothed ver-
31°. The total thickness is about 100 m, with layer rock types alter-
nating between shale and sand-shale mixture. a)
To describe the relationship between the rock properties of the
shaly sandstone reservoir rocks typically encountered in the Gulf of
Mexico and the elastic properties ␳ , VP, and VS, shale and laminated
sandstone-shale rock models are used. These models take advantage
of property trends derived from well logs 共see van der Burg, 2007,
for details兲. The reservoir-rock parameters inverted for are P-wave
velocity VP, vertical thickness h, and sand fraction SF. Here, we con-
centrate on VP and h. Gaussian distributions of these parameters are
assumed.
, , , , ,
Conventional inversion
In the first part of the comparative test, standard inversion is used
b)
to invert the stacked migrated image for the unknown reservoir-layer
parameters. A seismic-to-well match is done to derive the wavelet
from the migrated data using the detailed log information present at
well I, the first drilled well vertically penetrating the horizontal part
of the target. The top of Figure 12 shows the migrated data around
this well and the derived wavelet.
The prior mean values ␮共VP兲 and standard deviations ␴ 共VP兲 are
displayed, after flattening along the reference reflector, on the left
side of Figure 15, which summarizes the results of the old and new , , , ,
methods. Prior ␮共h兲 can be inferred from the interface positions;
␮共SF兲 is laterally constant.
Posterior ␮共VP兲, ␴ 共VP兲, and ␮共h兲 are depicted in the center of Fig- Figure 12. Good data quality on target level around well I on 共a兲 the
migrated substack and 共b兲 the 575-m common-offset gather, con-
ure 15. The seemingly rapid lateral changes in layer thicknesses re- taining more high-frequency information. Insets show 100 ms of de-
sult from the chosen way of plotting with much vertical exaggera- rived wavelets; a 90° phase rotation was applied to the migrated
tion. Posterior ␴ 共VP兲 are mostly smaller than the prior ␴ 共VP兲, indi- data.

Downloaded from [Link]

by CNOOC user
 R94 van der Burg et al.

sion from the original handpicked interface, which still fits the true the small extra traveltime in the target from having small-offset data
reflector quite well as seen from the near-stack migration image in while assuming zero offset, and the slightly different reflection coef-
the background. Also, the migration velocity model was smoothed ficients for small nonzero angles of incidence at the reflector.
in a trade-off between kinematic accuracy and dynamic stability
共van der Burg, 2007兲. Transforming the prior model from conventional inversion
Figure 13 also shows a range of source-receiver pairs that has
more than one reflection point on the reference reflector. This is the For ray-based inversion, the layer properties and thicknesses must
reflection-triplication mentioned earlier. For practical convenience, be specified along the raypath, which generally does not correspond
this area is avoided. Dealing with multivaluedness in principle can to the vertical direction along which conventional inversion is oper-
be incorporated into the forward-modeling kernel of ray-based in- ating. Figure 14 shows the situation for a portion of the reservoir
version. The chosen reflection-point range is indicated with an arrow model marked with a box in Figure 13.
along the reference reflector. The corresponding midpoint range is To obtain the prior model for the new method, a dip-dependent
12,700– 15,100 m, which includes the vertical well at 14,500 m. Fi- conversion of the prior model for conventional inversion must be
nally, the box indicates part of the reservoir 共enlargements are shown done. This conversion assumes that the target satisfies the applica-
in Figure 14兲. tion regime of 1D convolutional ray-based inversion. As a conse-
The lower panel of Figure 12 shows the area around the vertical quence of the plane-parallel layering assumption, normal-incidence
well on the 575-m common-offset gather and the wavelet derived for raypaths to the reference reflector are assumed to be straight through
ray-based inversion. For practical reasons, the derivation neglects the inversion target. The layer properties are evaluated along these
the spherical spreading and transmission losses in the reservoir zone, normal-incidence rays, starting from the reflection points on the ref-
erence reflector 共see also van der Burg, 2007; Fig-
ure 5.25兲.

Overburden amplitude correction

Dynamic ray tracing through the migration ve-
locity model to the reference horizon also yields
the laterally varying overburden effects CB / LB
needed for preprocessing the prestack unmi-
grated data in 1D convolutional ray-based inver-
sion.
The medium under the reference reflector is
chosen as a homogeneous half-space with known
elastic properties. With the medium properties
above the reflector also specified by the model,
the Zoeppritz unconverted P-wave reflection co-
efficient R共 ␪ 兲 at the reference reflector 共with ␪
⬇ 6°兲 is known exactly and can be divided out,
leaving the desired overburden amplitude effects
CB / LB in the calculated amplitudes.
Figure 13. Rays to the reference reflector in the migration velocity model; migration im- In the final amplitude correction applied to the
age is shown in the background.
traces from the common-offset gather, amplitude
variations faster than the lateral resolution on the
migration image were smoothed using an eighth-degree polynomial
fit.

Inversion results
After resampling to the grid used by standard inversion and apply-
ing a five-point moving-average filter, the VP and h estimates are
shown on the right side of Figure 15. Notice the decreased total pack-
age thickness compared to the conventional method, a result expect-
ed for the dipping part of the reservoir where the conventional meth-
od suffers from wavelet stretch. The anomalous depressions with a
peak in between, in the rightmost part of the reservoir above 14,800
m, correspond to a portion of low data quality on the offset gather.
From these observations, one should realize that the only place
Figure 14. Enlargement of the boxed area in Figure 13, showing the where a quantitative judgment of the inversion results can be made is
different evaluation directions and plotted with the prior layer posi- at the well location. This is the subject of the next section.
tions from each method.

Downloaded from [Link]

by CNOOC user
 Ray-based stochastic inversion R95

a) Prior b) Ray-based inversion c) Conventional inversion

Well II Well I Well II Well I Well II Well I
μ (VP)
Reference depth (m)

–100 2800
–80 2600
–60 2400
–40 2200
–20 2000
0 1800
200
–100
Reference depth (m)

–80 150

–60
100
–40
–20 50

0 0
14,000 14,200 14,400 14,600 14,800 15,000 14,000 14,200 14,400 14,600 14,800 15,000 14,000 14,200 14,400 14,600 14,800 15,000 σ (V )
P
Horizontal distance (m) Horizontal distance (m) Horizontal distance (m)

Figure 15. Overview of prior model 共a兲 and posteriors for ray-based 共b兲 and conventional 共c兲 inversion. Colors indicate VP 共m/s兲; mean values are
displayed on top, standard deviations below. For a clearer display, flattening is done along the reference reflector, and the vertical scale is exag-
gerated.

Comparison at well II
In the third part of the comparative test, the inversion results ob-
tained with the old and new methods are compared with the values
found at well II 共Figure 16兲. From the column of target layers, the
sandstone-shale mixture layers 共underburden 关UDB兴, 2, 4, and over-
burden 关OVB兴兲 are discerned easily using the gamma-ray and sonic
logs because of their high sand fraction. In the area around well II, it
is more difficult to discern sand-shale mixture layer 6 because of its
low sand fraction and, consequently, its low contrast with the sur-
rounding shale layers. Additional trends from the neutron log were
needed. After interpretation, the average P-wave velocity per layer
was determined from the blocked sonic log.
Figure 16 shows that the conventional method generally overesti-
mates layer thicknesses, whereas the 1D convolutional ray-based in-
version estimates are slightly better, with the values from well II
within one standard deviation from the estimated means. However,
for thin layers 2 and 4, the estimated thicknesses from the old method
are better 共but overestimated兲. The VP estimates are closer to the ac-
tual values using the new method.
Standard deviations are higher for the new method resulting from
the higher amount of noise on the offset gather compared to the near-
stack migrated section. However, the philosophy for full ray-based
inversion is to reduce these standard deviations by adding more mea-
surements 共offset gathers兲 into the inversion. Here, only one of 48
offsets was used. Moreover, an estimate with a larger standard devia-
tion does not necessarily have to be worse. For example, for the
P-wave velocity estimates of layer 2, the means are estimated about
the same. However, conventional inversion gives a misleadingly
small standard deviation; the true value falls well outside the error
bar.
The total package thickness of 86.5 m at well II is overestimated Figure 16. 共a兲 Rock column, showing true and estimated mean verti-
by the old method, as predicted by theory, to 99Ⳳ 3.5 m. The new cal thickness h of the shale 共Sh兲 and sand-shale mixture 共Ss/Sh兲 lay-
method somewhat underestimates the package thickness but re- ers at well II, with standard deviations for total h. UDB⫽under-
burden; OVB⫽overburden. 共b兲 Misestimates in h and P-wave ve-
mains within one standard deviation from the true value of 81⫾6.5 locity VP, including standard deviations 共dashed兲. Red curves are
m. conventional inversions; green curves are ray-based inversions.

Downloaded from [Link]

by CNOOC user
 R96 van der Burg et al.

CONCLUSIONS ACKNOWLEDGMENTS

Our new method for reservoir parameter estimation inverts The authors wish to thank Shell International E&P B.V. for sup-
prestack seismic reflection data before migration using stochastic in- porting this research financially and for permission to publish this
version along raypaths. The novelty in the technique is the combina- work. The field data set was made available by Shell Offshore Inc.
The reviewers are thanked kindly for their suggestions to improve
tion of ray tracing and stochastic inversion to use the original wave-
the paper. We are indebted to Paul Gary 共Shell E&P Co.兲 and Jaap
path and reflection-angle information contained in the prestack data
Leguijt 共Shell International E&P B.V.兲 for their continuous support
to estimate reservoir parameters, including uncertainties. It has be- in acquiring and inverting the field data.
come feasible to use ray tracing in a computationally intensive sto-
chastic scheme because of the increased processing power of com-
REFERENCES
puters and the computational speed and efficiency of the ray method.
The new ray-based stochastic method can be applied to estimate res- Black, J. L., K. L. Schleicher, and L. Zhang, 1993, True-amplitude imaging
ervoir parameters in a structurally complex subsurface with substan- and dip moveout: Geophysics, 58, 47–66.
Bleistein, N., J. Cohen, and J. Stockwell Jr., 2001, Mathematics of multidi-
tial lateral velocity variations and significant reflector dips. mensional seismic imaging, migration, and inversion: Springer.
Synthetic data tests show the distortion of the wavelet in the seis- Brown, R., 1994, Image quality depends on your point of view: The Leading
Edge, 13, 669–673.
mic migration image as a function of reflector dip and reflection an- Červený, V., 2001, Seismic ray theory: Cambridge University Press.
gle is an important effect that is not taken into account by conven- Chen, J., and G. T. Schuster, 1999, Resolution limits of migrated images:
Geophysics, 64, 1046–1053.
tional trace-inversion techniques. The new method operates in the Duijndam, A. J. W., 1987, Detailed Bayesian inversion of seismic data: Ph.D.
prestack unmigrated domain; therefore, it is unaffected by this mi- dissertation, Delft University of Technology.
Hertweck, T., C. Jäger, A. Goertz, and J. Schleicher, 2003, Aperture effects in
gration-induced wavelet stretch. 2.5D Kirchhoff migration: A geometrical explanation: Geophysics, 68,
The prestack data before migration that are inverted by the new 1673–1684.
Hubral, P., and M. Tygel, 1989, Analysis of the Rayleigh pulse 共short note兲:
method contain the original angle-dependent reflection information Geophysics, 54, 654–658.
needed for a good inversion for reservoir parameters. On the con- Lecomte, I., 2006, Illumination, resolution, and incidence-angle in PSDM: A
tutorial: 76th Annual International Meeting, SEG, Expanded Abstracts,
trary, conventional trace inversion techniques operate on migrated 2544–2548.
共sub兲stacks where angle-dependent reflection information is sacri- Leguijt, J., 2001, A promising approach to subsurface information integra-
tion: 63rd Annual Conference and Exhibition, EAGE, Extended Abstracts,
ficed for a better signal-to-noise ratio with respect to reflector posi- L35.
tioning. Levin, S. A., 1998, Resolution in seismic imaging: Is it all a matter of per-
spective?: Geophysics, 63, 743–749.
When applied to normal-incidence data, the new method inverts Mora, P., 1987, Nonlinear two-dimensional elastic inversion of multioffset
along raypaths that are perpendicular to the reflectors, the direction seismic data: Geophysics, 52, 1211–1228.
——–, 1989, Inversion⳱ migrationⳭ tomography: Geophysics, 54,
that offers optimal resolution for discerning the layering in the reser- 1575–1586.
voir assuming an isotropic subsurface and unconverted primary Newman, P., 1990, Amplitude and phase properties of a digital migration
process: First Break, 8, 397–403.
P-waves. Oldenburg, D. W., T. Scheuer, and S. Levy, 1983, Recovery of the acoustic
Obviously, the lateral resolution on the prestack data before mi- impedance from reflection seismograms: Geophysics, 48, 1318–1337.
Pratt, R. G., 1999, Seismic waveform inversion in the frequency domain, Part
gration is not as good as on migrated data that are focused, but the 1: Theory and verification in a physical scale model: Geophysics, 64, 888–
loss of resolution can be compensated in several ways, e.g., by accu- 901.
Pratt, R. G., and R. M. Shipp, 1999, Seismic waveform inversion in the fre-
rate modeling of edge diffractions. quency domain, Part 2: Fault delineation in sediments using crosshole
In its current implementation, the forward modeler of the new data: Geophysics, 64, 902–914.
Sambridge, M., and K. Mosegaard, 2002, Monte Carlo methods in geophysi-
method only handles isotropic subsurfaces and primary P-wave re- cal inverse problems: Reviews of Geophysics, 40, 1–29.
flections. A future study should investigate the impact on the inver- Schleicher, J., M. Tygel, and P. Hubral, 1993, 3-D true-amplitude finite-off-
set migration: Geophysics, 58, 1112–1126.
sion results of anisotropy in the overburden and converted shear Simmons, J. L., and M. M. Backus, 1994, AVO modeling and the locally con-
waves. A more robust extension of the method could then include verted shear wave: Geophysics, 59, 1237–1248.
Smith, G. C., and P. M. Gidlow, 1987, Weighted stacking for rock property
such effects in forward modeling. estimation and detection of gas: Geophysical Prospecting, 35, 993–1014.
Results from the Gulf of Mexico field data study show that the Tarantola, A., 2005, Inverse problem theory and methods for model parame-
ter estimation: Society for Industrial and Applied Mathematics.
thickness and P-wave velocity estimates obtained with the simpli- Toxopeus, G., J. Thorbecke, K. Wapenaar, S. Petersen, E. Slob, and J.
fied new method are better than those obtained with conventional in- Fokkema, 2008, Simulating migrated and inverted seismic data by filter-
ing a geologic model: Geophysics, 73, no. 2, T1–T10.
version, in the sense that the new estimates are more often within one Tygel, M., J. Schleicher, and P. Hubral, 1994, Pulse distortion in depth migra-
standard deviation from the desired values. Hence, a good result can tion: Geophysics, 59, 1561–1569.
van der Burg, D. W., 2007, Ray-based stochastic inversion of pre-stack seis-
be obtained by working properly with a tiny fraction of the data. mic data for improved reservoir characterisation: Ph.D. dissertation, Delft
Only 2% of the available prestack data were used with 1D convolu- University of Technology.
van der Burg, D., A. Verdel, and K. Wapenaar, 2004, Ray-based stochastic in-
tional ray-based inversion. However, we expect that adding multi- version: 74th Annual International Meeting, SEG, Expanded Abstracts,
offset data would constrain the inversion operation better. 1607–1610.
——–, 2005, Ray-based stochastic inversion for reservoir parameters using
It would be interesting to investigate the potential of reusing rays 1D convolutional forward modeling: 75th Annual International Meeting,
calculated on the diffraction grid for preserved-amplitude Kirchhoff SEG, Expanded Abstracts, 1441–1444.
——–, 2007, Ray-based stochastic inversion: Improving reservoir character-
migration in ray-based stochastic inversion to save computing time ization in a Gulf of Mexico field: 77th Annual International Meeting, SEG,
and to interweave inversion with migration. Expanded Abstracts, 1775–1779.
van Riel, P., and A. J. Berkhout, 1985, Resolution in seismic trace inversion

Downloaded from [Link]

by CNOOC user
 Ray-based stochastic inversion R97

by parameter estimation: Geophysics, 50, 1440–1455. Wright, J., 1987, The effects of transverse isotropy on reflection amplitude
Veeken, P. C. H., and M. Da Silva, 2004, Seismic inversion methods and versus offset 共short note兲: Geophysics, 52, 564–567.
some of their constraints: First Break, 22, 47–70. Young, G. B., and L. W. Braile, 1976, Acomputer program for the application
White, R., and R. Simm, 2003, Good practice in well ties: First Break, 21, 75– of Zoeppritz’s amplitude equations and Knott’s energy equations: Bulletin
83. of the Seismological Society of America, 66, 1881–1885.

Downloaded from [Link]

by CNOOC user