Tomographic velocity inversion

Part 3, Velocity Analysis & NMO

Learning objectives

  • Set up a tomographic inverse problem G Δs=Δt\mathbf{G}\,\Delta\mathbf{s} = \Delta\mathbf{t}, where G\mathbf{G} holds the ray lengths in each cell, Δs\Delta\mathbf{s} is the slowness update and Δt\Delta\mathbf{t} the travel-time residuals
  • Explain why tomography solves for the full model simultaneously instead of iterating per-CMP picks
  • Describe the straight-ray and bent-ray (iterative) variants and when each is appropriate
  • Identify the resolution limits and common failure modes (uneven illumination, regularization trade-offs)

Conventional velocity analysis (Sections 3.1-3.5) picks V(t_0)V(t\_0) at each CMP independently. That works when lateral variation is smooth. When it is not (foothills, reefs, karst, buried river valleys), neighboring picks can be inconsistent, and the per-CMP approach fails to reconstruct the full velocity field. The answer is tomography: simultaneously invert the travel-time residuals everywhere in the survey for the velocity model everywhere.

1. The inverse problem

Discretize the subsurface into cells (the model grid) with per-cell slowness s=1/Vs = 1/V. For each measured ray, its travel-time residual (observed minus predicted) is

Deltat_textray=sum_textcellsell_textray,cellcdotDeltas_textcell\\Delta t\_{\\text{ray}} = \\sum\_{\\text{cells}} \\ell\_{\\text{ray, cell}} \\cdot \\Delta s\_{\\text{cell}}

where ell\\ell is the ray length inside a given cell. Writing this for every ray gives a linear system

G Δs=Δt\mathbf{G}\,\Delta\mathbf{s} = \Delta\mathbf{t}

where mathbfG\\mathbf{G} is the N_textraystimesN_textcellsN\_{\\text{rays}} \\times N\_{\\text{cells}} ray-path length matrix, Deltamathbfs\\Delta\\mathbf{s} is the slowness update vector, and Deltamathbft\\Delta\\mathbf{t} is the travel-time residual vector. This is the Amathbfx=mathbfbA\\mathbf{x} = \\mathbf{b} shape from Section 0.7, and it is mixed-determined: overdetermined in well-illuminated cells, underdetermined elsewhere.

2. The figure: fitting the times is not recovering the model

The figure uses a 16 by 16 grid of 100 m cells with a fast 4 by 4 block, 2700 m/s in 2000 m/s, at its centre. By default ten sources at the surface shoot to ten receivers along a line 1600 m deep, a transmission geometry like a crosswell survey turned on its side (the Well to well button shows the upright version), and give 100 straight rays. Coverage is uneven: centre cells are crossed by up to 15 rays, while 24 cells near the side edges are crossed by only one. The rays are straight and fixed, so mathbfG\\mathbf{G} never changes and each iteration is one step of a linear solver:

  1. Compute predicted travel times t_textpred=mathbfG,mathbfst\_{\\text{pred}} = \\mathbf{G}\\,\\mathbf{s} along the fixed straight rays through the current model.
  2. Compute residuals Deltat=t_textobs−t_textpred\\Delta t = t\_{\\text{obs}} - t\_{\\text{pred}}.
  3. Back-project with SIRT: each cell moves by Deltas_j=dfrac1N_jsum_iell_ij,dfracDeltat_isum_kell_ik2\\Delta s\_j = \\dfrac{1}{N\_j}\\sum\_i \\ell\_{ij}\\,\\dfrac{\\Delta t\_i}{\\sum\_k \\ell\_{ik}^2}, the average over the N_jN\_j rays that cross it of each ray’s residual, spread along the ray in proportion to its length in each cell.
  4. Repeat from the updated model.

Move the iteration slider, or run 20 iterations, and compare the recovered model in (b) with the truth in (a), the error map in (c) and the two curves in (e).

Tomography: rays through a velocity gridMany rays through many cells → solve for cell velocities (least-squares)

In 20 iterations the RMS residual falls from 35.7 ms to 1.47 ms, yet the model error falls only from 175 m/s to 129 m/s. Over its 16 cells the recovered block averages 2285 m/s against a true 2700 m/s, 41 % of its contrast, and it is smeared up and down along the rays, with fast lobes above and below it in (c). Fitting the travel times is not the same as recovering the model: watch the error, not only the residual. Add the well-to-well rays and the same 20 iterations recover 77 % of the contrast, and the model error falls to 66 m/s, because horizontal rays see the block from the side. Resolution comes from the spread of ray angles through each cell.

3. Straight-ray vs bent-ray

The widget uses straight rays: shot and receiver are connected by a line regardless of velocity. This is:

  • Exact for a constant-velocity model.
  • Approximately correct when velocity contrasts are small and refraction is mild.
  • Wrong when contrasts are strong, e.g., salt body, steep velocity gradient.

Production tomography uses bent rays: ray-trace through the current velocity model, using Snell’s law or the eikonal equation. At each iteration the rays are re-traced through the updated model, so mathbfG\\mathbf{G} itself changes and the problem becomes nonlinear. This is more accurate and more expensive.

4. Resolution and illumination

A cell with many rays crossing at diverse angles is well-illuminated: tomography resolves it crisply. A cell with only one or two rays at similar angles is poorly illuminated: the solution there is ambiguous, and the inversion either leaves it at the starting value or over-fits noise.

Common illumination gaps:

  • Below the deepest ray. No rays go below the deepest turning point; the model below is unconstrained.
  • Shadow zones. Strong velocity contrasts bend rays away from certain cells.
  • Survey edges. Corners and edges see only one-sided illumination.

5. Regularization

Because of the poorly-illuminated cells and noise in travel-time picks, the raw least-squares update is unstable. Tomography solves a regularized problem:

min_Deltamathbfs;∣mathbfG,Deltamathbfs−Deltamathbft∣_22+lambda2,∣mathbfL,Deltamathbfs∣_22\\min\_{\\Delta\\mathbf{s}}\\; \\|\\mathbf{G}\\,\\Delta\\mathbf{s} - \\Delta\\mathbf{t}\\|\_2^{2} + \\lambda^2\\,\\|\\mathbf{L}\\,\\Delta\\mathbf{s}\\|\_2^{2}

where mathbfL\\mathbf{L} is a smoothness operator (discrete Laplacian or gradient) and lambda\\lambda is the regularization weight. In the figure mathbfL=h,mathbfD\\mathbf{L} = h\\,\\mathbf{D}: mathbfD\\mathbf{D} takes first differences between neighbouring cells and hh = 100 m, so the penalty is lambda2h2∣mathbfD,Deltamathbfs∣_22\\lambda^2 h^2 \\|\\mathbf{D}\\,\\Delta\\mathbf{s}\\|\_2^{2} and its lambda\\lambda is a plain number. A low lambda\\lambda fits the data closely and gives a noisy-looking model; a high lambda\\lambda gives a smooth model that fits the data less well. In the figure, choose least squares with both surveys and 3 ms of pick noise. With lambda=0\\lambda = 0 the residual after 40 iterations is 0.30 ms, below the noise, and the model error is 65 m/s; lambda=1\\lambda = 1 halves the error to 30 m/s; lambda=3\\lambda = 3 leaves the residual at 4.83 ms, above the noise, and the error rises again to 58 m/s. Picking lambda\\lambda is the key art of tomographic inversion.

6. Inputs to production tomography

  • First-break picks: refraction tomography for near-surface velocity.
  • Reflection residuals: residual moveout in CMP or migrated gathers, projected through the overburden model.
  • Well-to-seismic ties: checkshot or VSP velocity constraints added as equations.
  • Walkaway VSPs: dense tomographic data around a well.

Most production velocity models come from combining refraction-tomography (near surface) and reflection tomography (deeper) iteratively with migration: the output of tomography feeds migration, and migrated gathers feed residual tomography. The loop closes after 2-4 iterations.

7. Tomography vs FWI

Tomography uses travel times only. FWI (Section 6) uses the full waveform: amplitudes, phases, everything. FWI can resolve finer structure but needs a good starting model (usually from tomography) or it falls into local minima. Production workflow: tomography first, then FWI started from the tomography model to add the detail travel times cannot resolve.

The one sentence to remember

Tomography solves G Δs=Δt\mathbf{G}\,\Delta\mathbf{s} = \Delta\mathbf{t} for the whole velocity model simultaneously from ray paths and travel-time residuals, with regularization to tame poorly-illuminated cells; it produces the velocity model every other velocity-dependent algorithm downstream depends on.

Part 3 closes here

You can now: compute RMS from interval velocities and invert via Dix, apply NMO and understand stretch, pick velocities from semblance, account for VTI with the η parameter, refine with residual velocity and higher-order moveout, and build a full depth-velocity model V(x,y,z)V(x, y, z) by tomography. Part 4 turns to multiple attenuation (SRME, Radon, adaptive subtraction), the final clean-up before imaging.

References

  • Yilmaz, Ö. (2001). Seismic Data Analysis (2 vols.). SEG.
  • Virieux, J., Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics, 74, WCC1.
  • Tarantola, A. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49, 1259.
  • Pratt, R. G. (1999). Seismic waveform inversion in the frequency domain, Part 1. Geophysics, 64, 888.
  • Etgen, J., Gray, S. H., Zhang, Y. (2009). An overview of depth imaging in exploration geophysics. Geophysics, 74, WCA5.
  • Bishop, T. N., et al. (1985). Tomographic determination of velocity and depth in laterally varying media. Geophysics, 50, 903.
  • Woodward, M. J., Nichols, D., Zdraveva, O., Whitfield, P., Johns, T. (2008). A decade of tomography. Geophysics, 73, VE5.

This page is prerendered for SEO and accessibility. The interactive widgets above hydrate on JavaScript load.