Reading Rawlinson and Spakman (2016) on resolution tests in seismic tomography

Open N. Rawlinson and W. Spakman, “On the use of sensitivity tests in seismic tomography”, Geophysical Journal International 205, 1221–1243 (2016), next to this page. An open copy is in the Utrecht University repository. Section 2 links synthetic tests to the resolution matrix in a few lines of algebra. Section 3 then runs fourteen figures of experiments on one synthetic data set. By the end you should be able to say, for any published resolution test, which part of the resolution matrix it samples and which common mistake it might contain.

Sandbox Physics · 2026-10-30

Section 2.1 and 2.2: three equations to keep open

Eq. 6 is the linearized system, data equal the matrix times the model plus an error term; Section 4.19 of the MIT 12.201 notes writes the same system for straight rays, with path lengths in each block. Read the sentence after Eq. 6: the error collects picking noise, error of the forward theory, linearization error and the error of the chosen parametrization. Most later arguments are about which of these terms a given test can see.

Eq. 7 is the cost function, data misfit plus a penalty weighted by a trade-off parameter, and Eq. 9 is the generalized inverse that minimizes it. The common misreading is to treat resolution as a property of the rays alone. Eq. 9 shows the inverse depends on the data covariance and on the regularization, so the resolution matrix of Eq. 13 changes whenever the damping changes.

Eq. 13 splits the solution into the resolution matrix times the true model plus the inverse applied to the error. Eq. 14 gives the model covariance from the data covariance alone. Keep the two terms apart: resolution says how the truth is blurred, covariance says how noise is spread.

\mathbf d=\mathbf G\mathbf m+\boldsymbol\epsilon\ \ (6),\qquad \mathbf G^{-g}=\big(\mathbf G^{\mathsf T}\mathbf C_d^{-1}\mathbf G+\alpha^2\mathbf D^{\mathsf T}\mathbf D\big)^{-1}\mathbf G^{\mathsf T}\mathbf C_d^{-1}\ \ (9),\qquad \hat{\mathbf m}=\mathbf R\,\mathbf m+\mathbf G^{-g}\boldsymbol\epsilon\ \ (13)

Section 2.3: why a synthetic test samples the resolution matrix

Eqs. 15 and 16 are the heart of the paper. Make synthetic data from a test model with the same matrix, invert with the same generalized inverse, and the output is the same resolution matrix times the test model, plus the inverse applied to any synthetic noise.

So a noise-free test measures resolution and nothing else, a noise-only test measures error propagation, and a noisy test mixes the two. When a paper shows a noisy checkerboard, check which question it meant to answer.

\mathbf d_s=\mathbf G\mathbf m_s+\boldsymbol\epsilon_s\ \ (15),\qquad \hat{\mathbf m}_s=\mathbf R\,\mathbf m_s+\mathbf G^{-g}\boldsymbol\epsilon_s\ \ (16)

Section 2.5: four of the eight implications to mark

Point (ii): a test model in the null space is not recovered at all, and one in the row space is recovered perfectly. Here the authors cite Lévêque and co-authors (1993), whose idealized experiment recovered a fine checkerboard perfectly while a coarse one, with the same paths, was almost entirely in the null space. A good checkerboard therefore cannot prove good resolution at other scales; tests can only expose a lack of it.

Point (iii): a single spike returns one row and column of the resolution matrix, and well-separated spikes return many at once. Point (iv): the test must use the matrix of the real inversion. Tracing new rays through the test model gives a different matrix and tests something else. Point (viii): using the recovered model itself as a test input gives back nearly the same model when the resolution matrix is close to idempotent, so such a test hides artefacts.

One reading note: Section 3.6 refers to this last argument as point vii, but in the list it is point (viii). Point (vii) is the superposition argument of Eq. 19.

Figures 1 to 5: checkerboard against spikes

Figure 1 sets up every later experiment. Panel (a) is the true model, Gaussian anomalies up to plus or minus 20 per cent; (b) the rays through it; (c) the model recovered after 20 nonlinear iterations, with the RMS travel-time residual falling from 18.4 to 0.35 s; (d) the rays through the recovered model. Coverage is good in the north-west and degrades to the south-east, where it becomes nearly one-directional. All tests in Figures 2 to 5 reuse the rays of panel (d).

Figure 2 is a fine checkerboard and Figure 3 is a set of sparse spikes of the same size and amplitude. Compare them side by side. The checkerboard suggests the south-east is simply unresolved; the spikes show that most anomalies there are partly recovered but stretched, with aspect ratios up to 5 or 6 to 1, and show the direction of the stretch. The misreading to avoid is to treat a blank checkerboard region as no information.

Figure 5 checks the spike test against the real resolution matrix. With 3844 unknowns the authors could compute the matrix by Cholesky decomposition in about 12 minutes; the spike test took about 8 seconds, and its output is essentially the same.

Figures 6 to 14: nonlinearity, scale, noise and display

Figures 6 and 7 trace rays through the test models instead of reusing the observed rays. First arrivals bend towards fast regions, the coverage changes and the checkerboard smears far more along the diagonal. This is point (iv) made visible. Figure 8 estimates the nonlinear error term and shows that halving the anomaly amplitudes reduces it.

Figure 9 repeats the tests with anomalies about four times larger, and recovery improves: in this set-up larger structures are easier to resolve. Figure 10 feeds the recovered model back in and gets it back almost unchanged, the warning of point (viii). Figure 11 varies damping and smoothing; a synthetic test must use the regularization of the real inversion. Figure 12 adds Gaussian noise of 1 and 5 s; without regularization the 5 s case overfits and invents structure. Figure 13 shows the error from a parametrization that cannot represent the input, and Figure 14 shows that a colour scale alone can change the impression of resolving power.

The Summary ends with six recommendations: use a sparse set of discrete spikes in place of the usual tight checkerboard; reuse the coverage of the observation-based model; test several scale lengths; be careful with test models that mimic the recovered one; test with realistic noise; and, where possible, assess the effect of the parametrization error.

Worked example

A damped four-cell grid is not idempotent

Point (viii) depends on how close the resolution matrix is to idempotent, that is, on whether applying it twice equals applying it once. Take the four-cell grid of chapter six, with two row rays and two column rays of 1 km per cell, unit data errors and damping with a weight of one.

The matrix of the data term has eigenvalues 4, 2, 2 and 0. With damping each eigenvalue becomes itself over itself plus one, so the resolution matrix has eigenvalues 0.8, 0.667, 0.667 and 0, and every diagonal entry is 0.533. A matrix is idempotent only if its eigenvalues are all 0 or 1, so this one is not.

Run a spike of one unit in the top-left cell. The recovered image is 0.533 there, 0.200 in each neighbour and minus 0.133 in the opposite corner. Use that image as a new test input, as in Figure 10. The output is 0.382, 0.160, 0.160 and minus 0.062: it does not return the input. The near-complete recovery in the paper’s Figure 10 reflects the authors’ set-up and regularization, which they describe as at least quasi-idempotent. The damped four-cell grid shows how far another set-up can sit from that.

\mathbf R=\big(\mathbf G^{\mathsf T}\mathbf G+\mathbf I\big)^{-1}\mathbf G^{\mathsf T}\mathbf G:\ \ \mathrm{eig}=\tfrac45,\tfrac23,\tfrac23,0;\qquad \mathbf R\,\mathbf e_1=(0.533,\,0.200,\,0.200,\,-0.133),\quad \mathbf R^2\mathbf e_1=(0.382,\,0.160,\,0.160,\,-0.062)

A checklist for the next tomography paper

Find the test models: checkerboard, spikes, or the recovered model itself, and at how many scale lengths. Check that the synthetic data used the rays and the matrix of the real inversion, and the same damping, smoothing and parametrization. Note whether the test was noise-free, noise-only or noisy, and what noise level. Look at the colour scale and where it saturates. Finally, ask whether the paper shows resolution and model uncertainty separately, since each needs its own test.

An experiment to think through

Predict, change, explain

In the tomography Lab, set the true model to Checkerboard and then to Spike, with Stations all around. Predict which pattern shows the direction of smearing more clearly near the edge of the ray coverage.
Switch the true model between Checkerboard and Spike, and drag the green probe cell to an edge cell in each case.
Explain the difference with point (iii) of Section 2.5 and with Figures 2 and 3 of the paper.

Invert ray delays and probe the resolution

The Labs use flat, isotropic, elastic or acoustic rock with made-up but realistic speeds; none of them reproduces a real survey. Each Lab states its model and limits in its notes and can export the data behind its readouts.

Try explaining it yourself

1. A paper shows a checkerboard test with 3 per cent Gaussian noise and calls the recovered region “well resolved”. Which two effects does the figure mix, and what would separate them?

Read an explained answer

By Eq. 16 the output is the resolution matrix times the checkerboard plus the inverse applied to the noise. A noise-free test isolates the first; a noise-only inversion isolates the second.

2. Why can a perfectly recovered checkerboard still hide unresolved large-scale structure?

Read an explained answer

Recovery depends on which part of the model lies in the row space of the matrix. In the case Rawlinson and Spakman cite from Lévêque and co-authors (1993), a fine checkerboard was fully recovered while a coarse one lay almost entirely in the null space with the same rays. Tests at several scale lengths are needed.

3. The full resolution matrix took about 12 minutes and the spike test about 8 seconds. How large is the ratio, and why does the spike test still give nearly the same information in Figure 5?

Read an explained answer

About 90 times. Well-separated spikes each pick out one column of the resolution matrix with little overlap, so one inversion returns many columns at once.

\frac{12\cdot 60\ \mathrm s}{8\ \mathrm s}=90

Sources to read alongside

MIT OpenCourseWare, 12.201 Essentials of Geophysics (Fall 2004), lecture notes, Chapter 4: Seismology
Section 4.6 separates the elastic wave equation into P and S waves and gives their speeds (Eqs. 4.34 and 4.37). Sections 4.14 to 4.16 derive Snell’s law from Fermat’s principle (Eqs. 4.91 to 4.96) and show that the ray parameter is the slope of the travel-time curve (Eq. 4.98). Section 4.19 notes that in travel-time tomography the matrix entries are the ray lengths in each block.
N. Rawlinson and W. Spakman, On the use of sensitivity tests in seismic tomography, Geophysical Journal International 205, 1221–1243 (2016)
Section 2 writes the damped least-squares solution and its resolution matrix (Eqs. 7 to 14) and shows that a noise-free synthetic test samples the same resolution matrix (Eq. 16). Section 2.5 lists eight consequences. Section 3 compares checkerboard and spike tests on a 3844-unknown surface-wave example and ends with six recommendations.