Tomography · Part A: the method, on a planted truth

Earthquakes and seismometers sit in rock whose speed we do not know. Can the arrival times tell us that speed, cell by cell — and which parts of the picture they draw can we believe?

Session 2, Location & relocation, held the velocity model fixed and solved for the earthquakes. At its end (A.11) it let the velocity of each flat layer become an unknown too, and the arrival times were enough to fix a handful of layer speeds. This session takes that step all the way: the unknown is the speed in every cell of a grid — hundreds or thousands of numbers. That is seismic tomography, and the picture it returns is what geologists read when they ask where a geothermal reservoir is, where the rock is fractured, or where it holds steam rather than water.

Two things change when the unknowns go from four to thousands.

• The derivative becomes the hard part. Every method in Sessions 1 and 2 needed the derivative of each predicted time with respect to each unknown. With thousands of unknowns, computing those derivatives one at a time is too slow. The answer is the adjoint, and it is the same idea as the backpropagation that trains a neural network (Session 5, Deep learning).
• The data no longer determine the answer. There are more cells than arrival times, and many cells are crossed by no ray at all. So the picture depends on what we add to the data — the regularisation — and we need tests that say which features to believe.

The questions this notebook answers

1. What is the forward problem — how does a velocity grid become arrival times? (A.1)
2. If the rays are straight lines, tomography is a linear problem. How is it solved, and what does it get wrong? (A.2–A.3)
3. How do we get the derivative of every arrival time with respect to every cell, cheaply? (A.4)
4. Given that derivative, how do we fit the model — Gauss–Newton, or the gradient descent that trains networks? (A.5)
5. The data cannot decide everything. What does the regularisation add, and how do we test what the picture can show? (A.6–A.7)
6. The earthquake locations are uncertain too. What happens when they are unknowns as well? (A.8)

Part A is synthetic, because accuracy needs a known answer. We plant a velocity model, compute the arrival times it would produce, add reading errors, and ask each method to find the model again. Every figure can then be scored against the truth. The world is two-dimensional — a vertical slice 24 km long and 10 km deep — so that every grid, ray and gradient can be drawn and every cell runs in seconds.

Part B is a separate notebook. It runs the same method in three dimensions on the 2016 Geysers picks, where nobody knows the answer.

A.1 · The forward problem: from a velocity grid to arrival times

The world

The slice is cut into square cells of side $h$ = 0.5 km. We store one number at every node (a corner of the grid): the slowness $f = 1/v$, in seconds per kilometre. Slowness rather than velocity, because travel time is slowness added up along the path: a path of length $\ell$ through rock of slowness $f$ takes $\ell f$ seconds.

The planted model has two parts:

• a 1D background $v_0(z)$ that grows with depth $z$, because rock deeper down is more compressed and faster. This is the kind of model Session 2 solved for.
• two anomalies: a slow body (red) and a fast body (blue). These are what tomography is for.

Twelve seismometers sit on the surface, and forty earthquakes lie between 2 and 8 km deep. We draw every model as the change from the 1D background, $\Delta v / v_0$, in per cent, with red for slow and blue for fast. That is the convention of seismic tomography figures.

Two kinds of wave. An earthquake sends out two body waves: the P (compressional) wave, which arrives first, and the S (shear) wave, which travels slower. Their speeds are $V_p$ and $V_s$, and their ratio $V_p/V_s$ is about 1.6–1.8 in crustal rock. The ratio is sensitive to what fills the pores — water raises it, steam or gas lowers it — which is why geothermal studies map it. The planted world here uses one wave type, so that every idea is seen once; Part B inverts P and S picks together, for $V_p$ and $V_s$.

grid spacing 0.5 km: 21 x 49 = 1029 nodes; 12 stations, 40 earthquakes
background speed 4.00 km/s at the surface to 6.50 km/s at 10 km

The same picture will be drawn many times, so one small function draws it: the velocity change in per cent, the stations as triangles, the earthquakes as dots.

The eikonal equation

How long does the first wave take to reach a point? Call that time $T(x, z)$. Picture the wave front: the set of points it has reached at one moment. In rock of slowness $f$ the front moves one kilometre in $f$ seconds. So the time grows fastest across the front, at the rate $f$ per kilometre:

$$ |\nabla T| \;=\; \sqrt{\Big(\frac{\partial T}{\partial x}\Big)^2 + \Big(\frac{\partial T}{\partial z}\Big)^2} \;=\; f(x, z). $$

This is the eikonal equation. Its solution with $T = 0$ at a source is the first-arrival time from that source to every point. It contains no rays — the ray is found afterwards, by walking downhill on $T$ — and that is why it handles any model, however strongly the rays bend.

Reciprocity saves work. The time from an earthquake to a station equals the time from the station to the earthquake. So we solve once per station, with the source at the station, and read every earthquake's time off that one map. Twelve solves give all 480 station–earthquake times.

Solving it on the grid

On the grid, replace each derivative by a difference with a neighbour. At node $(i, j)$ let

• $a$ be the smaller time of its two vertical neighbours (above and below), and
• $b$ the smaller time of its two horizontal neighbours (left and right).

The wave arrives from the earlier side, so those are the neighbours that matter (this choice is called upwind). The equation at the node becomes

$$ (T_{ij} - a)^2 + (T_{ij} - b)^2 \;=\; f_{ij}^2 h^2 . $$

Solve it for $T_{ij}$. Order the two so that $a \le b$. If the wave from $a$ alone arrives before $b$ — that is, if $a + f_{ij} h \le b$ — then only the vertical term counts and $T_{ij} = a + f_{ij} h$. Otherwise both count, and the quadratic gives

$$ T_{ij} \;=\; \tfrac12 \Big( a + b + \sqrt{2 f_{ij}^2 h^2 - (a - b)^2} \Big). $$

Every node depends on its neighbours, so we iterate: set $T = 0$ at the source and a large number elsewhere, update every node at once from its neighbours, keep the smaller of the old and new value, and repeat until nothing changes. The function below does that for all twelve stations together. Production codes (ADTomo, which Part B uses) visit the nodes in four alternating orders instead — the fast sweeping method (Zhao 2005, 10.1090/S0025-5718-04-01678-3) — and need fewer passes. Both solve the same equations and give the same answer.

12 travel-time maps in 6 ms

Here is one of the twelve maps: the time from the station at 11 km to every node, with the wave fronts as contours. Where the rock is fast the fronts are farther apart; through the slow body (red, left of the station) they lag behind and bunch together.

station 5 is at x = 11 km

How accurate is the grid?

In the background alone, where speed grows linearly with depth, the first-arrival time is known exactly: for $v = v_s + g z$, a source at the surface and a point at distance $r$ from it, $T = \frac{1}{g} \operatorname{arccosh}\!\big(1 + \frac{g^2 r^2}{2 v_s v}\big)$. Comparing the grid solution with it shows the price of the grid.

grid minus exact, below 2 km: mean 68 ms, largest 104 ms
as a fraction of the time: up to 9 %

spacing 0.500 km: grid minus exact, mean  68 ms, largest 104 ms
spacing 0.250 km: grid minus exact, mean  43 ms, largest  65 ms
spacing 0.125 km: grid minus exact, mean  27 ms, largest  39 ms

The grid is slow by several tens of milliseconds, up to nine per cent of the time. The four-neighbour stencil can only step along rows and columns, so it cannot follow a curved wave front exactly, and the error is born at the source, where the front is most sharply curved, and carried outward from there. The second cell shows it shrinks only slowly as the grid is refined: halving the spacing removes well under half of it.

That error is larger than a good pick's reading error, so it must not be ignored on real data: Part B removes it by computing the times in the starting model both on its grid and on a grid four times finer, and subtracting the difference from the data. Here we plant the data with the same solver that inverts them, so the grid error is in both and cancels. The planted test then measures the inversion alone. (Making synthetic data with the same code that inverts them is called an inverse crime; Part B does not commit it.)

The data

Each earthquake's time is read off each station's map by bilinear interpolation: a weighted average of the four nodes around the earthquake. Those weights form a sparse matrix $W$ (40 earthquakes × 1029 nodes), with four non-zeros per row. Then we add a reading error of 20 ms, typical of a good automatic pick. (The origin times are set to zero here; A.8 makes them unknown.)

480 arrival times, 0.63 to 4.45 s
residual in the 1D background: rms 45 ms (reading error 20 ms)

The background leaves residuals well above the reading error: the anomalies are in the data. The question for the rest of the notebook is how to turn those residuals back into the picture.

A.2 · Straight rays: tomography as a linear problem

What we solve for

We do not solve for the slowness itself but for the logarithm of its ratio to the background, one number per node:

$$ f \;=\; f_0\, e^{m}, \qquad \text{so} \qquad m \;=\; \ln\frac{f}{f_0} \;\approx\; \frac{f - f_0}{f_0} \;\approx\; -\frac{\Delta v}{v_0}. $$

Here $f_0 = 1/v_0$ is the background slowness and $m$ is a vector of $N$ = 1029 numbers. For a small change, $m$ is the fractional change of slowness: $m = 0.05$ is 5 % slower. A slow body has $m > 0$. Working with a fraction puts every node on the same scale, whatever its depth, and the exponential keeps the slowness positive whatever value $m$ takes.

The data we fit are the residuals $r$: the observed times minus the times in the background, $r = t^{\text{obs}} - t(f_0)$. One number per station–earthquake pair, 480 of them.

Linearise

For a small $m$, each residual is a weighted sum of the $m$ values:

$$ r_p \;\approx\; \sum_n J_{pn}\, m_n , \qquad \text{or} \qquad r \approx J\, m . $$

• $p$ counts the picks (station–earthquake pairs) and $n$ the nodes.
• $J_{pn} = \partial t_p / \partial m_n$ is the Jacobian: how much pick $p$'s time changes per unit change of $m_n$. It has 480 rows and 1029 columns.

This is exactly Session 2's linearisation, with 1029 unknowns instead of four.

The straight-ray Jacobian

Suppose each wave travels on the straight line from the earthquake to the station. Then its time is slowness added up along that line, $t_p = \int_{\text{line}} f\, d\ell$. Cut the line into short steps of length $\Delta\ell$. At each step the slowness is the bilinear average of the four nodes around it, so the time is a weighted sum of node slownesses:

$$ t_p \;=\; \sum_n G_{pn}\, f_n, \qquad G_{pn} = \text{the length of ray } p \text{ credited to node } n. $$

At the background ($m = 0$), changing $m_n$ changes $f_n$ at the rate $\partial f_n / \partial m_n = f_{0,n}$, so the Jacobian is $J_{pn} = G_{pn} f_{0,n}$: the time the ray spends near node $n$ in the background. This was the first form of seismic tomography (Aki & Lee 1976, 10.1029/JB081i023p04381). Most entries are zero, because each ray passes near only a few dozen of the 1029 nodes, so $G$ is stored as a sparse matrix.

G: 480 picks x 1029 nodes, 24004 non-zeros (4.9 %)
nodes no ray passes near: 389 of 1029

Two pictures of $G$. On the left, one row: the nodes a single ray touches, coloured by the length credited to each. On the right, the number of rays that pass near each node — the ray coverage. Below the deepest earthquakes and at the far corners no ray passes, so the data say nothing there at all.

More unknowns than data: damped least squares

There are 480 equations and 1029 unknowns, and hundreds of nodes appear in no equation. Least squares alone has infinitely many solutions. Session 1, Regression & uncertainty, met the same problem in a small form and fixed it with ridge regression: add a penalty that prefers small coefficients. Here we add two penalties.

• Damping prefers a small change from the background: the term $(m_n / d)^2$ for every node. The scale $d$ = 0.1 says a 10 % change at one node costs as much as one pick that misses by its reading error.
• Smoothing prefers neighbours to agree: the term $(\Delta m / (s\,h))^2$ for every pair of neighbouring nodes, where $\Delta m$ is their difference and $h$ the spacing. The scale $s$ = 0.03 per km says a change of 3 % per kilometre costs as much as one such pick.

Write $D$ for the matrix that takes all the neighbour differences per km. It has one row per pair of neighbouring nodes — every horizontal pair and every vertical pair, 1988 rows — and one column per node, 1029; each row holds a $+1/h$ and a $-1/h$. The whole objective is

$$ \Phi(m) \;=\; \Big\| \frac{r - J m}{\sigma} \Big\|^2 \;+\; \Big\| \frac{D m}{s} \Big\|^2 \;+\; \Big\| \frac{m}{d} \Big\|^2 . $$

Each term is a sum of squares, so $\Phi$ is the squared length of one long stacked vector, and minimising it is one ordinary least-squares problem with more rows:

$$ \begin{bmatrix} J/\sigma \\ D/s \\ I/d \end{bmatrix} m \;\approx\; \begin{bmatrix} r/\sigma \\ 0 \\ 0 \end{bmatrix}. $$

The damping rows are exactly Session 1's ridge. The scales $s$ and $d$ are hyperparameters: the data do not choose them, and A.6 shows what changing them does.

Solving it without forming $J^\top J$: LSQR

Session 1 solved least squares with the normal equations, $(J^\top J)\, m = J^\top r$. Here $J^\top J$ would be a dense 1029 × 1029 matrix; in three dimensions, with a million nodes, it could not even be stored. LSQR (Paige & Saunders 1982, 10.1145/355984.355989) is an iterative method that reaches the same answer using only two operations: multiply a vector by the matrix, and multiply a vector by its transpose. Keep that in mind: products with $J^\top$ will turn out to be exactly what the adjoint of A.4 delivers.

D: 1988 neighbour pairs x 1029 nodes
smoothing scale s = 0.03 per km, damping scale d = 0.1

Solve, and turn the answer back into a velocity change: $\Delta v / v_0 = e^{-m} - 1$. To score a model against the truth we use one number throughout: the correlation between the found and the true $\Delta v / v_0$, over the nodes inside the box the earthquakes span (2–8 km deep, 3–21 km along). 1 is a perfect pattern; 0 is no relation. It ignores amplitude, so we also print the largest slow and fast values.

straight rays: correlation 0.94, -6.4 to +4.9 %
misfit: rms 21.0 ms (reading error 20 ms)

Both bodies are found, smoothed and weaker than the truth. Damping pulls every node toward zero and smoothing spreads each body, so the recovered amplitude is lower. That is a general rule worth remembering for every tomographic image you will read: its amplitudes are lower bounds.

A.3 · Rays bend

A straight line is the path only in uniform rock. Our background speeds up with depth, so the quickest path spends more of its length in the faster rock below: it leaves the earthquake heading more sideways than the straight line, and turns to climb steeply under the station, where the rock is slowest. The true ray is found from the travel-time map: start at the earthquake and walk downhill on $T$, against its gradient, until you reach the station. The time decreases fastest along that path, so it is the path the first wave took.

Over these short distances the rays bend only a little, but enough to pass through different nodes than the straight line does. The straight-ray matrix credits the time to the wrong nodes. How much that matters depends on how far the rays travel and how strong the velocity contrasts are.

To use the true paths we need the true Jacobian: how the time computed by the eikonal solver changes when one node's slowness changes. That is the subject of the next section.

A.4 · The adjoint: the derivative of a solver

We want the Jacobian of the eikonal solver: how each arrival time changes when one node's slowness changes. There are three ways to get it.

Way 1: finite differences

Change one node's slowness by a tiny amount $\epsilon$, solve again, and divide the change in every arrival time by $\epsilon$. That gives one column of the Jacobian, so the full Jacobian needs one extra solve per node: here 1029 solves, and in a three-dimensional grid of a million nodes, a million solves per iteration.

one solve (12 stations): 5.9 ms, 80 passes over the grid
finite differences for every node: 6 s

A few seconds is fine for this slice, but the cost grows with the number of nodes, and it is needed again at every step of the fit. We keep finite differences for one job only: checking the faster method on a few nodes.

Way 2: automatic differentiation through the iterations

The solver is a loop of ordinary arithmetic, so a library like PyTorch could record every operation and run the chain rule backwards through them (reverse-mode automatic differentiation, Session 5). It would work, but it would store the travel-time maps of every pass — the number of passes printed above, and many more in 3D — and it would differentiate the path the iteration took, when only the answer it reached matters.

Way 3: the adjoint — differentiate the equations, not the iterations

The iteration does not matter; only its result does. The result satisfies one equation per node. Write it as a residual that is zero at the solution:

$$ R_n(T, f) \;=\; (T_n - a_n)^2 + (T_n - b_n)^2 - f_n^2 h^2 \;=\; 0 \quad \text{for every node } n \text{ except the source,} $$

where a term is kept only if that neighbour is upwind (A.1), and $a_n$, $b_n$ are the times of the upwind neighbours. At the source the equation is simply $T = 0$, which does not depend on $f$. Now follow three steps.

Step 1 — perturb. Change the slowness by a small $\delta f$. The times change by some $\delta T$, and every equation must still hold. To first order:

$$ \sum_k \frac{\partial R_n}{\partial T_k}\,\delta T_k \;+\; \frac{\partial R_n}{\partial f_n}\,\delta f_n \;=\; 0 . $$

Step 2 — name the pieces. Call $A$ the matrix of $\partial R_n / \partial T_k$. Row $n$ has at most three non-zeros: $2(T_n - a_n) + 2(T_n - b_n)$ on the diagonal, and $-2(T_n - a_n)$ and $-2(T_n - b_n)$ at the two upwind neighbours. And $\partial R_n / \partial f_n = -2 f_n h^2$. So

$$ A\, \delta T \;=\; 2 f h^2 \odot \delta f , $$

where $\odot$ multiplies element by element (and the source's entry on the right is zero). The derivative of every time with respect to every slowness is therefore $\partial T / \partial f = A^{-1}\, \mathrm{diag}(2 f h^2)$.

Step 3 — ask for what we need. We rarely need the whole Jacobian. Usually we need the gradient of one number — the misfit $L$ — with respect to all slownesses. The chain rule gives

$$ \frac{\partial L}{\partial f} \;=\; \Big(\frac{\partial T}{\partial f}\Big)^{\!\top} \frac{\partial L}{\partial T} \;=\; 2 f h^2 \odot \lambda, \qquad \text{where } \lambda \text{ solves } A^\top \lambda = \frac{\partial L}{\partial T} . $$

The vector $\lambda$ is the adjoint field. Getting it costs one sparse linear solve per station, whatever the number of nodes. That is the whole trick.

Two remarks make it concrete.

• $A$ is easy to solve. Each node depends only on neighbours with earlier times, so sorted by time, $A$ is triangular, and back-substitution solves it in one pass. (A sparse solver finds that on its own.) In the adjoint the order runs backwards, from the latest times to the source.
• This is reverse-mode automatic differentiation. The forward pass computes $T$; the backward pass receives $\partial L / \partial T$ and returns $\partial L / \partial f$. That is exactly the job of a backward function in PyTorch. So we write the solver as one custom torch.autograd.Function, and PyTorch's own automatic differentiation does everything else (the interpolation, the misfit, the penalties), as backpropagation did for the networks of Session 5.

Now the solver and its adjoint become one PyTorch operation. forward runs the eikonal solver of A.1; backward runs the adjoint above. Nothing else is needed: PyTorch calls backward whenever a gradient has to pass through the solver.

Check it against finite differences

A derivative written by hand must be checked. Compute the gradient of the misfit at the background model by the adjoint, then compare a few of its entries with finite differences, which need nothing but two extra solves each.

misfit and its gradient at all 1029 nodes: 34 ms, 5.7 times one solve
node at x  7.0 km, z 3.0 km:  adjoint     -3458.54   finite difference     -3458.54
node at x 16.0 km, z 4.0 km:  adjoint      1703.41   finite difference      1703.41
node at x 10.0 km, z 1.5 km:  adjoint      -873.81   finite difference      -873.81
node at x 12.0 km, z 6.0 km:  adjoint      -299.19   finite difference      -299.19
node at x  7.0 km, z 0.0 km:  adjoint      -274.79   finite difference      -274.79

They agree to every printed digit: the adjoint gives the exact derivative of the discrete solver, all 1029 entries together for the cost of a handful of solves — the forward solve plus twelve sparse linear solves, one per station — whatever the number of nodes. The last node is a station. In its own map the time there is fixed at zero, but the other eleven maps pass through it, and the two methods agree there too.

What one arrival time is sensitive to

Take a single pick, station 5 to earthquake 3, and ask for the gradient of its time alone: set $\partial L / \partial T$ to that earthquake's interpolation weights in station 5's map and zero everywhere else. Multiplied by $f_0$, the answer is one row of the true Jacobian. It is called the pick's sensitivity kernel. Next to it is the same row from the straight-ray matrix of A.2.

kernel sums: bent 1.666 s, straight 1.638 s
times in the background: exact 1.636 s, straight line 1.638 s, eikonal solver 1.666 s

Both show station 5 to earthquake 3 (the star). The adjoint's kernel follows the bent ray: it rises steeply under the station, as the traced rays of A.3 did. It is also not a thin line but a fan that widens away from the station, about two kilometres across at the earthquake. The four-neighbour grid passes time only along rows and columns, so the time at a node depends on a spread of nodes between the vertical and the diagonal. The kernel is the exact derivative of the solver we use, grid and all.

Each kernel sums to its own travel time. Raising every slowness by the same small fraction raises the time by that fraction, whatever the path, so the kernel adds up to the time itself.

The three printed times say what each approximation costs. The exact first arrival is the fastest: by Fermat's principle the true ray is the path of least time, so no other path — the straight line included — can be faster. Here the bent ray wins by only a few milliseconds, because the rays bend little. The grid solver's time is slower than both, by tens of milliseconds: that difference is the grid error of A.1, not the path.

The gradient of the whole misfit

The gradient we checked above is the sum of all 480 kernels, each weighted by its residual. It points, at every node, in the direction that increases the misfit fastest. So its negative, drawn below, is the first step a gradient-descent method would take from the background.

The first gradient already shows the two bodies in place, but smeared along the rays and strongest near the stations, where every ray passes. Turning that direction into a model is the job of the fitting methods in A.5.

A.5 · Fitting the model

We now have everything needed to fit: the forward solver, its exact derivative, and the regularised objective of A.2. There are two ways to use the derivative. Both appeared before, Gauss–Newton in Session 2 and gradient descent with Adam in Session 5.

Gauss–Newton: relinearise and solve, again and again

Straight rays gave one linear problem, solved once. With bent rays the Jacobian depends on the model, because the rays move when the model changes. So we repeat:

1. solve the eikonal equation in the current model and compute the residuals $r$;
2. compute the Jacobian $J$ in the current model;
3. solve the damped least-squares problem of A.2 for a step $\delta m$, with the penalties on the total model $m + \delta m$;
4. take the step and go back to 1.

This is Session 2's Gauss–Newton (Geiger's method), with 1029 unknowns. Step 2 needs the whole Jacobian, one row per pick. The adjoint gives a row for one right-hand side $\partial L/\partial T$; we need 40 per station, one per earthquake. The matrix $A$ is the same for all 40, so we factor it once per station and reuse the factors: twelve factorisations for the whole Jacobian.

Gauss-Newton, 8 iterations in 0.3 s
rms misfit per iteration (ms): 45.3 20.6 20.6 20.6 20.6 20.6 20.6 20.6 20.6
bent rays:     correlation 0.97, -6.7 to +5.6 %
straight rays: correlation 0.94, -6.4 to +4.9 %

The misfit drops to the reading error in the first iteration and then stays: the rays bend only slightly, so the problem is nearly linear. Compare the two models with the truth. The straight-ray model puts the slow body too shallow and too weak, because it credits each delay to the nodes on the straight line rather than on the ray. The bent-ray model puts both bodies at their true depths, with a higher correlation. In this gentle slice straight rays already do fairly well; the difference grows with the velocity contrast and with the distance the waves travel.

Gradient descent with Adam, and early stopping

Gauss–Newton needs the whole Jacobian at every step. For a 3D model with a million nodes and a million picks that matrix is enormous. Session 5 trained networks with millions of weights using only the gradient of the loss, which the backward pass gives at the cost of about one more forward pass. The same works here: the Eikonal function of A.4 makes the misfit an ordinary PyTorch loss, and Adam can minimise it.

We run Adam with no regularisation at all, to see what the optimiser does on its own. A random fifth of the picks is held out as a validation set, as in Session 5: the model never fits them, and their misfit tells us when it starts fitting noise.

300 Adam steps in 4 s
validation misfit is lowest at step 68: 1.07 (training 0.78)
at step 68: correlation 0.63, -12.7 to +20.2 %
the correlation with the truth is highest at step 24: 0.69
at step 299: correlation 0.34, -39.8 to +188.8 %

The training misfit falls all the way. The held-out misfit falls, reaches a minimum, and rises again: after that point the model is fitting the reading errors of the training picks, which do not repeat in the held-out ones. The correlation with the truth (blue) peaks earlier still and then falls steadily. The held-out picks share their stations and earthquakes with the training picks, so they notice overfitting late: a useful guide, not a perfect one. Spikes appear first at the nodes just below the deepest earthquakes, which only a few picks touch, and by the end they cover the model.

Stopping at the validation minimum is early stopping, and it acts as a regulariser. Gradient descent changes first what the data constrain strongly — large regions crossed by many rays — and only later the fine detail that a few rays touch. Stopping early leaves that detail at the background, much as damping does. But it is a weaker and less controlled regulariser than the penalties of A.2: the stopped model is rougher than the Gauss–Newton model, and nothing in it says which structure was left out.

Now give Adam the same penalties Gauss–Newton used. Then both minimise the same objective, and should find the same model.

Adam + penalties, step 299: correlation 0.97, -6.7 to +5.6 %
Gauss-Newton:               correlation 0.97, -6.7 to +5.6 %
largest difference between the two models: 0.02 %

With the penalties in the objective, the optimiser no longer decides what the model looks like: Adam arrives at the Gauss–Newton model. It needed 300 steps where Gauss–Newton reached the noise level in one iteration (A.5), because Adam uses only the gradient while Gauss–Newton also uses the curvature. The objective decides the answer; the optimiser decides how fast you get there. In the rest of this notebook we use Gauss–Newton, because it converges in a few iterations on a problem this size. Gradient methods (Adam, or L-BFGS, which estimates the curvature from past gradients) take over when the Jacobian is too large to store.

A.6 · The regularisation decides part of the picture

The scales $s$ and $d$ of A.2 were chosen, not measured. Scale both penalties together by one strength factor and fit the same data three times: with the setting used so far, three times weaker, and three times stronger. For each, print the misfit next to the reading error and the correlation with the truth. Print one more number, which A.7 explains in full: the false structure — the rms $\Delta v / v_0$ the same settings find when the planted truth has no anomalies at all.

strength x0.33: misfit 18.6 ms,  correlation 0.97, -9.4 to +8.3 %,  false structure 0.81 %

strength x1.00: misfit 20.6 ms,  correlation 0.97, -6.7 to +5.6 %,  false structure 0.28 %

strength x3.00: misfit 26.8 ms,  correlation 0.89, -3.3 to +2.5 %,  false structure 0.09 %

Three different pictures, from the same data.

• Weaker penalties fit the data better than the reading error of 20 ms allows. The extra fit is noise, drawn into the model as blobs that are not in the truth; the null test shows it, with about three times more false structure.
• Stronger penalties give a smoother, weaker image that no longer fits the data to the reading error: real signal is left in the residuals, and the amplitudes fall by half.
• The middle setting fits the data at about the reading error. That is the usual rule for choosing the strength: the strongest regularisation that still fits the data as well as the noise allows (the discrepancy principle). It needs an honest estimate of the reading error.

Session 1 chose ridge's strength by cross-validation; that works here too, with held-out picks as in A.5. Either way, the image you see is the data and a choice, and a paper that shows a tomographic image should say what the choice was.

A.7 · How we know what to believe: planted tests

On real data there is no truth to compare against. But the geometry — where the stations and the earthquakes are — is known, and the geometry decides what any inversion can see. So we plant a truth of our choice in the real geometry, make its data with the same reading error, invert them with exactly the same settings, and compare. Three kinds of truth answer three questions.

1. The null test — the background alone, no structure at all. Everything the inversion draws is false. Its size is the level a real feature must clearly exceed.
2. The anomaly test — realistic structure: the two bodies of A.1. You have already run it: the fit of A.5 is this test. It asks whether such structures come back, and how strongly.
3. Checkerboards — cells of ±6 %, 4 km wide and 4 km tall, and a thin one with cells only 1.5 km tall. Where the pattern comes back, the data can resolve features of that size; the thin one tests resolution in depth.

null test                   found rms 0.28 %

checker 4 km x 4 km         found rms 1.99 %, correlation 0.65, -5.0 to +5.4 %

thin checker 4 km x 1.5 km  found rms 0.81 %, correlation 0.12, -3.6 to +4.3 %
anomaly test (A.5)          found rms 2.74 %

Read the three columns in turn.

• The null test comes back almost blank. Its false structure, a fraction of a per cent, is what the reading errors alone produce with these settings. The anomaly test (A.5), measured the same way (rms over the same box), is about ten times larger, so its anomalies are not noise.
• The 4 km checkerboard comes back in its top row of cells, 0–4 km deep, where rays from every earthquake cross under the stations — at about the right places and with lower amplitude. The second row, 4–8 km, comes back only in part, and nothing comes back below the earthquakes.
• The thin checkerboard is lost, apart from a trace of its top layer just under the stations. The rays from deep earthquakes to surface stations are nearly vertical, so each one crosses several fast and slow layers, and their effects cancel in the travel time. Rays that all run in one direction cannot separate structure along that direction. This is the commonest blind spot of local earthquake tomography, where every station sits on top.

A planted test answers "could this geometry see a feature like this one?". It does not prove a real feature is right, and a checkerboard that comes back can still hide smearing along the rays. So use several tests, each shaped like the question you are asking. Part B does that for the Geysers geometry, with a test truth built to look like the real result.

Optional: the resolution matrix and the L-curve

Two tools appear in almost every tomography paper, and they answer the two questions of A.6 and A.7 without planting anything.

• The resolution matrix. For the linear problem of A.2 the estimate is a linear function of the data, and so of the true model: $\hat m = R\, m_{\text{true}}$ + (a noise term), with $R = \big(J^\top J/\sigma^2 + D^\top D/s^2 + I/d^2\big)^{-1} J^\top J/\sigma^2$, an $N \times N$ matrix. Row $n$ of $R$ shows how the estimate at node $n$ averages the true model: a sharp spike at $n$ means the node is resolved, a broad row shows the smearing, and a small diagonal means the node is barely resolved at all. A checkerboard test is $R$ applied to one particular truth. For millions of nodes $R$ is expensive, which is why 3D studies more often show planted tests.
• The L-curve (Hansen 1992, SIAM Review, doi:10.1137/1034115). Fit the data at many strengths and plot, on log axes, the misfit against the size of the penalty term. The points trace an L: weak regularisation fits the noise at a large penalty, strong regularisation pays in misfit. The corner is a common choice of strength when the reading error is not known well enough for the discrepancy principle.

A.8 · The earthquakes are unknowns too

So far the earthquakes sat at their true positions. On real data they do not: they were located in a velocity model, and that model is the thing we are trying to correct. An earthquake has three unknowns in this slice — $x$, $z$ and its origin time $t_0$ — so forty earthquakes add 120.

Locating in the background

First do what a catalogue does: locate every earthquake in the background model, starting from the middle of the slice. This is Session 2's Gauss–Newton (Geiger's method), one earthquake at a time. The derivatives of a time with respect to the earthquake's position are the slopes of the station's travel-time map at the earthquake, $\partial T / \partial x$ and $\partial T / \partial z$, and the derivative with respect to $t_0$ is 1. Each step is capped at 1 km, so a bad first guess cannot throw an earthquake out of the slice.

120 location unknowns: x, z and t0 for 40 earthquakes
located in the background: error mean 415 m, largest 1298 m
located in the true model:  error mean 108 m (the reading error alone)

The wrong velocity moves the earthquakes by hundreds of metres, several times more than the reading error alone does.

Fix the earthquakes where the background put them?

The simplest plan: keep those locations and origin times, and invert for the velocity as in A.5. The trouble is that part of what the residuals say has already been absorbed. Each earthquake moved to wherever the wrong model fits best, and that move explained away some of the travel-time anomaly the velocity should have explained.

events fixed at the background locations: correlation 0.57, -2.4 to +1.7 %
events fixed at the truth (A.5):            correlation 0.97, -6.7 to +5.6 %

Solve for both together — with parameter separation

So let the earthquakes move during the fit too. The difficulty: velocity and location trade off. The waves of an earthquake inside a slow body arrive late at every station. A later origin time explains much of that delay as well, and so does a move of the earthquake. One delay has two possible causes, and the fit must decide how much goes to the velocity and how much to the earthquake.

Take one earthquake $e$ and its twelve residuals $r_e$, one per station. Linearised, a velocity change $\delta m$ and a hypocentre change $\delta h_e = (\delta x, \delta z, \delta t_0)$ explain them together:

$$ r_e \;\approx\; J_e\,\delta m \;+\; H_e\,\delta h_e . $$

• $J_e$ is this earthquake's twelve rows of the Jacobian of A.3 (12 × 1029).
• $H_e$ is 12 × 3: how each of the twelve times changes when the earthquake moves in $x$, moves in $z$, or shifts its origin time. Its columns are $\partial t/\partial x$, $\partial t/\partial z$ and 1, as in Session 2.

The three columns of $H_e$ give every pattern of residuals that a move of this earthquake can make. Parameter separation (Pavlis & Booker 1980, 10.1029/JB085iB09p04801) gives all those patterns to the earthquake, and fits the velocity only to the rest.

The tool is a projection:

$$ Q_e \;=\; I - H_e \big(H_e^\top H_e\big)^{-1} H_e^\top . $$

$Q_e$ removes from a 12-vector the part that the columns of $H_e$ can make, and keeps the rest. So $Q_e H_e = 0$. Multiply the equation above by $Q_e$, and the earthquake's unknowns drop out:

$$ Q_e\, r_e \;\approx\; Q_e J_e\,\delta m . $$

A simple case shows what $Q_e$ does. Keep only the column of ones (the origin time). Then $Q_e$ subtracts the mean of the twelve residuals. A delay that is the same at every station is never blamed on the velocity, because a later origin time explains it just as well. The two position columns do the same for the patterns that a move of the earthquake makes.

Each iteration has three steps.

1. Locate every earthquake in the current velocity model (locate above).
2. Project: multiply each earthquake's residuals and Jacobian rows by its $Q_e$.
3. Step: take the velocity step from the projected system, with the penalties of A.2. Go back to 1.

The velocity step is the same as if we solved for the 1029 velocity nodes and the 120 location unknowns at once, but the system holds only the velocity. It is the same coupled problem as Session 2's A.11, where the unknown was a handful of layer speeds rather than 1029 nodes.

joint fit in 1 s
joint, parameter separation: correlation 0.93, -4.6 to +2.1 %
location error: background 415 m -> joint 211 m (mean)

The grey lines join each found location to the true one. Located in the background, the deeper earthquakes were pushed deeper still, by up to a kilometre. Fixed there, they keep those errors, and the velocity image keeps only about a third of its amplitude, blurred: the moves had already absorbed most of the anomaly.

The joint fit recovers much of both. The slow body comes back close to the fixed-at-truth image, the fast body weaker, and the mean location error halves. Fixing earthquakes where the wrong model put them keeps the pattern — the right places and signs — but loses most of the amplitude; solving for them together with the velocity recovers most of it.

Two limits remain, and Part B meets both on real data. Near an earthquake, velocity and location can never be fully separated, so anomalies at the earthquakes are the least certain part of any image. And an error shared by all earthquakes — every S pick late by the same amount, say — looks exactly like a slower model with deeper earthquakes, which no amount of data can tell apart.

A.9 · Optional: the adjoint-state method without the grid

A.4 derived the adjoint for the grid's equations. The same three steps work on the eikonal equation itself, and show what the adjoint field $\lambda$ is.

Write the eikonal equation as $|\nabla T|^2 = f^2$, and the misfit as $L = \tfrac12 \sum_p \big(T(\mathbf{x}_p) - t_p\big)^2$, a sum over the receivers $\mathbf{x}_p$ of one source. Call $\rho_p = T(\mathbf{x}_p) - t_p$ the residual at receiver $p$.

Perturb. Change $f$ by $\delta f$; then $T$ changes by $\delta T$ and the misfit by

$$ \delta L = \sum_p \rho_p\, \delta T(\mathbf{x}_p). $$

The eikonal equation must still hold, so to first order $2\,\nabla T \cdot \nabla \delta T = 2 f\, \delta f$, or

$$ \nabla T \cdot \nabla \delta T \;=\; f\, \delta f . $$

Multiply by a field $\lambda(\mathbf{x})$ and integrate over the whole domain:

$$ \int \lambda\, \nabla T \cdot \nabla \delta T \; d\mathbf{x} \;=\; \int \lambda f\, \delta f \; d\mathbf{x}. $$

Integrate the left side by parts, moving the derivative from $\delta T$ onto $\lambda \nabla T$ (the boundary term vanishes with the right condition on $\lambda$ at the edges):

$$ -\int \delta T\; \nabla \cdot (\lambda \nabla T) \; d\mathbf{x} \;=\; \int \lambda f\, \delta f \; d\mathbf{x}. $$

Choose $\lambda$ so that the left side is exactly $\delta L$. That requires

$$ -\nabla \cdot (\lambda\, \nabla T) \;=\; \sum_p \rho_p\, \delta(\mathbf{x} - \mathbf{x}_p), $$

the adjoint equation: a point source of strength $\rho_p$ at each receiver. Then

$$ \delta L \;=\; \int \lambda f\, \delta f\, d\mathbf{x}, \qquad\text{so}\qquad \frac{\partial L}{\partial f(\mathbf{x})} \;=\; \lambda(\mathbf{x})\, f(\mathbf{x}). $$

The adjoint equation is a transport equation along $\nabla T$, that is, along the rays. It carries each receiver's residual backwards along the rays to the source, and $\lambda$ is non-zero only where the rays went. The discrete $A^\top \lambda = \partial L/\partial T$ of A.4 is the same equation on the grid, solved from the latest times back to the source. This is the method of Leung & Qian (2006) and Taillandier et al. (2009), and the reason "adjoint" and "backpropagation" name one idea: the forward pass carries information from the source outwards, and the backward pass carries the error back the same way it came.

A.10 · Where a network enters

Session 5 showed that a network is a parameterised function fitted by the gradient of a loss. Two places in tomography use that idea directly.

• The model as a network. Instead of one number per node, write $f(\mathbf{x})$ as a small network of $\mathbf{x}$ and fit its weights. Nothing above changes: the adjoint gives $\partial L / \partial f$ at the nodes, and backpropagation carries it on to the weights. The network's shape then plays the role of the smoothing penalty.
• The solver as a network. EikoNet (Smith, Azizzadenesheli & Ross 2021, IEEE TGRS 59(12), 10685–10696) trains a network $T_\theta(\mathbf{x}_s, \mathbf{x})$ to satisfy the eikonal equation itself, with the loss $\big(|\nabla_{\mathbf{x}} T_\theta| - f(\mathbf{x})\big)^2$ at random points. It needs no grid, so it has no grid error, and its derivatives come from automatic differentiation for free. The price is training time and an accuracy that is only as good as the training. It is a physics-informed neural network, and the idea behind the location method of Session 2's HypoSVI.

The programs the field runs

• SIMULPS and its successors (Thurber 1983; Eberhart-Phillips 1986) made most local earthquake tomography of the 1980s–2000s: approximate ray bending, and velocity and hypocentres together with parameter separation, as in A.8.
• tomoDD (Zhang & Thurber 2003) adds differential times between nearby earthquakes, which sharpen the image where the earthquakes cluster.
• FMTOMO (Rawlinson & Sambridge) computes travel times with the fast marching method, an eikonal solver like ours, and is common in regional and teleseismic tomography.
• ADTomo (AI4EPS) is the method of this notebook in three dimensions: an eikonal solver with an adjoint written as a PyTorch function, Gauss–Newton with parameter separation, and L-BFGS as a backup. Part B runs it on the Geysers.

Takeaways

Machine learning

• The adjoint method gives the gradient of a misfit with respect to every model parameter for about the cost of one more forward solve, by solving the transposed linearised equations instead of differentiating through the solver's iterations.
• That is reverse-mode automatic differentiation: a solver with a hand-written adjoint becomes one torch.autograd.Function, and the library differentiates everything around it.
• A hand-written derivative is trusted only after it matches finite differences on a few entries, including the special ones such as a source.
• Gradient descent with early stopping regularises, but weakly and implicitly; with an explicit penalty the optimiser no longer decides the answer, only how fast it is reached.
• When there are more unknowns than data, the regularisation is part of the answer. Choose its strength so the data are fitted to their noise level, and measure what it invents with a test whose truth has no structure.

Seismology

• Travel-time tomography fits a velocity grid to arrival times. Straight rays make it linear; the eikonal equation handles bent rays, and its derivative comes from the adjoint.
• A tomographic image is smoothed and its amplitudes are lower bounds, because damping and smoothing pull toward the background.
• What an image can show is decided by the geometry: planted tests in the real geometry — a null test, checkerboards, structures shaped like the result — say which features to believe.
• With stations only at the surface, rays are steep, and structure layered in depth cancels along them: depth resolution is the common blind spot.
• Earthquakes located in the wrong model have already absorbed part of the anomaly; velocity and locations must be solved together, with parameter separation so the velocity fits only what moving the earthquakes cannot explain.

Part B runs the same method in three dimensions on the 2016 Geysers picks, with P and S waves, and asks which features of the field's velocity structure survive the tests of A.7.