Reflection tomography versus stacking velocity analysis

£1 LIEI I EIJFII 51C5 ELSEVIER Journal of Applied Geophysics35 (1996) 1- 13 Reflection tomography versus stacking velocity analysis Gualtiero Boehm *, Jos~ M. Carcione, Aldo Vesnaver Osserz'atorio Geofisico Sperimentale, P. O.Box 2011, 34016 Trieste, Italy Received 6 March 1995;accepted 5 September 1995 Abstract The travel-time inversion of reflected arrivals reconstructs the structure of the main interfaces with a precision comparable to the seismic wavelength. The resolution of the conventional stacking velocity analysis is lower, i.e. of the order of the seismic spread length. Furthermore, the stacking velocity field is defined in the time domain, and its conversion to the depth domain is not straightforward. Both methods require selecting: this is quite difficult and time consuming for the pre-stack reflected events, but simpler and inexpensive for the velocity spectra. There is thus a tradeoff between the two approaches in terms of costs and benefits. In this paper we compare the main features of the two methods by applying them to different synthetic models of increasing complexity. We modelled the related seismograms using the Fourier pseudo-spectral method. 1. Introduction Stacking velocity analysis is a fundamental procedure of digital seismic processing. Exploration geophysicists already exploited its basic principle (the multifoid coverage) in the era of analogue recording (Mayne, 1962, 1989). Velocity spectra on the other hand have been used for over two decades since their introduction (Taner and Koehler, 1969; Neideil and Taner, 1971; Sguazzero and Vesnaver, 1987). Initially, this method was the best tool for recovering the velocity distribution in the time domain and, indirectly, also in the depth domain. In current practical applications, it is the fastest and cheapest approach to building a macro-model of the velocity * Corresponding author. field. This estimate often serves as an initial approximation for more sophisticated inversion algorithms. The introduction of reflection tomography in seismic exploration is more recent. Bishop et al. (1985) tried a simultaneous inversion of the velocity field and the depth of reflectors; they pointed out the limited reliability in estimating the vertical velocity variations between the reflectors. Carrion et al. (1993a,b)) estimated separately the velocity field and the depth of the reflectors, detecting only the lateral variations of the velocity within each layer. This approach is robust and reliable, and we therefore adopt it in this paper. We compare these two methods in the following sections, trying to help the seismic analyst in his everyday dilemma: is it better to aim for a stable, quick and cheap result, or to carry out expensive and time-consuming work that may give the highest resolution but perhaps not the highest reliability? 0926-9851/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. SSDI 0926-9851(95)00025-9 2 G. Boehm et al./Journal ~{'Applied Geophysics 35 (1996) 1-13 2. Stacking velocity analysis and with the dip move-out correction, respectively) are other examples of provelocities. We will see that in a certain sense the velocity field obtained by reflection tomography is also a provelocity field. A detailed description of the main features of stacking velocity analysis may be found in the classic paper of Al-Chalabi (1974) and in several textbooks (e.g., Hubral and Krey, 1980; Hatton et al., 1986; Yilmaz, 1987). Here, we would like to briefly recall some assumptions and characteristics that may explain the resolution limits of the velocity spectra. Stacking velocity functions are a relationship between two parameters, the zero-offset travel time and the stacking velocity. These functions correspond to real physical conditions only if the medium is homogeneous and the reflector is parallel to the surface where the sources and receivers are located; otherwise, they lose any physical meaning and only approximate in a compact form the actual travel times of waves at different offsets. So, when local inhomogeneities and dipping reflectors are present, we cannot directly recover the reflector depth and the velocity in the medium from a few travel times, e.g. by applying the Dix formula (Dix, 1955). Loinger (1983) showed that the stacking velocity depends linearly on the lateral velocity variations, instead of just on the local velocity. Therefore, when velocity anomalies occur, it is necessary to apply some spatial filter to obtain mainly the correct low wave numbers (Lynn and Claerbout, 1982), or some additional tool, such as a tomographic analysis (Harlan, 1989). Stacking velocity is probably the best example of what Al-Chalabi (1994) calls a provelocity: a parameter of the processing sequence with the physical dimensions of a velocity, which assumes values related to the actual local velocity of seismic waves. Pre-stack and post-stack migration velocities (without 1p 1 2 3. Reflection tomography Reflection tomography is a much more recent and promising method for reconstructing the velocity field at depth, which should overcome most of the limits of the stacking velocity approach. First, since each reflected arrival is considered singly, no lowpass spatial filter is applied, as implied for the stacking velocity which is a moving average with a window length equal to the spread. Secondly, since each event is modelled separately by ray tracing, local inhomogeneities crossed by a single ray will influence the inversion result; this is not so for velocity spectra, where only some of them will be visible if the time shift they produce can significantly modify the moving average along the spread. A further important consequence of the individual nature of reflected arrivals in seismic tomography is that quite complex geometries, such as dipping horizons or sharp lateral discontinuities, do not violate any of the assumptions of the algorithm. In such cases reflection tomography will yield good results, unlike velocity spectra, as will be shown later. In their pioneering work, Bishop et al. (1985) proved that it is barely possible to recover the vertical velocity gradients between reflectors. Therefore, in the model parametrization adopted here (see also Carrion et al., 1993a,b), we choose pixels such that their base and top coincide with adjacent reflecting 30 DISTANCE(Km) ~3 ,~15 4 6 7 O 8 10 11 12 0 0 + Ix_ ~--~2- .I. "I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ,t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [ l I [ 1 . . . . ~ 2 . . . . I 3 . . . . I 4. . . . . I 5 . . . . I " A I § . . . . I 7 ' l ) F I 8 . . . . ] 9 . . . . 110 . . . . I 11 2 . . . . 12 Fig. I. Synthetic model. The circles at the top indicate the position of the sources. The zone outside the dotted line was introduced to eliminate wrap-around from the grid boundaries. G. Boehm et al. / Journal of Applied Geophysics 35 (1996) 1-13 3 estimating in turn the depth of the reflection points and the velocity field. The initial guess for the velocity distribution and reflector shapes can be quite different from the actual one. The initial guess for the reflector depth can be improved (as will be seen surfaces and with vertical sides. Then we force the inversion algorithm to estimate an average velocity between the interfaces and the lateral changes along the layers. The inversion procedure advances iteratively, by DISTANCE (km) 1 g 9 ..!..L..!...'_L.L.'..L !. !..'. ! . ' ' ::::i ~: :~ z.~ ~. 5 4 5 ! .'...' t ~. !..t..!..'..t 6 7 8 '~ ! '....L' ' ! L.!...'...L,!...',..L!...'....'...L'....L'.,!..t...'...!..' 10 11 ! .!...'.'' !.L!! iii~!i~iiiiiiiii~i~iiii~iiiiiiiiiiiiii:~iiiiiiTiiiiTiT~ii~iiiiii~ii~7~T~iiii~iii~:~:!:2 i i•i i i•••ii i i i i i !i•!ii i i i i•!i!i•i i i i i i i i•••ii i i i i i•i i••iiTiTi i i•i•i •••••••••ii•i i • _ ~ ~ : ~ 1Z [ o ~[ ~, ~ i ~ ~ 2 0 I ~ 5 4 5 6 0 1 Z 3 4 5 BSTANCE 6 0 1 ~ 3 4 5 6 7 B 9 10 11 12 7 8 9 10 11 12 7 8 9 10 11 12 7 8 9 10 1] 1Z 7 8 9 1£1 11 1Z (km) D~STANCE (kin) 0 0 1 ~ 3 4 5 O I Z 5 4. 5 6 DISTANCE (km) o ' .!..!..t L.!...!..L..L.t..t .! ..L '.. L.~..t../...!...! ".. !_.'....t !_.L.t . ' ' i~#~ii ~i i i!:,!!:/#~iiiiiiiiiiii:,~,iiiiiii:iiiiiiiiii:,L~iliiiiiiii:i!!!!!i ii~iiiiiii:iliiiiiii i 2 1 3 4 .5 L . L . L . L L !...!...L!...L.!../..,L.'../...!...'...L..L 6 7 • 9 L.!.. f .t..t. ~..! ' 10 11 0 12 models: with homogeneous layers ( t o p ) , w i t h a l a t e r a l v e l o c i t y gradient in the second layer (second), with a s m a l l v e l o c i t y (third) and with a large velocity a n o m a l y ( b o t t o m ) . T h e numbers indicate the value of the velocities in k m / s . F i g . :2. G e o l o g i c a l anomaly Z 6 4 G. Boehm et al./Journal of Applied Geophysics 35 (1996) 1-13 in a later example) by observing the pattern of the estimated reflection points as a function of the offset. Small offsets are less sensitive to velocity errors and are, therefore, always closer to their true position. They thus suggest how to modify the interfaces to get a better approximation of the true solution: drawing a curve passing through the small offsets is always a good step ahead. The velocity field can be improved in two steps: by observing the pattern of the estimated reflection points as a function of the offset and by using the travel time inversion algorithm. If the estimated points corresponding to large offsets are closer to the surface than the points corresponding to small offsets, the local velocity must be increased, and vice versa. In general, we adopt the criterion of minimum dispersion: when the true velocity field is found, the images in depth provided by all offsets must coincide, and the dispersion of the corresponding reflection points must be minimum. OFFSET:~ ~ to g ¢'4 0 i'~ "~ g g ~ g g ~ '~r t¢) (O P"- CO O) ~ ~ Once the positions of the reflectors have been updated, they are fixed and only the travel times are inverted. This yields a further improvement of the velocity estimate. Among the many algorithms presented in the literature (e.g., ART, SIRT, conjugate gradient and others) we preferred the dual tomography (Carrion, 1991) due to its fast convergence and low dependence on the initial guess of the model. In the example shown later, we will start the inversion process with a homogeneous medium as the initial guess, i.e. without assuming any a priori knowledge of the investigated area. Under this condition, most algorithms do not perform well and some do not converge at all. A further advantage of dual tomography is its robustness with respect to the null space (Carrion, 1989). Our theory does not consider several factors affecting the velocity of the seismic waves, such as vertical velocity variations between reflectors, intrinsic anisotropy, viscoelasticity and thin layering, to ~ g g g ~ g g ~ g ~ ~ g t'N £'1 £'1 ~ ~C -- ~ COMMON c~ c-J SHOT section ~ ~ r~) ~ ~ ~ ~T TRACE [15] Fig. 3. Common shot gather No. 15, where the source is located at the left of the anticline, at 2890 m from the left limit of the model. G. Boehm et al. / Journal of Applied Geophysics 35 (1996) 1-13 mention but a few of them. Reflection tomography also provides us with a provelocity field (Al-Chalabi, 1994) or, to use a more common term, the parameters of a macro-model. Below it will, however, be seen that this provelocity field is much closer to the true velocity field than that supplied by the velocity spectra. smooth anticline and, at the bottom, a third horizontal interface at 1.8 km depth. We introduced the first horizontal layer to validate the algorithms, since the delays of the reflected arrivals are a simple hyperbolic function of the offsets. The anticline has dipping flanks, with a maximum dip of 30 degrees. Furthermore, the flat top of the anticline is at a depth close to the first horizon, to simulate a thin layer (300 m). The third reflector is used to analyze possible velocity anomalies caused by the overburden. We performed four simulations, corresponding to different velocity fields (Fig. 2). In the first case, the three layers are homogeneous: the values of the velocities are shown in the figure and the unit is k m / s . In the next model, the second layer has a linear horizontal velocity gradient ranging from 2 k m / s (left) to 3 k m / s (right). Thus, we get a wide range of velocity contrasts across the first two interfaces, obtaining a condition that is quite far from the hypotheses of the RMS velocity formula. The third and fourth models, on the other hand, have velocity 4. Features of the models The basic features of the models considered in this paper are displayed in Fig. 1. We generated 63 shots, moving the seismic spread from left to right. The first shot is located at 1 km from the left limit of the model, and the distance between shots is 135 m. For each source, indicated by a circle, we placed 48 receivers, regularly spaced 90 m apart with a minimum offset of 135 m, so that the maximum offset is 2295 kin. The structure of the reflecting interfaces is the same in all the models described below: a horizontal layer at a depth of 0.45 kin, followed by a OFFSET:~ i~ ,~t ~r ~ co r~ ~ ~ co o3 o 5 ~ (~ I~ ~ ~ ~ ~" ~') co ~ ~ r~ 0o (39 co £~ ~ c~ ~ ~ ~ ~ c~ e ~ v .~' . . . . ~ -- ~- ~ ~ ~ ~ r~ r~ r-~ TRACE COMMON SHOT section ( 3 0 ) ]Fig. 4. Common shot gather No. 30, where the source is located at the left flank of the anticline, at 4915 m from the left limit of the model. 6 G. Boehm et al./Journal of Applied Geophysics 35 (1996) 1-13 anomalies (linearly interpolated) below the anticline. The lateral extension of these anomalies is similar to that of the seismic spread. The synthetic data was modelled by the Fourier pseudospectral method, briefly outlined in Appendix A. This technique is efficient and sufficiently accurate for the purposes of this work: we verified it by selecting the reflected arrivals in simple models and comparing them with the known analytical solutions. Absorbing strips were set up at the boundaries of the model, in order to eliminate wave-field wrap-around (see e.g., Kosloff and Kosloff, 1986). The dashed lines in Fig. 1 show the inner limits of the absorbing area. The wavelet at the source is a zero-phase Ricker wavelet, with a dominant frequency of 30 Hz and a cut-off at 60 Hz. The seismograms recorded the time variation of the pressure at the surface with a sampling rate of 4 milliseconds. Figs. 3 and 4 show two common shot gathers obtained at the left side and near the top of the anticline, respectively, in the first model. The main difference between the gathers is the response of the anticline, i.e. the third arrival. In Fig. 3, the signals scattered back do not show a simple hyperbolic pattern, since they come mainly from the left flank and the corners of the anticline. The small bars indicate the selected travel times obtained by an automatic procedure that determines the arrival time with the maximum signal amplitude within windows provided by the user. Since the signal-to-noise ratio is particularly favourable in the data considered, this simple approach is sufficiently precise for our purposes. Figs. 5 and 6 represent two common offset gathers, with a distance between sources and receivers of 135 and 2025 metres, respectively. The zoom corresponds to the anticline in the central part of the model. Noticeable features are compression of the arrivals at large offsets, and some "noise" just below the third reflection in the central part: the latter is probably due to a complex diffraction and reverberation effect produced by the edges of the CDP: 0 ~ o 0 o ,~,0 ._ ~ 1 1 1 1 : COMMON OFFSET section 1 3 5 Fig. 5. Common offset gather (135 m), ~:~ ~: ~ i: - shot G. Boehm et al. / Journal of Applied Geophysics 35 (1996) 1-13 CDP: (I~ I0 I0 It') I0 I0 I0 co shot COMMON OFFSET section 2 0 2 5 Fig. 6. Common offset gmher (2025 m). The zoom corresponds to the top of the anticline. Stocking 2ooo i l l l l l velocity 25co . . . . . . . , . . . . . . . . . Stocking (m/s) 2 ioo~ 30oo 36oo k ......... i,, Ill velocity 25o0 ........ ' ......... Stocking velocity (m/s) 3o0o ' ......... 55o0 I, 2OO0 hi 'l (m/s) 3OO0 t'£CO . . . . . . . . . ' . . . . . . b,'' ......... 3500 ''' 0 25 -025 -0.25 0 25 0.5O -O5O -0.5O 0.50 0.75 -0~5 -075 075 1.00 100 1.25 125 11.~ 150 175 175 1.00 L25 i -1.25 oQ -1.50 150 175 - ,,,i . . . . . . . . . i . . . . . . . . . i .......... 2OOO 250O 30OO C.M,P. n. 80 ,,~r . . . . . . . . . i . . . . . . . . . i . . . . . . . . . 200) 2500 3000 5~() C.M.P.n. 90 ,,,, ......... 200~ ~ ......... 2500 , ......... 300) , 3500 C.M.P.n. 100 Fig. 7. Three different velocity spectra: from left to right they correspond to the left, the flank, and the centre of the anticline, respectively. Dark shading indicates high coherency values (close to 1), while the white background corresponds to low values (close to 0). 8 G. Boehm et al./Journal q[Applied Geophysics 35 (1996) I 13 oscillations corresponding to the flanks of the anticline, where we also obtained an anomalously high apparent velocity in the second layer and a low apparent velocity in the third layer. The errors are mainly due to the dip effect (Levin, 1971): multiplying the estimated velocity by the cosine of the angle of dip of the flanks yields a value that is closer to the true one. The velocity anomaly produces artifacts at the flanks of the anticline (see e.g., Loinger, 1983). Fig. 9 shows four stages of the tomographic inversion corresponding to Model 1. We start with a fiat model and a velocity distribution quite far from the true one, indicated by the dotted lines: this produces a noticeable dispersion of the estimated reflection points (top) and a large mismatch between the positions of the interfaces corresponding to the initial guess and those provided by the first iteration. The pattern of the reflection points suggests that the velocity of the first layer must be increased: this (second step) produces an exact reconstruction of the first layer, where the dispersion is minimum and the reflection interface coincides with the velocity boundary. Similarly, the same steps are then applied to the lower horizons. The final iteration (lower picture) gives a perfect reconstruction of the velocity field and the reflector positions. The interval velocity field of Model 2, obtained by conventional stacking velocity analysis, is shown in Fig. 10. The lateral gradient in the second layer, above the flank of the anticline, is fairly well-resolved. However, as in Model 1, the central part is anticline and the second interlace just above it. The same effect can be noticed at other offsets in the shot gather of Fig. 4. 5. Comparison of the two methods Fig. 7 shows three conventional velocity spectra corresponding to different CMP gathers (80, 90 and 100, respectively) from the left flat zone to the flank of the anticline. Well-defined peaks can be observed in the first and third velocity spectra: their interpretation is straightforward and leads to a correct estimate of the interval velocities. The second gather on the other hand displays two peaks at 1.2 s, which introduce an ambiguity. They are caused by the flat slopes of the flanks of the anticline, and the resolution depends on the cable extension. Naturally both are valid, but they cannot be represented simultaneously by a single stacking velocity, unless a dip move-out correction is applied. This phenomenon highlights the limited spatial resolution of a stacking velocity analysis: it cannot properly resolve two close but distinct events. Fig. 8 shows the interval velocity field corresponding to the first model, obtained by selecting both the conventional and the continuous velocity spectra along the synthetic profile. In the first layer and in the flat zone of the anticline, the velocity is perfectly recovered; there are, however, prominent C.M.P 50 0.~ ii ill lil 75 iIi lOO rlrl!l 125 rlllll 15o 175 200 iPplkrll ill llllpl 225 250 Jplll, ,r 275 509 rllllPlll 325 ,~1 ,,, 350 ,I,I, , 0,00 0.50 0,50 %~1.08 1,00 B m 50 75 100 125 150 175 200 225 250 275 500 325 350 Fig. 8. Interval velocity field obtained by the Dix formula, corresponding to Model I. G. Boehm et al./Journal of Applied Geophysics 35 (1996) 1-13 lateral velocity gradient was supplied by the initial iteration output (Fig. 11, upper picture): on the left the downward dipping alignments of the estimated points are evident, while on the right, there is an opposite trend. This suggests that the velocity was too high on the left and too low on the right. Fig. 12 shows the interval velocities estimated from the data of Model 3. There are pronounced badly resolved: the contrasting dips of the anticline surface produce velocity anomalies that are superimposed onto the lateral variations. Despite this fact the final tomographic image is very good, as can be appreciated in the lower picture of Fig. 11. We emphasise that the velocity gradient was estimated starting with a homogeneous initial model, without any a priori assumption. A hint to consider a possible (i<,,,) O~'ANC~ 1 t t ' ' I .... 0 j. . . . . . . . . . . . . . 0 2 i .... ~-. . 3 t .... . . . . i 5 p .... i .... . . . . . . . 9 . 6 i .... . . . . 7 r .... . . . . 8 t .... . . . . 9 r .... . . . . 10 t .... . . 11 t .... 2 3 4 5 1 2 3 4 $ 6 0 . 777-/77;777;77;77777;777UT,-77Z77UT.-L-777-;I' 1 12 7 fl 9 10 11 12 7 8 9 10 11 12 DISTANCE[ ( k i n ) 0 °iLL.'., . 6 _ , i _ , . _,' . _~~_~;,;". __;i '4_ _.._..__ , ' '_____i_''_;'_t° .....~_t~_":_~ ~t-- ~ ......................... I~' ,~/,'~.--. 0 1 2 3 4 5 0 1 2 3 4 5 S 7 8 ! 10 8 9 10 11 12 11 12 oL~rANCE 0<") 6 ...................... 2'~ 0 7 ................... ~ " .,_;_._,_;_,_;_to ~ 2 71-77-;71-77-;7~-771-7~ ~i-7~-;7~-77-;7~-7~-,'7~-;71-7 ;7777-7777;777; 1 l 3 4 6 5 7 IS i 10 11 i Ilzzt 10 Ill'l 12 DISTAI~E (kin) 0 ~/11 rl 1 2 3 4 5 6 [ l 1 1 1 1 1 L I l l l l t l l l l l l l l l l l l l Z r ...... 2 7 III 8 I I ' l l l l l km/~ 11 12 I t l l l / O ............ 2 0 1 2 3 4. 5 IS 7 8 9 10 11 12 Fig. 9. Four stages of the travel time inversion procedure corresponding to Model 1 (see text for details). G. Boehm et a l . / Journal q['Applied Geophysics 35 (1996) 1-13 10 C.M.P. 75 50 100 125 150 175 200 22b 250 275 300 350 325 0.00 0.00 0.50 ¸ 0.50 %.~_ 1.00 F~ 1,00 1,50 1.50 545 7'6 100 125 150 17,5 200 225 2,50 275 300 325 ,35(3 F i g . 10. Interval v e l o c i t y field obtained by the D i x formula, c o r r e s p o n d i n g to M o d e l 2. velocity oscillations in the range between 3.1 and 3.5 k m / s due to the interference of signals originating from the flanks of the anticline and the velocity anomaly below the anticline. The tomographic image (Fig. 13) indicates that the algorithm has succeeded in delineating the boundaries of the velocity anomaly, but has overestimated the value of the velocity (3.3 k m / s versus a true value of 3.2 k m / s ) . This error is due to the model discretization, since, in the inver- sion process, we chose long vertical pixels whose upper and lower boundaries coincided with selected reflected events. Therefore, since it is not possible to detect vertical velocity variations between the selected horizons, the inversion algorithm produces a local average velocity. Model 4 is similar to Model 3; the difference is that the anomaly below the anticline is wider. Its length is now similar to that of the acquisition DISTANCE (km) 0 2 Ilillllll[lllllllllll 3 , s o Jrlll'l'lblllllllllllll ............................ 1 2 3 2km/@ 4 5 7 8 9 ~o [llllllllllllliio ~ ............................. 6 7 B 9 ~2 I 10 11 12 DISTANCE(kin) 0 1 ,lltlLl*, 2 3 I I I Illlt 4 II I I I I l l 5 6 II 0 kIm I I/ s b i l l 7 B II 2I I I I I l l l 9 ill 10 11 12 ~111 tll I IIll~ T j 2 km/s 3 km/s I E3 2 ,,,iJ, 1 2 3 4 5 6 7 B 9 10 F i g . 11. T w o stages o f the travel t i m e inversion, c o r r e s p o n d i n g to M o d e l 2. 11 12 G. Boehm et al./Journal of Applied Geophysics 35 (1996) l-13 11 C.M.P. 50 75 lOO 125 15o 175 200 225 25(3 275 300 ,325 350 50 75 100 125 150 175 203 225 250 275 300 325 350 o.oo g,5o .~ 1.00 F- 1.50 Fig. 12. Interval velocity field obtained by the Dix formula, corresponding to Model 3. 0 1 ~E 2 3 DISTANCE (kin) 5 6 7 4 8 9 10 11 12 2.5 k m / s [l_ [ If' 1 2 ' ' l ' f ~ l l l r ' l l ' r I I l ' I ' ' l ' 3 4- 5 6 '~rl 7 ' ' ' 8 I l l 9 10 11 12 Fig. 13. Result of the travel time inversion, corresponding to Model 3, C.M.P. 50 75 1gO 125 150 175 L::K]O 225 250 275 303 ,325 550 (3.03 0.00 0...50- 0.50 %c~ 1 . 0 0 1,00 I-- 1.50 ¸ 1.50 50 75 100 125 150 175 200 225 250 275 300 325 350 Fig. 14. Interval velocity field obtained by the Dix formula, corresponding to Model 4. 12 G. Boehm et al. /Journal of Applied Geophysics 35 (1996) 1-13 DISTANCE (krn) 1 2 I=l~lllitlllt 3 IIIlllllllt 4 5 4 5 6 7 8 I I I I t l l l l l l l l l l l 9 10 11 12 I I I t l l l t l l t l ' l l l ' l l 0 T t~ "~ 2 1 2 a 8 7 8 9 10 11 12 Fig. 15. Result of the travel time inversion, corresponding to Model 4. spread. The interval velocity field (Fig. 14) contains more oscillations than the previous one (Fig. 12). The tomographic inversion (Fig. 15) achieves a much better reconstruction of the model then a conventional velocity analysis (Fig. 14), in particular, the spurious oscillations at the flanks and at the centre of the anticline have disappeared. 6. Conclusions Stacking velocity analysis provides a reliable image of the earth at depth only if the lateral velocity variations are gradual and if the dip of the reflectors is small. In all other cases, i.e. for most geological features of practical interest, reflection tomography is the proper tool to use. The main drawback of the travel time inversion is the selecting of reflected arrivals. This operation requires significant additional effort from the seismic analyst, especially if the signal-to-noise ratio is low and the angular coverage of the rays is high. The related increase in processing cost is not too high compared with other algorithms that are now standard practice, such as pre-stack migration, and negligible with respect to that of a dry exploration well. pressure, c(x,z) is the wave velocity, s(x,a,t) is the source and t is the time variable. The spatial derivatives are computed by using the Fourier pseudospectral method, where the transformed pressure P is computed by the Fast Fourier Transform (FFT). In particular, the algorithm used here is based on a vectorized version of the mixedradix FFT (Temperton, 1983). The steps of the calculation corresponding to the first partial derivative are as follows: Fvr FFT 1 3p p ~ P--+ ikP --+ - 3y (2) where v represents either x or = and k is the wavenumber. The method is very accurate up to the Nyquist wavenumber, which corresponds to a spatial wavelength of two grid points. This means that if our source is band-limited, the algorithm is free of numerical dispersion provided that the grid spacing is chosen such that D _< Cmin/(2fna×), where .~1~ is the cut-off frequency of the source and Cmin is the minimum phase velocity in the mesh. The time integration of Eq. (1) is carried out with the following second-order (staggered) differentiation scheme: q,,+,/2 = q , - , / 2 + d t ( D p " - s") (3) and Appendix A. The modelling technique p , + 1 = p,, + dtp,,+ 1/2 The pressure wave equation for a constant density medium is: where q = a p / O t and the time variable becomes t = n dt, with n a natural number. The differential operator D is given by: ~2p ~2p OX 2 -~- OZ 2 1 OZp C2 a t 2 + S (1) = where (x,z) is the position vector, p(x,z) is the 2[ 02 + C~2 t ] (4) (5) G. Boehm et al./ Journal of Applied Geophysics 35 (1996) 1-13 References AI-Chalabi, M., 1974. An analysis of stacking, RMS, average and interval velocities over a horizontally layered ground. Geoph. Prosp., 22: 458-475. Al-Chalabi, M., 1994. Seismic velocities - a critique. First Break, 12: 589-604. Bishop, T.N., Bube, K.P., Cutler, R.T., Langan, R.T., Love, P.L., Resnick, J.R., Shuey, R.T., Spindler, D.A. and Wyld, H.W., 1985. Tomographic determination of velocity and depth in laterally varying media. Geophysics, 50: 903-923. Carrion, P., 1989. Robust constrained tomography. Boll. Geof. Teor. Appl., XXXI: 167-200. Carrion, P., 1991, Dual tomography for imaging complex structures. Geophysics, 56: 1395-1376. Carrion, P,, Boehm, G., Marchetti, A., Pettenati, F. and Vesnaver, A., 1993a. Reconstruction of lateral gradients from reflection tomography. J. of Seism. Expl., 2: 55-67. Carrion, P., Marchetti, A., Boehm, G., Pettenati, F. and Vesnaver, A., 1993b. Tomographic processing of Antarctica's data. First Break, 11: 295-301. Dix, C.H., 1955. Seismic velocities from surface measurements. Geophysics, 20: 68-86. Hatton, L., Worthington, M.H. and Makin, J., 1986. Seismic data processing. Blackwell, Oxford. Harlan, W., 1989. Tomographic estimation of seismic velocities from reflected raypaths. Expanded Abstracts of the 59th SEG Annual Meeting, Dallas. 13 Hubral, P. and Krey, T., 1980. Interval velocities from seismic reflection time measurements. SEG, Tulsa. Kosloff, R. and Kosloff, D., 1986. Absorbing boundaries for wave propagation problems. J. Comp. Phys., 63: 363-376. Levin, F.K., 1971. Apparent velocity from dipping interface reflections, Geophysics, 36:510-516. Loinger, E., 1983. A linear model for velocity anomalies. Geoph. Prosp., 31:98-118. Lynn, W.S. and Claerbout, J. 1982. Velocity estimation in laterally varying media. Geophysics, 47: 884-897. Mayne, W.H., 1962. Common reflection point horizontal data stacking technique. Geophysics, 27: 927-938. Mayne, W.H., 1989. 50 years of geophysical ideas. SEG, Tulsa. Neidell, N.S. and Taner, M.T. 1971. Semblance and other coherency measures for multichannel data. Geophysics, 36: 482497. Sguazzero, P. and Vesnaver, A., 1987. A comparative analysis of algorithms for seismic velocity estimation. In: Bernabini et al. (Editors) Deconvolution and Inversion, pub. Blackwelk pp 267-286. Taner, M.T. and Koehler, F., 1969. Velocity spectra - Digital computer derivation and applications of velocity functions. Geophysics, 34: 859-881. Temperton, C., 1983. Fast mixed radix real Fourier transform. J, Comp. Phys., 52: 340-350. Yilmaz, O., 1987. Seismic data processing. SEG, Tulsa.