GEOPHYSICS, VOL. 64, NO. 1 (JANUARY-FEBRUARY 1999); P. 222–229, 12 FIGS.
3-D preserved amplitude prestack depth migration
on a workstation
Philippe Thierry∗ , Gilles Lambaré∗ , Pascal Podvin∗ ,
and Mark S. Noble∗
is incorporated in this kinematic approach, migration-based
amplitude variation with offset (AVO) analysis cannot be
performed (Beydoun et al., 1993; Tura et al., 1997). Moreover,
complex structure imaging is often unsuccessful when based
only on first-arrival traveltimes (Geoltrain and Brac, 1993). Until now, the extension of eikonal solvers to amplitude and multipath computation have not been successful. On the contrary,
ray tracing does not suffer such limitations. It is now possible to
improve 3-D depth migration since the last generation of dynamic ray-tracing algorithms based on the wavefront construction (WFC) method (Vinje et al., 1993, 1996a, b; Lambaré et al.,
1996; Lucio et al., 1996) is perfectly adapted to this purpose.
In addition to computational efficiency, the possibility of obtaining quantitative information on reflectors (namely, reflection coefficient or impedance contrast) is another major advantage of ray-based migration methods. Theoretically, migration
is reformulated in the frame of linearized seismic inversion.
In particular, it was shown that ray+Born (Beylkin, 1985) or
ray+Kirchhoff (Bleistein, 1987) linearized inversion translated
into rather simple quantitative migration formulas. Jin et al.
(1992) reconciled such direct inversion approaches with the
general stochastic approach (Tarantola, 1987; Beydoun and
Mendes, 1989) and introduced an iterative minimization together with extension to multiparameter inversion (Forgues
and Lambaré, 1997). Here, we define preserved amplitude
prestack depth migration (PAPsDM) as the result of a single iteration of ray+Born inversion for a single parameter (Thierry
and Lambaré, 1995). We already demonstrated that PAPsDM
could reasonably recover the amplitude of reflectors (Thierry
et al., 1996; Operto et al., 1997a, b).
We present an algorithm for 3-D PAPsDM together with an
application to real data. We show that 3-D PAPsDM can be
done with a relatively low computational cost (all results presented in this paper were computed on a workstation). Such a
performance is attained with the use of embedded interpolations that were discussed and calibrated for the 2-D and 3-D
cases in Operto et al. (1997a, b) and Thierry et al. (1999).
ABSTRACT
We present an algorithm based on the ray+Born approximation for 3-D preserved amplitude prestack depth
migration (PAPsDM) of seismic reflection data. This
ray+Born inversion scheme allows the quantitative recovery of model perturbations. The Green’s functions
are estimated by dynamic ray tracing in 3-D heterogeneous smooth velocity fields with a wavefront construction (WFC) method. The PAPsDM algorithm was implemented on a single-processor Sun SPARC 20 workstation.
Special attention was paid to CPU efficiency and memory requirements. We present an application on a 3-D
real marine data set (13 Gbytes). About one week of
CPU time is needed to obtain a migrated image of 7 ×
1 × 1 km.
INTRODUCTION
Three-dimensional prestack depth migration (PsDM) is certainly the most accurate approach for imaging laterally heterogeneous media. Because of extensive CPU and memory
requirements for 3-D seismic applications, most of the present
3-D PsDM algorithms are based on ray theory (Kirchhoff migration). In fact, it is now accepted that ray theory provides an
excellent compromise between precision and computational
efficiency in 3-D heterogeneous media.
The main difficulty of ray-based migration is the computation of traveltimes in the target for all shot and receiver
positions. In this perspective, a finite-difference solution of
the eikonal equation for computing first-arrival traveltimes
(Vidale, 1988; Podvin and Lecomte, 1991) was an early significant breakthrough. It provides extremely fast algorithms for
kinematic migration (Reshef, 1991; Mufti et al., 1996; Noble
et al., 1996). Yet, since no reliable amplitude information
Presented at the 66th Annual International Meeting, Society of Exploration Geophysicists. Published on Geophysics Online, October 8, 1998.
Manuscript received by the Editor August 4, 1997; revised manuscript received April 9, 1998.
∗
Ecoles des Mines de Paris, ARMINES-GEOPHY, 35 rue Saint Honore, 77305 Fontainebleau Cedex, France. E-mail: philippe.thierry@geophy.
ensmp.fr, gilles.lambare@geophy.ensmp.fr, pascal.podvin@geophy.ensmp.fr, mark.noble@geophy.ensmp.fr.
°
c 1999 Society of Exploration Geophysicists. All rights reserved.
222
3-D PAPsDM
In this paper, we first recall the preserved amplitude migration approach in the specific case of a 3-D marine acquisition. Practical aspects of the PAPsDM algorithm, including
interpolation strategies, are then detailed in the context of our
application. Three million seismic traces were migrated into a
small target (7 × 1 × 1 km). All parameters, including computation times, are provided for future comparisons and improvements. We show a residual data section, namely, the difference
between an observed shot gather and a synthetic one computed with 3-D ray+Born modeling. Finally, we show that 3-D
PAPsDM improves the image significantly when compared to
2.5-D PAPsDM.
THEORETICAL ASPECTS
Preserved amplitude migration is based on the ray+Born
approximation (Cohen et al., 1986). In this approximation, the
model is split into a reference model (background) and a perturbation model δm. The perturbation of Green’s function (δG)
corresponding to the reflected scalar wavefield is related linearly to the model perturbation by
δGcal (r, ω, s) = ω2
Z
C(δm) =
1
2 LSR
The solution δm minimizing the cost function (2) is well
known as
δm = (B † QB)−1 B † QδGobs ,
dωQ|δGobs − δGcal (δm)|2 ,
(2)
ω
where L , S, R are, respectively, the line, shot, and receiver numbers, and Q is a weighting function equivalent to a covariance
matrix in the data space.
(3)
where † denotes a transposed conjugated operator. The operator (B † QB) is called the Hessian. In classical seismic inversion (Tarantola, 1984; Beydoun and Mendes, 1989), this operator is discretized and becomes a huge matrix that cannot
be inverted numerically. Consequently, iterative gradient minimization generally is introduced. To improve this local minimization, Jin et al. (1992) proposed to diagonalize the Hessian by choosing a judicious local weighting function Q in the
cost function (2). Such an approximate Hessian allows for an
efficient quasi-Newton iterative minimization of the cost function (Lambaré et al., 1992; Thierry and Lambaré, 1995). The
first iteration is equivalent to the Beylkin migration formula
(Beylkin, 1985). For our marine case, the final 3-D PAPsDM
expression is
δm(x0 ) =
dx δm(x)A(r, x, s) eiωT (r,x,s) , (1)
where s, r, and x respectively denote the source, receiver,
and diffractor point positions, ω is the angular frequency, and
A(r, x, s) = A(r, x)A(x, s) and T (r, x, s) = T (r, x) + T (x, s)
are, respectively, the amplitude of the Born operator and the
two-way traveltime computed in the reference model (Figure 1).
The linear relation (1) can be inverted in the frame of the
stochastic inverse problem theory (Tarantola, 1987). Consider
a standard 3-D marine data set (δGobs ) with (ideally) parallel
acquisition lines where each line consists of a number of shots
recorded with moving streamers (Figure 1). Let us introduce
the weighted ℓ2 cost function C associated with a model perturbation δm(x),
XZ
223
XXX
1
max
[|q(x0 )|]min L S R
× E(r, x0 , s)δGobs (r, T (r, x0 , s), s),
(4)
with
¯
¯
1 |q|1L1S1R ¯¯ ∂(q) ¯¯
E(r, x0 , s) =
,
(2π )2 A(r, x0 , s) ¯ ∂(L , S, R) ¯
(5)
where 1L, 1S, and 1R denote, respectively, the increments in
line, shot, and receiver numbers, and vector q = pr + ps is the
sum of the slowness vectors (Figure 1).
The migration operator E compensates for the two-way amplitude (1/A), accounts for the acquisition geometry through
the Jacobian |∂(q)/∂(L , S, R)|. The Jacobian operator can be
written as
¯ ¯
¯
¯
¯ ∂(q) ¯ ¯ ∂(pr + ps ) ∂(r, s) ¯
¯=¯
¯
¯
¯ ∂(L , S, R) ¯ ¯ ∂(r, s) ∂(L , S, R) ¯,
(6)
where the paraxial quantities (∂(pr )/∂(r), ∂(ps )/∂(s)) are computed during ray tracing (Farra and Madariaga, 1987; Lambaré
et al., 1996) and the matrix (∂(r, s)/∂(L , S, R)) is inferred from
the acquisition geometry. It can be evaluated by finite differencing r(L, S , R) and s(L, S , R) (Clochard et al., 1997) or more
simply approximated by constant mean values, as was done in
our application.
PRACTICAL ASPECTS
FIG. 1. Ray+Born approximation. Vector q is the sum of the
slowness vectors pr = ∇x0 T (x0 , r) and ps = ∇x0 T (x0 , s).
Because of the size of a 3-D data set and target, 3-D prestack
depth migration is one of the most expensive processing steps.
Practical implementation requires the use of simplifications.
As discussed in Thierry and Lambaré (1995), we use embedded interpolations to limit the number of rays traced from the
surface and the number of operations per image point in the
target.
Accounting for the preserved amplitude term, our 3-D
PAPsDM algorithm [equation (4)] requires 11 parameters [the
traveltime T , the amplitude A, the slowness vector p, the
Jacobian operator ∂(p)/∂(x), where x denotes the surface position of the ray shooting point]. The computation of these 11 parameters is done in a smooth background velocity field with the
dynamic ray-tracing algorithm developed by Lucio et al. (1996)
224
Thierry et al.
and based on the WFC method proposed by Vinje et al. (1993).
This algorithm uses a uniform ray-density criterion, providing
an optimal compromise between computational efficiency, accuracy, and robustness, even in the case of multipathing and
caustics.
When the reference velocity field does not exhibit very
strong lateral variations, surface interpolation of Green’s functions may be reasonably proposed to reduce the number of
maps computed by ray tracing. We chose to compute all the
maps during a preliminary stage and to store them on disk
(Thierry and Lambaré, 1995). The required maps are read during the migration and partially stored in memory to take advantage of the vicinity of receiver positions to limit the number
of input/output operations.
Computation of the migration operator is an expensive operation. If T and E are supposed to have slow variations over
the target zone, it appears reasonable to calculate them over
a sparse regular target grid and then use linear interpolation
onto the fine target grid. Finally, the contribution of each trace
to the final migrated image is limited to a region around the
common-midpoint position (Lambaré et al., 1992; Thierry and
Lambaré, 1995).
The precision of all these approximations was discussed and
calibrated in structurally simple cases (as in our 3-D application) and also in case of complex synthetic models (Operto
et al., 1997a; Thierry et al., 1999).
Preprocessing and data management
Data preprocessing (including spiking deconvolution, waterbottom multiple removal, and filtering) and velocity analysis
were done by Elf GRC.
In order to deal with limited disk resources and to manage
this 13 Gbyte data set, we developed software for a very simple
FIG. 3. The 29 navigation lines of the 3-D real data set. Note
the irregularity in the acquisition.
APPLICATION
We present an application of the 3-D PAPsDM scheme to a
real 3-D marine data set acquired in the North Sea and provided by Norsk Hydro. This data set consists of 29 survey lines
with about 500 shots per line (see the common first offset section in Figure 2). The acquisition consisted of one source and
two streamers. Each streamer had 120 receivers spaced 25 m
apart (about 3 million traces). The mean distance between cables was 100 m. Figure 3 shows the distribution of shot positions
for the 29 navigation lines and Figure 4 shows receiver positions for the first shot of lines 25 and 26. This real acquisition
geometry, substantially different from an ideally periodic pattern, is accounted for in the migration scheme.
FIG. 4. Streamers positions for the first shot of the lines 25 and
26. The feathering goes up to 500 m.
FIG. 2. Common first offset section for line 16 after preserved amplitude preprocessing. In accordance with our
approach, no trace equalization was performed.
3-D PAPsDM
data compression/decompression algorithm. Our goal was to
store data in a moderately compressed form (a compression
rate on the order of three was sufficient) that would still allow
fast, flexible, and direct access from the migration/inversion
codes. A compressed seismic data format was adapted from
the standard IEEE-754 32 bits floating point encoding norm,
simply using fewer bits for both mantissa and exponent. For our
data set, we adopted a (very conservative) encoding on 15 bits
(exponent on 5 bits, mantissa on 9 bits), ensuring a relative
precision better than 0.1%, whereas the dynamic range was
on the order of 4.0e + 9. Last but not least, simply accounting
for the applied mute yielded an almost supplementary 30%
gain in file size. As a result, each seismic line originally occupying about 500 Mbytes was compressed onto a single file
roughly one-third times smaller. Both computational and precision penalties caused by the decompression filter was found
to be insignificant.
225
done in a coarsely sampled moving cube (5.0 × 2.6 × 2.0 km)
centered around each shotpoint, with its top at 1490 m and
sampled at 101 m in the three directions.
Computation time for each ray tracing was 290 s on a Sun
SPARC 20 workstation. These initial maps (317 Mbytes) were
stored on disk, and their computation took a total time of
20 hours. The second step consisted of a linear interpolation
of the initial 500 × 500 m surface grid onto one at 100 × 100 m
(Figure 7a). Interpolations were performed at constant offset
to be sensitive only to lateral variations of the velocity field.
Finally ray parameter maps were interpolated by cardinal cubic B-splines onto the desired target in about 4 hours
Velocity model and ray tracing
A velocity model was provided in the form of a velocity grid
with 500 m in-line (i.e., in the east-west direction), 25 m crossline (in the north-south direction), and 200 m vertically. In the
ray-tracing code used for 3-D PAPsDM, model properties must
be represented by cardinal cubic B-spline weights (Lucio et al.,
1996). The original model was converted into B-spline representation by fitting the slowness field (Operto et al., 1997b). As
there were only gentle cross-line velocity variations, the spacing
of the B-spline nodes was chosen to be 500 m in both in-line
and cross-line directions, while the original vertical spacing was
kept unchanged (Figure 5).
Because of the gentle structure of the macro model, we decided to compute the ray parameters maps using three embedded interpolation steps. First, shooting was performed from a
coarse surface grid (Figure 6) with 36 × 7 points and a spacing
of 500 m in both the x- and y-directions. Each ray shooting was
FIG. 6. Initial ray shotpoint grid. The spacing of this surface
grid was 500 m in both the x- and y-directions, and there were
36 × 7 points.
FIG. 5. Several slices in the 3-D heterogeneous velocity model. This very simple model was provided by Elf GRC
after conventional velocity analysis.
226
Thierry et al.
(Figure 7b). The advantage of this three-step procedure is that
once the initial ray tracing has been done, ray parameter maps
can be interpolated easily and rapidly onto any desired target.
a)
Results
b)
FIG. 7. (a) Linear interpolation of maps from a 500 × 500 m
surface grid onto a 100 × 100 m surface grid. Interpolations
were performed at constant offset in order to be sensitive only
to lateral variations of the velocity field. (b) The ray parameter
maps were interpolated by cubic B-splines onto the desired
coarsely sampled target.
As a first application, we imaged a 7 × 1 × 1 km volume
using all the 29 navigation lines. The bin size was 25 × 25 × 5 m.
This application took one week on a Sun SPARC 20 workstation
using only 48 Mbytes of RAM.
Figure 8 shows a vertical slice where tilted blocks appear
separated by faults. Figures 9 and 10 allow us to see the lateral
variations of the structure making it easier to understand the
3-D geometry of the target.
From the migrated image, we can compute synthetic traces
using the 3-D ray+Born expression (1). Figure 11 shows the
real, the calculated, and the residual data for shot 201 of navigation line 16. In calculated data, major reflected events were
recovered with a slight time mismatch because an imperfect
velocity model was used. Thus this imperfection essentially explains the high energy in residual data.
As a second application, we imaged a vertical in-line section
using seven seismic lines (numbers 22–28) with the 3-D algorithm, and only one line (number 25) with the 2.5-D algorithm.
Figure 12 show the 3-D and 2.5-D migrated sections of the
velocity perturbation. Computation times on a Sun SPARC 20
workstation were, respectively, 4 hours and 12 minutes, including the Green’s function computation. Profiles are plotted at
the same scale. Three-dimensional PAPsDM greatly improved
the general quality of the image. The in-plane hypothesis fails
both because of the strong lateral drift of the streamers (Figure 4) and the lateral variation of the geological structure (the
major faults are not oriented orthogonally to the lines).
CONCLUSION
We have demonstrated that 3-D PAPsDM for a small target
volume is affordable even on standard workstations if some
FIG. 8. Vertical slice in the 3-D migrated cube. Some majors faults appear in the tilted block area.
3-D PAPsDM
227
FIG. 9. Horizontal section in the 3-D migrated cube, clearly showing the major faults in the tilted block area.
FIG. 10. Sections in the 3-D migrated cube showing the lateral variation of the structure.
judicious simplifications are introduced. We have shown that
application to 3-D real data greatly improves the 2.5-D results
even in case of small variations in the velocity model.
The interpolation effects as well as the quantitative recovery
of model perturbations was detailed in Thierry and Lambaré
(1995) and Thierry et al. (1999). The main advantages of
interpolations are in the significant decrease of the required
memory space and in a slight decrease of the computation
time.
The model is certainly not very complex in terms of velocity field. However, it has been demonstrated that this method
could be applied successfully to complex cases (Marmousi 2-D
and Overthrust 3-D) (Operto et al., 1997a, b).
With such an approach, computation time for PAPsDM
is slightly greater than the one needed for kinematic PsDM
(which cannot handle any reliable amplitude information), but
it introduces a practical method for 3-D migration-based AVO
analysis (Tura et al., 1997) since it can be transformed easily
into ray+Kirchhoff common-offset migration.
This efficient algorithm could also be used in the near future
as a basis for migration velocity analysis (Liu and Bleistein,
1995; Tieman, 1995).
ACKNOWLEDGMENTS
This work was partly funded by the European Commission and Norwegian Research Council in the framework of
the JOULE projects “3-D Asymptotic Seismic Imaging” (contract JOU2-CT93-0321) and “Reservoir Oriented Delineation
Technology” (contract JOF3-CT95-0019). We thank Norsk
Hydro for providing the 3-D data set, and C. Hanitzsch and
A. Tura (Elf GRC) for preprocessing it. The paper benefited
from very fruitful discussions with Stephane Operto (École des
Mines de Paris), Jan Pajchel (Norsk Hydro), and Henri Calandra (Elf Aquitaine).
REFERENCES
Beydoun, W. B., and Mendes, M., 1989, Elastic ray-born ℓ2 -migration/inversion: Geophys. J., 97, 151–160.
228
Thierry et al.
FIG. 11. Real data (a), calculated data (b), and residuals (c)
for one streamer of shot gather 201. Considering the low signal/noise ratio of data and the inacurracy of the velocity model
used for migration, a low-pass filter ([0, 14, 20, 30] Hz) had to
be applied to the data to allow for significant comparisons.
Beydoun, W. B., Hanitzsch, C., and Jin, S., 1993, Why migrate before
AVO? A simple example: 55th Ann. Mtg., Eur. Assn. Expl. Geophys., Extended Abstracts, B044.
Beylkin, G., 1985, Imaging of discontinuities in the inverse scaterring
problem by inversion of a causal generalized radon transform: J.
Math. Phys., 26, 99–108.
Bleistein, N., 1987, On the imaging of reflectors in the earth: Geophysics, 52, 931–942.
Clochard, V., Nicoletis, L., Svay-Lucas, J., Mendes, M., and Anjos, L.,
1997, Interest of ray Born modeling and imaging for 3-D walka-
ways: 59th Ann. Mtg., Eur. Assn. Geoscientists Engin., Expanded
Abstracts, P035.
Cohen, J. K., Hagin, F., and Bleistein, N., 1986, Three-dimensional Born
inversion with an arbitrary reference: Geophysics, 51, 1552–1558.
Farra, V., and Madariaga, R., 1987, Seismic waveform modeling in
heterogeneous media by ray perturbation theory: J. Geophys. Res.,
92, 2697–2712.
Forgues, E., and Lambaré, G., 1997, Parametrization study for acoustic
and elastic ray+born inversion.: J. Seis. Expl., 6, 253–277.
Geoltrain, J., and Brac, S., 1993, Can we image complex structures with
first-arrival traveltime?: Geophysics, 58, 564–575.
Jin, S., Madariaga, R., Virieux, J., and Lambaré, G., 1992, Twodimensional asymptotic iterative elastic inversion: Geophys. J. Internat., 108, 575–588.
Lambaré, G., Lucio, P. S., and Hanyga, A., 1996, Two-dimensional multivalued traveltime and amplitude maps by uniform sampling of ray
field: Geophys. J. Internat., 125, 584–598.
Lambaré, G., Virieux, J., Madariaga, R., and Jin, S., 1992, Iterative
asymptotic inversion of seismic profiles in the acoustic approximation: Geophysics, 57, 1138–1154.
Liu, Z., and Bleistein, N., 1995, Migration velocity analysis: Theory and
an iterative algorithm: Geophysics, 60, 142–153.
Lucio, P. S., Lambaré, G., and Hanyga, A., 1996, 3-D multivalued travel
time and amplitude maps: Pageoph, 148, 449–479.
Mufti, I. R., Pita, J. A., and Huntley, R. W., 1996, Finite-difference
depth migration of exploration-scale 3-D seismic data: Geophysics,
61, 776–794.
Noble, M., Marsset, B., Missiaen, T., and Versteeg, W., 1996, Near surface 2-D and 3-D data processing—beyond stack: 58th Ann. Mtg.,
Eur. Assn. Expl. Geophys, Expanded Abstracts, M033.
Operto, S., Lambaré, G., Podvin, P., and Thierry, P., 1997a, 3-D
preserved amplitude prestack imaging of the Overthrust model:
59th Ann. Mtg., Eur. Assn. Geoscientists Eng., Expanded Abstracts,
A043.
——— 1997b, CPU efficient ray+Born inversion for complex velocity
fields: 59th Ann. Mtg., Eur. Assn. Geoscientists Eng., Expanded
Abstracts, P033.
Podvin, P., and Lecomte, I., 1991, Finite difference computation of
traveltimes in very contrasted velocity model: A massively parallel
approach and its associated tools: Geophys. J. Internat., 105, 271–284.
Reshef, M., 1991, Prestack depth imaging of three-dimensional shot
gathers: Geophysics, 56, 1158–1163.
Tarantola, A., 1984, Inversion of seismic reflexion data in the acoustic
approximation: Geophysics, 49, 1259–1266.
——— 1987, Inverse problem theory: Methods for data fitting and
model parameter estimation: Elsevier.
Thierry, P., and Lambaré, G., 1995, 2.5D true amplitude migration
on a workstation: 65th Ann. Mtg., Soc. Expl. Geophys., Expanded
Abstracts, 156–159.
Thierry, P., Lambaré, G., Podvin, P., and Noble, M., 1996, 3D prestack
preserved amplitude migration: Application to real data: 66th Ann.
Mtg., Soc. Expl. Geophys., Expanded Abstracts, 555–558.
Thierry, P., Operto, S., and Lambaré, G., 1999, Fast 2-D ray+Born
migration/inversion in complex media: Geophysics, 64, 162–181.
Tieman, H. J., 1995, Migration velocity analysis: Accounting for the
effects of lateral velocity variations: Geophysics, 60, 164–175.
Tura, A., Hanitzsch, C., and Calandra, H., 1997, 3-D AVO migration/inversion of field data: J. Seis. Expl., 6, 117–125.
Vidale, D., 1988, Finite-difference calculation of travel time: Bull. Seis.
Soc. Am., 78, 2062–2076.
Vinje, V., Iversen, E., Astebøl, K., and Gjøystdal, H., 1996a, Estimation of multivalued arrivals in 3-D models using wavefront
construction—Part I: Geophys. Prosp., 44, 819–842.
——— 1996b, Estimation of multivalued arrivals in 3-D models
using wavefront construction—Part II: Tracing and interpolation:
Geophys. Prosp., 44, 843–858.
Vinje, V., Iversen, E., and Gjøystdal, H., 1993, Traveltime and amplitude estimation using wavefront construction: Geophysics, 58,
1157–1166.
3-D PAPsDM
FIG. 12. (a) 2.5-D PAPsDM section of line 25 (y = 5.81 km). (b) 3-D PAPsDM result for the same section using
seven navigation lines (numbers 22–28). The two sections are plotted with the same gain.
229