Geofisica Colombiana p P 5 -1 1 0121 - 2974N° 5 d l c i e mb r e de 2001 Bogota 0 C. ISS N Interval velocity determination by downward continuation of the traveltime function: Paraxial ray approximation LUIS MONTES V,DES Profesor asociado Departamento de Geociencias - Facultad de Ciencias - Universidad Nacional de Colombia. RESUMEN Present amos un rnetodo para estirnar las velocidades de intervalo, la profundidad y la geometrfa de los reflectores en modelos 3D que constan de un apilado de capas hornogeneas e isotr6picas con velocidades y densidades arbitrarias, separadas por interfaces suaves. El tiempo de transite de cualquier rayo reflejado en una interfaz particular y registrado en la vecindad de un rayo zero-offset se expresa mediante una funci6n referida a un sistema de coordenadas centrado en el rayo; tal funci6n se estima en la superficie superior del modelo. La funci6n tiempo de transite de reflexi6n asociada a cada superficie reflectora se deterrnina en la superficie superior en la vecindad del rayo central. La geometria de la superficie limitante superior de una capa particular y el tiernpo de transite estirnado sobre la misma permite calcular la velocidad de intervalo de la capa en cuesti6n y la geometria de la interfaz limitante inferior. Con la velocidad de intervale y la geometrfa de las interfaces limitantes, se estima la funci6n de tiempo de transite del siguiente reflector sobre la interfaz limitante inferior. En este paso se simula el posicionamiento de las fuentes y los detectores sobre la superficie anterior de la pr6xima capa subyacente. El proceso se repite recursivamente en las capas mas profundas hasta obtener la soluci6n completa sin conocimiento previo, excepto el obtenido en las capas superiores y la funci6n tiempos de transite de cada superficie reflectora. Se desarrollaron programas de computador que expresan el algoritmo del metodo y, posteriorrnente, se probaron con datos sinteticos, suministrando velocidades de intervalo y profundidades de los reflectores con errores considerados aceptables. PALABRAS CLAVE: TRAZAMIENTO DE RAYOS, TEORIA DEL RAVO, PROBLEMA INVERSO, CONTINUACION HACIA ABAJO, VElOCIDAD DE INTERVALO ABSTRACT We present a method to estimate interval velocities, reflector depths and geometries in 3D models consisting of a pile of isotropic and homogeneous layers of any velocities and densities separated by smooth interfaces. The travel time of a ray reflecting on a particular interface and registered in the vicinity of a zero-offset ray is expressed by a function referred to a ray-centered coordinated system, fnnction which is estimated at the uppermost surface of the model. The reflection travel time function associated to each reflecting surface is determined at the superior surface in the neighborhood of the reference ray. The geometry of the upper limiting surface of a particular layer and the travel time function estimated on this interface allow to calculate the interval velocity of the layer and the geometry of the bottom limiting interface. With the interval velocity and geometry of the two limiting interfaces of the layer, the travel time function of the following reflector is estimated at the bottom interface. This step simulates positioning the sources and detectors on the anterior surface of the next subjacent layer. The procedure is repeated recursively at deeper layers getting the complete solution without a priori knowledge but the upper deter- mined layers and the estimated travel time functions of each reflecting surface. Computer's programs expressing the algorithm of the method were developed and tested with synthetic data, providing the interval velocities and reflector's depths with errors considered acceptable. Key WORDS: RAY TRACING, RAY THEORY, INVERSE PROBLEM, DOWNWARD CONTINUATION, INTERVAL VELOCITY INTRODUCTION The interval velocity is a basic parameter in seismic data processing, specially in depth migration and interpretation of seismic sections. In the high-frequency range the propagation of seismic wave is well de- scribed by paraxial ray theory. In order to estimate the interval velocity using a travel time inversion method we recourse to the mathematical formalism developed by Bortfeld (1989). The Bortfeldss originally ideas were presented for earth models consisting of homogeneous lay- ers. Using dynamic ray tracing, Hubral et al. (1992) extended posteri- orly the validity of the formalism to layered media with lateral inho- mogeneity. The parameter describing a paraxial ray are its distance to GEOFislCA COLOMBIANA, 5, DICIEMBRE DE 2001 the central ray and the deviation of its slowness vector from that of the central ray. At any point of a paraxial ray, the parameters are linearly dependent on those at its initial point. The dependency is described by the propagator matrix, which is written as a product of many ray-seg- ment propagator matrices, each one associated to each layer. After Kahn (1987) the ray-segment propagator matrix through one layer is decom- posed in a product of matrix operators containing the parameters of the layer, and to determine the parameters of that layer is necessary to know the travel time function on its anterior surface. At the uppermost sur- face the travel time function of each reflecting interface can be esti- mated measuring the travel times of several rays in the neighborhood LUIS MONTES VIDES of the central ray. Based on the anterior ideas a method to estimate the interval velocities and reflector's depth was developed. Test on syn- thetic data showed that the algorithm representing the method provides results with errors considered acceptables. The solutions obtained are valid only in the second order approxi- mation of the travel times function and its range of applicability corre- sponds to the hyperbolic dynamic correction. TRAVELTIME FUNCTION According the formalism developed by Bortfeld (1989), the earth model consists of a pile of homogeneous layers of any velocities and densities. The layers are separated by arbitrarily curved but smooth interfaces. The uppermost interface, called anterior surface, is assumed to be plane and contains sources and detectors. The reflecting interface is called the pos- terior surface. One normal ray to the reflecting interface, called central ray, intersects the anterior surface at the origin of the xyZ coordinate system and the posterior surface at the origin of the x'y'z' coordinate system. In the vicinity of the central ray, any transmitted ray is described by its initial position X and slowness 13 vectors, and by its final posi- tion X' and slowness p' vectors. Representation is by two- component vectors. The two-component position and slowness vectors are obtained by projecting the three-component vectors into the tangent plane to the anterior and posterior surface at the intersection point with the central ray. The vectors X and 13 of any transmitted ray determine the vectors X' and p'. The first-order approximation of X and 13 sets a linear relationship: A, fl, ~ and D are 2x2 Jacobian matrices and Po is the slowness vector of the central ray. Similar equations apply for reflected events. According Bortfeld, in the second-order approximation, the travel time of a transmitted ray from any source position at the anterior surface X to any receiver position x' at the reflecting interface is given by: t(x,x') = to - Po ·x++x'· Dfl-Iy' +1..i· B-1 AX -i· B-1i' 2 - - - Two transmitted paraxial rays (one starting at X and X" other at) combine into a reflected one, with the travel time function t(x, X") given by the sum of the travel times-of the two transmitted rays: t(x, X") = to - 2Po + (x + x") + + (x + x") . D -I C 1.. (x + x") + 1.. (x" - x) . B-1 A 1.. (x" - x) - -2 2 - -2 where to is the two way travel time of the central ray. Expression (3) is called the parabolic approximation of the travel time along the paraxial ray, and is not the best approximation. It is known a long time ago that in simple layered media, seismic near- vertical reflections are better approximated by hyperbolic rather than parabolic travel times curves. Ifwe square (3) and retain only its terms up to the second-order in X and x', the hyperbolic approximation of the travel time for the paraxial ray is obtained (Ursin, 1982). ( ")2 _ I" 2t X,X -(to-2po'z(x+X)) + 2to(1(X + X")· o' c t (X + X") + t (X" - X)· B-1 A t (X" - X)/ (4) Equations (3) and (4) contains the same nine unknown parameters: one in to' two in Po' three in ll-IA and three in D -I ~ (because the matrices are symmetric), in consequence (3) and (4) can be deter- mined registering at least nine travel times for rays in the vicinity of a central ray at uppermost surface. For each reflector interface in the model, a travel time function like equation (4) will be estimated at the uppermost surface. INTERVAL VELOCITY Considering the travel time function determined at the anterior surface of a layer and known the geometry of that interface ( K N-I)' we will show how the interval velocity of the layer can be determined, uniquely. The total transfer matrix, which describes the ray propagation in the N" layer, can be decomposed after Kahn (1987) in a product of matrices: (I) [AN BN] [ I a~VNSNI] eN DN - V;;IKN 1 [~ o ][ L ~] (5) -I 1 -IL -COS({3N)V;; KN-I [COS({3N) ~]with L = 'L({3N) = 0 (6) (2) The capital letters in the matrix relation (6) are 2x2 matrices: e.g., 1 is the unity matrix, Q is the zero matrix and K N is the matrix of curvature of the N" interface. The N" interval velocity is VN ' and SN the path length of the central ray in the layer. The angle between the normal vector of the (N - I )th interface at the intersection point with the central ray and the slowness vector of the central ray is f3 N . The ma- trices on the right side of the transfer matrix are responsible for the following transfers of positions and slowness vectors, from right to left: (I) from the (N - I )th interface into the tangent plane, (2) from the tangent plane into the plane perpendicular to the central ray, (3) from the plane perpendicular to the central ray at the (N - I )th interface into a perpendicular plane to the central ray at the N'" interface, and finally from the tangent plane into the N" interface, i.e. into the reflector. !1t N is the one way travel time of the central ray through the N" layer. In this step, besides K N-I ' we assume to know the interval velocity of the superior layer VN-I and the angle of incidence of the paraxial ray aN _I at the anterior surface. Multiplying the matrices in expression (5), we obtain the following equations: (3) GEOFislCA COLOMBIANA, 5, DICIEMBRE DE 2001 INTERVAL VELOCITY DETERMINATION BY DOWNWARD CONTINUATION OF THE TRAVELTIME FUNCTION, I'ARAXIAL RAY APPROXIMATION AN ='L({3N)-COS({3N)SN'L-1({3N)KN-1 EN = VNSN'L-'({3N) c, =V;/KN'L({3N)-V;;'COS({3N) [SNKN + L]'L-' ({3 N )KN-' DN =[SNK N +!JI-I({3N) Determining the inverse of !l N from (8) and multiplying by (7) we obtain !1~'AN' that we recalled U: equation (12) is equivalent to the following three equations: 1 1 U22 = -- - -COS({3N )KN_, 22 VNSN VN ' At the anterior surface of the layer the Snell's law is used: COS\{3N)= 1-Sin2(aN_1){ VN }2 VN_1 where aN_I is the angle of incidence and {3N is the angle of trans- mission of the central ray, see figure I. By definition the interval velocity is: v - 2SN N - (17) /)"t N Solving simultaneously the equations (13), (14), (15), (16) and (17), see Appendix A, the interval velocity of the layer VN is determined, except in case of a spherical reflector as was pointed by Krey (1989).This step permits to know also the three others unknowns SN' /)"t N' fJ N ' and in con sequence to determine A Nand !l N . GEOMETRY OF POSTERIOR SURFACE To determine the geometry of the posterior surface is necessary to know SN' /)"t N' fJ N and VN' i.e., which were solved the anterior step. To calculate the inverse of D N Equation (10) is used: D~' =I({3N)[SNKN -rr' the product of (18) by (9) gives the matrix GEOFISICA COLOMBIANA, 5, DICIEMBRE DE 2001 (7) (8) T(Sn, Gn) ~N.J (9) anterior surface L-----i--~ (10) (12) posterior surface (13) (14) Figure J. Left: A seismic system and the travel lime function of rays reflected on the posterior surface, sources and detectors on the anterior surface. Right: The central ray cross the anterior surface with angles of incidence and transmission. (15) ( 16) (19) Reorganizing terms in (19) VN D-;" CN = 'L({3 N )[S N KN + Ir' {KN 'L({3 N) - COS({3N)[S N KN + I]'L-' ({3 N)KN-d (20) then 'L({3N )[S NKN + Ir' KN 'L({3N) = VNDN'CN +cos({3N)KN-1 (21 ) to find: [SNKN +Ir'KN = 'L-1{({3N)VNDN1C.N + COS({3N)KN-I }'L-' ({3N) (22) term simplified by recalling ( I 8) Q = I.-' ({3N) {V N D'N' c, + Cos({3 N)KN-' }I.-1 ({3N) to get: (23) LUIS MONTES VIDES ~---I------7 T3 ~ __ -+-__ -;? T2 ~----t-----? T1 r 3 V 3 Figure 2. The model is analyzed like three seismic systems, each one with its own central ray and each reflector like its posterior surface. There are travel time functions associated to each system and estimated on the anterior surface (flat). Finally the matrix of curvature of the posterior surface in the vicinity of the central ray K N is obtained: K = Q{I - S Q}-l_N _ _ N_ Now, we can know numerically the transfer matrix which describes the ray propagation through the N" layer, substituting the terms VN' SN' /).t N' f3N' K N in equation (5). TRAVELTIME FUNCTION IN DEPTH It will be necessary to introduce different central rays in the model, a central ray for each "piece" of seismic system. The first niece of seis- mic system is defined between the uppermost surface (zero interface) and the first reflecting surface, the second between the first one and the second one, and so on. Each piece is considered a seismic system with a central ray perpendicular to the corresponding reflecting surface at the midpoint position (see figure 2). . The determination of the travel time function in depth is interpreted like if sources and detectors were moved from the uppermost surface, following the raypaths, to the anterior surface of the deepest layer of the seismic system (see figure 3). We will show how to determine the travel time function at the anterior surface of the deepest layer from the knowledge of the travel time functions at the upper interfaces of the seismic system. The propagator matrix describes the ray propagation of the paraxial rays through the system, and can be decomposed in a product of many- ray segment propagator matrices. In a system with N layers the propa- gator matrix is written as the product of the N individual propagator matrices of each "piece": B ]k~N-l[A B]_N II -] -] D . C D_N ]~l -] -] (25) If we have solved the system composed by the upper N - Ilayers, i.e., the transfer matrix for each superimposed layer is known, then (25) is rewritten: (26) where the matrix in the middle of (26) propagates the ray in the Nth layer (assumed unknown), and the left most matrix in (26) corresponds to the product of the N - Ipropagator matrices of the upper N - I solved layers, i.e., it has been assumed that the transfer matrix for each superimposed layer is numerically known. The transfer matrix posses the property of simplecticity (Borfeld & Kempert, 1990), this property means that the inverse of the propagator matrix of the system composed of N - Iupper layers is: (27) the super index means the matrix transposed, and by the same argu- ment the inverse of the propagator matrix of the N" layer is: [~ !l]-I [D' - B ]__ N _N !2 -~N A~ (28) (24) using (27) and (28) we transform (26) in: [ [2~ -C'_N [ D~ -C'_N !In.][ [2' - !l'] D -C' A'-n· _ _ (29) su~dIPoKsurfa \~ ~0 layers : Traveltime An . ~functionin tenor de th surfa~cs*\- -/G* P Unknown layer poste~ V ~ <»: surfa~~-' ~ Figure 3. If the sources and detectors were moved from the superior surface to the ante- rior surface following the ray paths, then the travel time function is determined at the top of the deeper system. GEOFf SICA COLOMBIANA, 5, DICIEMBRE DE 2001 INTERVAL VELOCITY DETERMINATION BY DOWNWARD CONTINUATION OF THE TRAVELTIME FUNCTION, PARAXIAL RAY APPROXIMATION factorizing D' at left of (31) D~CN =[DnkC -CnkD'][(Ank -BnkC D-1)D'r1 (32) The product of symetric matrices is conmutative, then D~CN =[DnkC -CnkDl]Dl[(Ank -BnkC D-1)J1r1 (C33) Dt)CN =[Dm(C' D-')-Cnd[Ank -Bw(C D-t)r1 (34) by symmetry D-1 ~ = C' D-' then D~' ~N = [Dnk (D-I~) - ~nk HArlk - ~nk (D-1 ~)]-I (35) Following a similar process we demonstrate that ~~IAN = U2nk [~-I A] - ~nk HAnk - ~nk [~-I An (36) This step makes possible to know the matrices u; AN and !2~1~N in depth. The two-way time ToN and the slowness vector jj oN complete the knowledge of the travel time function in depth, these last two parameter are calculated by ray tracing through the model, as it is shown in the next procedure. RAY TRACING In this step, the interfaces limiting the N upper layers and the their interval velocities are known. The paths of rays traveling across those layers can be calculated (see figure 4). In the vicinity of the central ray, the Nth interface is approximated by the second order polynom Z = + (x, y) . K N (x, y) , where K N = (k11 ' k12, kn) is the matrix ~f curvature. Central 1'1 ray r3 Paraxial rays v. Figure 4. Ray tracing of the c;,nl~~l and the paraxial rays through the last solved layer. The new coordinate system (I. ,.I ) has rotated an angle e GEOFislCA COLOMBIANA, 5, DICIEMBRE DE 2001 According figure 5, the normal ray to the N + I interface is refracted at the N surface in the point: (x, y, z) = -rl + £...r + r2 (37) where SI' SZ are the respective path lengths of the paraxial and central ray in the Nth layer, ~ = (0,0, S I) and £...r = (£...rl ' /sr, ,£...r3) are known. Therefore (37) is equivalent to equation (38). dr\ r2 laW N Normal ~~ector (x, y, z) II II Figure 5. The ray '2 is normal to the next interface in depth, been refracted in (x, y, z). a, f3: incidence and transmission angles. Substituing X, Y in the third expression of (38) we obtain: S;[sinz f3(~ kll cos' 8 +~ k22 sin' 8 + k'2 sin8cos8)] + S2[sinf3{k11 cos8fuj + k2Z sin8&z + 'sz(sin8fuj + COS8&2)} - cosf3] + [S, -11r" + ~(kl,fujz + k22&Z2 + 's2fuj&2)] = a (39) which has the form AS; + BS Z + C = O. With S2 known, we calculated: the position (x, y, Z), the normal vector in this point - v· by N = 1'11:::1 and the incidence angle in the Nth interface by cos aN = - f,"1 . Using the Snell's law, the refracted angle f3N+I of the paraxial ray is determined. Now is possible to estimate the time used by the paraxial ray to travel through the upper known layers and the remaining time to reach the posterior interface of the next unknown layer, i.e, we estimate the travel time function at the anterior surface of the unknown layer, to be solved. APPLICATION TO SYNTHETIC DATA Synthetic data were generated in several models using the ray tracing pro- gram Anis ray_3D, developed by Costa et al, at the UFPa in Brazil (1993), but only two of them are shown here. Two programs were developed, the first one to estimate the travel time function at the uppermost interface of the seismic system, and the second one to determine the interval velocities by the estimation of the travel time in depth as was explained in before. LUIS MONTES VIOES o -05 -I ~-15 r:r ·2 -25 -3 -5 XIKml Vjm/s ] Vlm/s ] Erro[m/s] Erro[m/s] 0 0 0 0 200 200 200 200 400 400 400 400 GOO GOO GOO GOO 800 800 L- 800 LJ 800 K 1000 000 E 000 K 000O£ O£ N N N N '1200 200 200 200 1400 400 ) 400 400 11 GOO GOO 600 600 1800 800 800 800 2000 000 000 000 0 5000 0 5000 0 50 0 50 Figure 6. On top the model used to test the algorithm. On bottom interval velocities and erros obtained in two points at the superior surface. Besides, the second program also calculates the matrices of curvature of the interfaces an the points where the rays intercept them. Figure 6 shows on top the first model used to generate the travel time and in bottom the interval velocity model obtained by the pro- gram, in two different positions fo the model. The two leftmost graph- ics in figure 6 shows the errors in estimating the interval velocities, being fewer than 50 mls. A second model was used to show the program estimating the points where the rays reach the interfaces; the model is shown in top, figure 7, and in bottom the final position reached by the rays calculated by the program. In the third interface the program srnoothes it, satisfying the restriction impossed on the surfaces of the models. The others points show a well fitted with the interfaces. CONCLUSIONS The second order approximation of the travel time and the paraxial ray theory make possible to develop an inversion method to estimate the interval velocities in models with no structural complexity or strong dipping interfaces. The method is a tool to estimate interval velocity and reflector depths, based on post stacked seismic data. Due to the paraxial ray theory difficulties could occur when the positions of sources and detectors exceed the validity of the approximation, being valid only for small values of distance to the zero-offset ray position. Its range of applicability correspond to the hyperbolic dynamic correction. Because 10 :~--- -----------0----5 5 ·5 Figure 7. On top the second model and on bottom the interfaces and the points where the algorithm estimates final positions for normal rays. the input data are travel time measurements on seismic traces, a good quality on seismic data is impossed. Nevertheless the method provides the inverse solution without a priori knowledge of the model, just using travel time measurements on the uppermost surface. BIBLIOGRAPHY BORTFELD,R. 1989: Geometrical ray theory: Rays and travel times in seismic sys- tems. Geophysics, 54, 342-349. BORTFELO,R., and KEMPER, M. 1991: Geometrical ray theory: Line foci and point foci in the anterior surface of seismic systems (second order approximation of the travel times). Geophysics, 56, 806-811. COSTA, 1., SCHOENBER, M. and URBAN, 1. 1993: 3D Ray tracing through no-flat homogeneous anisotropic layers. 3rd International Congress of the Brazilian Geophysical Society. Expanded Abstracts, vi, 47-52. HUBRAL, P, SCHLEICHER,J. and TYGEL, M. Three-dimensional paraxial ray prop- erties - I Basic relations: 1. Seis. Expl., I, 265-279. KAHN, S. U. 1987: Die Aufdeckung der Mehrdeutigkeiten in der Reflexionsseimik: Ph.D. tesis, Tech. Univ. Clausthal. KREY, T. 1989: The non uniqueness of the determination of interval velocities from moveout velocities. Geophysics, 54, 1209-1211. URSrN, B. 1982: Quadratic wavefront and travel time approximations in inhomoge- neous layered media with curved interfaces. Geophysics, 47, 1012-1021. ACKNOWLEDGEMENTS This paper is the result obtained in the Montes's Ph. D. thesis at the Federal University of Para, Brazil. The project was supported by a fel- lowship from Coordenacao de Aperfeicoamento de Pessoal de Nivel Superior (CAPES). APPENDIX DETERMINING THE INTERVAL VELOCITY To determine the interval velocity in the N" layer is necessary to solve the set of simultaneous equations, represented by the matrix equa- tion B-1 A: GEOFlslCA COlOMBIANA, 5, QICIEMBRE DE 2001 INTERVAL VELOCITY DETERMINATION BY DOWNWARD CONTINUATION OF THE TRAVELTIME FUNCTION, PARAXIAL RAY APPROXIMATION (A I) (A2) 1 UI2 = ---COS ([3N)KN-I12 VNSN ' the Snell's law at the anterior surface of the Nth layer (A3) COS2([3N)=I-Sin2(aN_1){ VN }2 VN_1 with angles aN-l of incidence and f3 N oftransmission, and the inter- val velocity (M) 2SN VN =& (AS) N In order to solve (AI), (A2), (A3), (A4) and (AS) we will consider five different cases: First: when k N-I,ll = k N-I,12 = k N-1,22 = 0 , i.e, the anterior sur- face is a flat one, or when kN-1,12 = kN_1,22 = 0 and kN_1,11 '" 0, sub- stituting (AS) in the reduced equation (A3) (M) Second: when k N-1.I2 * 0 we use equations (A2) and (A4) to get [ k ]2[ 22_ N-1,l2 ,2 VN VN - u I - Sin (0. N - I ) V 2 - 12 ~I and finally (A 7) VN=[kUI2 ]2+(Sin~(~N-I)J (A8) N-1.12 N-I Third: when k N-l,12 = k N-I,ll = 0 and k N-I,22 '" 0 we substitute (AS) in (AI) _ 1 [ 1 ' 2 ) V; ]ull - -- 1- - Sin (aN-1 -0- VNSN VN'_I to get (A9) V - I1tNu11 + sin 2 (aN_1) N- 0 (AID) 2 VN'_I Fourth: when kN-1,22 * kN_1,1I * 0 and kN-1,12 = 0 we re- write equation (AI) COS 2 ([3N) VNSNKN_1,11 (A It) KN_I,II and equation (A3) GEOFislCA COLOMBIANA, 5, DICIEMBAE DE 2001 ~ 1 COS(f3N) KN-1,22 VNSNKN-1,22 VN (AI2) Subtracting (A 12) from (A 1I) and factorizing the factor VN SN we obtain V S ( - kN_1,ll J+ ~N-l,11 = N NUll Ul2 kN-1,12 kN-1,12 1 ' 2( ) VN-Sin aN-1 -- VN_1 (AI3) including (AS) in (Al3) and reorganizing tenns,the interval velocity is V -N- (A 14) r: k-I N-I - N-I,n Fifth: when kN-1,22 = kN_1,11 * 0 and kN-1,12 = O,i,e"thean- terior surface is a spherical one, From (A3) VNU22 KN-1,22 (AIS) and using (A4) in (A IS) and squaring the resulting expression we have 1 '2( )[VNJ2-( 1 VNUnJ2-Sin a -- - ---- N-l VN_1 kN-1,22S N kN-1,22 and substituting (AS) in (A 16) (A 16) (A 17) The expression (A 17) is reorganized in terms of VN (A l8) The expression (AI8) has the form a(V[~) +b(V; )+C = 0 and its solution is given by V~ = (2a)-1 (-b ±.J b2 - 4ac), in conse- quence the square of the interval velocity can be known but there is an ambiguity to determine the interval, i.e., the unique solution does not exist for this particular situation. The result shown here was first re- ported by Kahn (1987) and later by Krey (1989). 11