Imagine a boat releasing impulsive sounds waves at equal intervals of time. The pressure waves propagate down to the seabed and then deeper into the earth: they are reflected by both the seabed and from underground layer boundaries. Seismic data are produced by towed receivers which record the amplitudes and times of the reflected signals. Our goal is to build a fast and robust algorithm for finding the sound speeds ("seismic velocities") inside the earth from these data. Careful measurements of the seismic velocity are part of obtaining accurate seismic images showing layers and cracks inside the earth, and are used to determine geologic features, such as the location of oil trapped at the sides of a salt dome, shown here by the red color. 
A location of a point A inside the earth can be described in either of two coordinate systems:
Depth coordinates, on the left above, give a location of a point A inside the earth in terms of its position on the top and depth below. In contrast, time coordinates, shown on the right, think of each point A under the earth's surface as corresponding to a pair (x0, t0): if you imagine a wave starting at point A, x0 is where it first hits the earth's surface, and t0 is the time when it does so. While it may be counterintuitive, seismic data are most naturally first recorded in time coordinates rather than depth coordinates.
The two main approaches to seismic imaging produce models for the underground velocities: the first is a process called time migration, which takes seismic data in time coordinates, and produces images and timemigration velocities, which are an averaged velocity of a particular type. The second, depth migration, takes seismic data in depth coordinates and produces seismic images in depth coordinates.
Time migration

Depth migration



Adequate for

mild lateral velocity
variation 
arbitrary velocity
variation 
Implementation
requires 
seismic data

seismic data
+ velocity model 
Produces images in

time coordinates

depth coordinates

For the horizonally constant case the relationship between the seismic velocities and the time migration velocities vm(x, t) was derived by C.H. Dix in 1955:
In order to use depth migration, we would like to understand how time migration velocities relate to true seismic velocities in the general case of horizontally variable velocity. The Dix velocities vDix(x, t)computed from the time migration velocities vm(x, t) by the formula above serve as a more convenient input. Ultimately, we are faced with an inverse problem:
Nonetheless we have attempted a regularized reconstruction for a short enough interval of time. (The seismic data are typically acquired only up to 2 seconds and up to 5 km in depth.) We have found two ways to solve this PDE.
We have shown that these schemes work because of the following reasons
Once we obtain Q we obtain the seismic velocity v in the time coordinates. The next step is to convert the seismic velocity to depth coordinates. To do this we have developed an efficient, Dijkstralike timetodepth conversion algorithm. It solves the Eikonal equation with an unknown righthand side: it does this by systematically building the velocity field by coupling the Eikonal equation with an orthogonality relation. Its motivation and a building block was the fast marching method. This algorithm is considerably faster and more robust than existing techniques. However, because we are required to solve two coupled equations simultaneously in order to build the seismic velocity in depth coordinates, a very subtle issue arises as to whether one can maintain causality in systematically building the solution. After an exhausting six months, the answer turns out to be "yes".
We generalize the PDE and our finite difference numerical scheme to 3D, and test our numerical techniques on a collection of synthetic examples, demonstrating that we are able to restore the seismic velocity quite accurately. Results are compared with the standard Dix estimate, and demonstrate that the Dix estimate might differ qualitatively from the original velocity while our correction gives a significant and qualitative improvement to the Dix estimate. Our test on the smoothed Marmousi data confirms the effectiveness of the proposed approach.
Symmetric Gaussian anomaly.
Top: the exact velocity. Middle: the input: the Dix Velocity. Bottom: the reconstructed velocity and the rays. 
Asymmetric Gaussian anomaly. Top: the exact velocity. Middle: the input: the Dix velocity. Bottom: the reconstructed velocity and the rays. 
3D Gaussian anomaly. Top row: xzplane; Middle row: yzplane; Bottom row: xyplane at 2.55 km in depth. 1st column: the reconstructed velocity; 2nd column: the exact velocity; 3rd column: the Dix velocity. 
3D archshaped Gaussian anomaly. Top row: xzplane; Middle row: yzplane; Bottom row: xyplane at 2 km in depth. 1st column: the reconstructed velocity; 2nd column: the exact velocity; 3rd column: the Dix velocity. 