This invention is relevant to seismic data processing in the field of geophysical exploration for petroleum and minerals. The general seismic prospecting method involves transmission of elastic, or "seismic," waves into the earth and reception of reflected and/or refracted waves at the earth's surface (or, occasionally, in a wellbore) via geophones, hydrophones, or other similar devices (hereinafter referred to collectively as "geophones"). The elastic waves may be generated by various types of sources, dynamite and hydraulic vibrators being particularly common. As these waves propagate downward through the earth, portions of their energy are sent back to the earth's surface by the acts of reflection and refraction which occur whenever abrupt changes in impedance are encountered. Since these impedance changes often coincide with sedimentary layer boundaries it is possible to image the layers by appropriate processing of the signals returned to geophone groups.
Many methods of seismic data processing in use today require calculations using traveltimes to determine information regarding subsurface geology. "Traveltime" means generally the amount of time a seismic signal takes to travel from a seismic source to a subsurface reflection point to a seismic receiver. For example, the concept of migration is well known in the art. Simplistically, raw seismic data as recorded are not readily interpretable. While they show existence of formation interfaces, raw data do not accurately inform the interpreter as to the location of these interfaces. Migration, also called imaging, is the repositioning of seismic data so that a more accurate picture of subsurface reflectors is given. In order to perform the migration calculations, the seismic velocities (hereinafter referred to as "velocities") of the subsurface at a multitude of points must first be determined or approximated. These velocities are often estimated using traveltime information. Current methods of computing the traveltimes necessary to perform three-dimensional depth migration and associated velocity analyses are inefficient and/or potentially error-prone when applied to the complicated velocity models typically encountered.
There are currently two general methods of determining the grid of traveltimes needed to migrate data: the well-known ray tracing methods; and the more recently proposed methods which seek a direct solution of the eikonal equation. As is known to one skilled in the art, the eikonal equation is a form of the wave equation for harmonic waves, valid only where the variation of properties is small within a wavelength, otherwise termed "the high-frequency condition." In rectangular coordinates, the eikonal equation is as follows: EQU (.differential.t/.differential.x).sup.2 +(.differential.t/.differential.y).sup.2 +(.differential.t/.differential.z).sup.2 =S.sup.2 (x,y,z) (1)
where the coordinate axes are x, y, and z; t is the traveltime; and S is the slowness, which is the inverse of velocity. This equation relates the gradient of traveltime to the velocity structure.
Seismic ray tracing methods are applied to determine traveltimes in most applications used commercially today. Ray tracing equations are linear, ordinary differential equations derived by applying the method of characteristics to the eikonal equation, a technique known to those skilled in the art. Ray tracing allows the determination of arrival times throughout the subsurface, by following raypaths from a source location, which raypaths obey Snell's law throughout the subsurface volume for which the velocity distribution is known. Traveltimes along the rays are then interpolated onto a three dimensional grid of the subsurface.
The ray equations may be solved with shooting methods ("shooting") or with bending methods ("bending"), as are well known to those skilled in the art. Shooting formulates the ray tracing equations into an initial-value problem, where all ray direction and position components are defined at the source location at time zero. Then the equations are recursively solved to trace the rays throughout the medium. Bending is based on Fermat's principle, which states that the seismic raypath between two points is that for which the first-order variation of traveltime with respect to all neighboring paths is zero, and attempts to locate a raypath between two points by determining a stationery traveltime path between them. It formulates the ray tracing equations into a two-point boundary value problem. Shooting is generally more efficient computationally than bending; however, both approaches present difficulties and potential inaccuracies when used to compute the gridded traveltime fields required by three-dimensional depth migration. Three-dimensional depth migration typically requires robust grids of traveltime for high quality images. As used herein, the term "robust" means a process which reliably generates accurate grids of traveltimes regardless of velocity model complexity.
Complications arise in shooting and bending calculations because each ray is computed independently of all others, and because small changes in the angle of incidence may lead to large changes in ray direction. This complication is a manifestation of the nonlinearity of the problem of ray tracing. Thus, with shooting methods, it is difficult to obtain the uniform ray coverage throughout the model that will admit traveltime interpolation onto a grid. Also, rays may cross, indicating the natural multivaluedness of wave propagation. However, when interpolating traveltimes onto a grid, one may not interpolate between different branches of a multivalued wavefront, hence a "smart" interpolation, that is, an interpolation geared to adjust its parameters based upon conditions at various locations within the grid, is necessary for a multivalued wavefront. Unfortunately, a "smart" interpolation is a difficult task in three dimensions.
Ray shooting itself is quite efficient, but a robust application using it to compute grids of traveltimes may not be efficient. While the bending method will determine a raypath to any point in the grid, it is subject to local, instead of global, minimization, resulting in traveltime errors. Bending is also inefficient for determining three dimensional grids of traveltimes.
An alternative to ray tracing is to solve the three dimensional eikonal equation directly, using numerical techniques. This method was proposed in Vidale, J. E., "Finite-difference Calculation of Traveltimes in Three Dimensions," Geophysics, vol. 55, pp. 521-526 (1990), Vidale succeeded in calculating a grid of traveltimes with a finite-differencing scheme based on the rectangular version of the eikonal equation as set forth above in Equation (1), which scheme makes a plane-wave approximation, a simplifying assumption which introduces error into the calculation. An alternative to Vidale's method was proposed in Podvin, P. and LeComte, I., "Finite-difference Computation of Traveltimes in Very Contrasted Velocity Models: a Massively Parallel Approach and its Associated Tools," Geophys. J. Int., vol. 105, pp. 271-284 (1991), which also makes a plane-wave approximation, which approximation is not used in the method of the present invention. This alternative constitutes a finite-difference approximation of Huygens principle. Huygens principle holds that every point on the boundary of a "cell," wherein a plurality of these cells define a subsurface grid, acts as a secondary source emitting an impulse at the moment it is reached by the first wave arrival.
These algorithms do not vectorize or parallelize at the loop levels. This means that these algorithms are not set up to contain "independent" nested iterations of calculations or "loops". A vectorizable algorithm contains single unnested do loops whose iterations may be computed in any order (thus they are independent). A vectorizable and loop parallelizable algorithm contains doubly nested do loops, whose iterations are independent of one another. The inner loop vectorizes and the outer loop parallelizes. Vector loops may be executed very efficiently on any one processor of a vector computer, while parallel loops may be executed efficiently by multiple processors on a parallel computer. A CRAY.TM. Y-MP is an example of such a machine.
The two-dimensional graph theory approach of Moser, T. J., "Shortest Path Calculation of Seismic Rays," Geophysics, vol. 56, pp. 59-67 (1991), extends to three dimensions; however, this approach, also a Huygens principle method, is also inefficient.
Van Trier, J. and Symes, W. W., "Upwind Finite-Difference Calculation of Traveltimes," Geophysics, vol. 55, pp. 521-526 (1991), drew upon fluid dynamics technology and presented a two dimensional upwind (meaning only traveltimes at points in the direction of the negative traveltime gradient, with respect to the current calculation point, are used in the determination of the traveltimes) finite-difference traveltime calculation technique that fully vectorizes, and Popovici , A. M., "Finite-difference Traveltime Maps," Stanford Exploration Project Report vol. 70, pp. 245-256 (1991), presented the extension of the Van Trier and Symes method to three dimensions, while conceding that the algorithm was unstable. The three dimensional algorithm radially extrapolates the three components of the slowness vector (where slowness is the inverse of velocity), i.e., .differential.t/.differential.r, (1/r).differential.t/.differential..theta., and (1/[r sin .theta.]).differential.t/.differential..phi., in a spherical coordinate system, where
t=traveltime; PA0 r=the radius from the origin point of the spherical coordinate system to the point of interest; PA0 .theta.=the angle that the z-axis makes with the radial line connecting to the point of interest; and PA0 .phi.=the angle from the x axis to a projection of the radial line onto the x-y plane.
While the efficiency of Popovici's algorithm makes it attractive, it is, as Popovici conceded, unstable. It is particularly unstable and inaccurate when applied to the highly variable velocity fields encountered in three-dimensional depth migration applications. Popovici, "Stability of finite-difference traveltime algorithms," Stanford Exploration Project Report , vol. 72, pp. 135-138 (1991), presented a geometrical derivation of the Courant-Friedrichs-Lewy (CFL) stability condition for the two dimensional method. However, no one has been able to improve the Popovici three-dimensional solution so as to remove its inherent instabilities.