Back to Article
Reverse Time Migrations (RTM) and LS-RTM
Download Notebook

Exercise6: From velocity model to RTM and LS-RTM in JUDI

This notebook is a start-from-scratch imaging exercise for a seismic imaging class.
Students will go through the basic workflow:

  1. define or load a true velocity model,
  2. build a smooth migration/background model,
  3. generate synthetic seismic data,
  4. form an RTM image,
  5. run a few iterations of least-squares RTM (LS-RTM).

The goal is to distinguish clearly between:

  • data generation,
  • adjoint imaging,
  • inverse imaging.

The notebook is inspired by the JUDI inversion notes and the LS-RTM example script, but is rewritten as a guided class exercise that starts from a velocity model rather than from pre-generated data.

Learning goals

After this exercise, you should be able to explain:

  • what role the background velocity model plays in imaging,
  • how RTM relates to the adjoint of the Born operator,
  • why LS-RTM is an inverse problem and not just a single migration step,
  • why synthetic data generation and inversion are usually based on different models: a more detailed true model and a smoother migration model.

1. Background: imaging vs inversion

Let (F(m)) denote the nonlinear forward modeling operator, where (m) is squared slowness.

Synthetic data generation

For a true model (m_{}), synthetic data are generated by solving the wave equation: \[ d_{\text{true}} = F(m_{\text{true}})\,q, \] where (q) is the source wavelet.

Linearized modeling

To do migration and LS-RTM, we choose a smooth background model (m_0) and linearize around it: \[ d_{\text{Born}} \approx J(m_0)\,\delta m, \] where (J(m_0)) is the Born/Jacobian operator and \[ \delta m = m_{\text{true}} - m_0. \]

RTM

RTM applies the adjoint of the linearized modeling operator: \[ \delta m_{\text{RTM}} = J(m_0)^\top d. \] This is fast and useful, but it does not solve a least-squares problem.

LS-RTM

LS-RTM solves the inverse problem \[ \min_{\delta m}\ \frac{1}{2}\|J(m_0)\delta m - d\|_2^2. \] So LS-RTM iteratively improves the image by repeatedly applying forward and adjoint operators.

In the JUDI docs, RTM is described as the adjoint of the linearized modeling operator, while LS-RTM is implemented as a least-squares problem with matrix-free operators.

2. Packages

Run the next cell once.
If needed, uncomment the Pkg.add(...) lines the first time you use this notebook.

In [1]:
using Pkg
# Pkg.add(["JUDI", "PyPlot", "IterativeSolvers", "LinearAlgebra", "Statistics", "Random"])

using JUDI
using PyPlot
using IterativeSolvers
using LinearAlgebra
using Statistics
using Random

3. Define or load a true velocity model

This cell builds a simple synthetic 2D velocity model so the notebook is self-contained.
For your class, you can replace vp_true with your own velocity model array.

Convention used here - velocity is in km/s - grid spacing is in meters - time is in milliseconds

The model below has: - a slow near-surface layer, - deeper faster layers, - two embedded anomalies.

That gives students something simple but nontrivial to image.

In [6]:
# Grid
shape   = (301, 201)              # (nx, nz)
spacing = (10.0f0, 10.0f0)        # meters
origin  = (0.0f0, 0.0f0)

# True velocity model (replace this with your own model if desired)
vp_true = 1.5f0 .* ones(Float32, shape)
vp_true[:, 70:end]  .= 2.0f0
vp_true[:, 130:end] .= 2.5f0

# Add a high-velocity lens
for ix in 90:150, iz in 95:135
    if ((ix-120)/24)^2 + ((iz-115)/16)^2 <= 1
        vp_true[ix, iz] = 2.9f0
    end
end

# Add a low-velocity pocket
for ix in 185:245, iz in 120:155
    if ((ix-215)/24)^2 + ((iz-138)/14)^2 <= 1
        vp_true[ix, iz] = 2.15f0
    end
end

# Quick look
figure(figsize=(7,4))
imshow(vp_true', cmap="jet", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
colorbar(label="Velocity (km/s)")
xlabel("x (m)")
ylabel("z (m)")
title("True velocity model")
tight_layout()

Question 1

Why is this called the true model in a synthetic experiment?
What information would be unavailable in a real field-data problem?

Write 1–2 sentences here.

4. Build a smooth migration/background model

For migration and linearized inversion, we usually do not use the fully detailed true model as the background model.
Instead, we smooth it to obtain a kinematically reasonable but less detailed model.

This is important because the Born operator is defined around a background model.

Below, we use a simple repeated 5-point averaging smoother so the notebook stays lightweight.

In [9]:

function smooth2d(v; nrepeat=40)
    out = copy(v)
    tmp = copy(v)
    nx, nz = size(v)

    for _ = 1:nrepeat
        tmp .= out
        for ix in 2:nx-1, iz in 2:nz-1
            tmp[ix, iz] = (
                out[ix, iz] +
                out[ix-1, iz] + out[ix+1, iz] +
                out[ix, iz-1] + out[ix, iz+1]
            ) / 5.0f0
        end
        out .= tmp
    end
    return out
end

vp_mig = smooth2d(vp_true; nrepeat=60)

figure(figsize=(12,4))
subplot(1,2,1)
imshow(vp_true', cmap="jet", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("True velocity")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

subplot(1,2,2)
imshow(vp_mig', cmap="jet", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("Smooth migration velocity")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

tight_layout()

We now convert velocity to squared slowness and define:

  • the true model used for data generation,
  • the migration model used for RTM and LS-RTM,
  • the reflectivity / perturbation (m = m_{} - m_0).
In [11]:
m_true = 1.0f0 ./ vp_true.^2
m_mig  = 1.0f0 ./ vp_mig.^2
dm     = m_true .- m_mig

model_true = Model(shape, spacing, origin, m_true; nb=80)
model_mig  = Model(shape, spacing, origin, m_mig;  nb=80)

figure(figsize=(7,4))
imshow(dm', cmap="gray", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
colorbar(label="δm")
xlabel("x (m)")
ylabel("z (m)")
title("Model perturbation: δm = m_true - m_mig")
tight_layout()

Question 2

Why do we smooth the velocity model before building the Born operator?
What would happen if the background model tried to contain every reflector exactly?

Write 2–3 sentences here.

5. Acquisition geometry and source wavelet

We next define a simple surface acquisition geometry with: - sources near the surface, - receivers near the surface, - a Ricker wavelet.

This part closely mirrors the style used in your earlier imaging notebook, where students manually build the geometry before generating seismic data.

In [14]:
# Source positions
nsrc = 15
xsrc = range(0.0f0, (shape[1]-1)*spacing[1], length=nsrc)
ysrc = 0.0f0 .* xsrc
zsrc = 20.0f0 .* ones(Float32, nsrc)

# JUDI expects arrays-of-arrays for separate sources
xsrc, ysrc, zsrc = convertToCell.([xsrc, ysrc, zsrc])

# Receiver positions
nrec = 201
xrec = range(0.0f0, (shape[1]-1)*spacing[1], length=nrec)
yrec = 0.0f0
zrec = 20.0f0 .* ones(Float32, nrec)

# Time axis
record_time   = 2500.0f0      # ms
sampling_rate = 4.0f0         # ms

src_geom = Geometry(xsrc, ysrc, zsrc; dt=sampling_rate, t=record_time)
rec_geom = Geometry(xrec, yrec, zrec; dt=sampling_rate, t=record_time, nsrc=nsrc)

# Wavelet
f0 = 0.015f0                  # 15 Hz = 0.015 kHz because time is in ms
wavelet = ricker_wavelet(record_time, sampling_rate, f0)
q = judiVector(src_geom, wavelet)

# Plot model + acquisition
figure(figsize=(8,4))
imshow(vp_true', cmap="jet", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
scatter([xs[1] for xs in xsrc], zsrc, c="white", s=18, label="Sources")
scatter(xrec, zrec .* 1.0f0, c="black", s=8, label="Receivers")
xlabel("x (m)")
ylabel("z (m)")
title("Acquisition geometry on true model")
legend(loc="upper right")
tight_layout()

figure(figsize=(6,2.5))
plot(wavelet)
title("Ricker wavelet")
xlabel("Time sample")
ylabel("Amplitude")
tight_layout()

6. Generate synthetic seismic data

We will generate two data sets:

  1. Full-wave synthetic data \[ d_{ ext{true}} = F(m_{ ext{true}})\,q \] This is the more realistic data set.

  2. Born / linearized data \[ d_{ ext{Born}} = J(m_0)\,\delta m \] This is the consistent data set for the RTM and LS-RTM derivations below.

In class, it is often useful to show both: - d_true is what the full physics produces, - d_born is what the linearized theory assumes.

In [16]:
Random.seed!(1234)

opt = Options(IC="isic", dt_comp=1.0f0)

# Nonlinear forward operator on the true model
F_true = judiModeling(model_true, q.geometry, rec_geom; options=opt)

# Linearized/Born operator around the migration model
F_mig = judiModeling(model_mig, q.geometry, rec_geom; options=opt)
J = judiJacobian(F_mig, q)

# Synthetic data
d_true = F_true * q
d_born = J * vec(dm)

println("Number of shots: ", nsrc)
println("Model size     : ", prod(shape))
In [17]:
data_extent = [xrec[1], xrec[end], 1f-3*record_time, 0]

figure(figsize=(12,6))
for i in 1:4
    subplot(2,4,i)
    imshow(d_true.data[i], cmap="gray", aspect="auto", extent=data_extent)
    title("Full data, shot $i")
    xlabel("x (m)")
    ylabel("t (s)")

    subplot(2,4,4+i)
    imshow(d_born.data[i], cmap="gray", aspect="auto", extent=data_extent)
    title("Born data, shot $i")
    xlabel("x (m)")
    ylabel("t (s)")
end
tight_layout()

Question 3

Compare one full-wave shot record and one Born shot record.

What features look similar, and what differences do you notice?
Why might the Born data be preferable for a first LS-RTM exercise?

Write 2–3 sentences here.

7. RTM: adjoint imaging

The simplest migration image is obtained by applying the adjoint of the Born operator: \[ \delta m_{\mathrm{RTM}} = J^Td. \]

Below, we compute: - an RTM image from Born data (img_rtm_born), which is the most consistent test, - an RTM image from full-wave data (img_rtm_true), which is useful for discussion.

In [20]:
img_rtm_born = reshape(adjoint(J) * d_born, shape)
img_rtm_true = reshape(adjoint(J) * d_true, shape)

# Mute the top few grid points to suppress strong acquisition imprint
img_rtm_born[:, 1:10] .= 0.0f0
img_rtm_true[:, 1:10] .= 0.0f0

figure(figsize=(12,4))
subplot(1,2,1)
imshow(img_rtm_born', cmap="gray", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("RTM image from Born data")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

subplot(1,2,2)
imshow(img_rtm_true', cmap="gray", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("RTM image from full data")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

tight_layout()

8. LS-RTM: least-squares inversion

We now solve the inverse problem \[ \min_x \frac{1}{2}\|J_p x - d_p\|_2^2, \] where we optionally include simple preconditioning/muting operators:

  • Mr: a model-side top mute,
  • Ml: a data-side mute.

This follows the same operator-based idea as the JUDI LS-RTM examples: build matrix-free operators and then solve the least-squares problem iteratively.

In [22]:
Mr = judiTopmute(model_mig; taperwidth=10)
Ml = judiDataMute(q.geometry, rec_geom)

Jp = Ml * J * Mr
dp = Ml * d_born

niter = 12
x_lsrtm = zeros(Float32, prod(shape))

lsqr!(x_lsrtm, Jp, dp; maxiter=niter)

img_lsrtm = reshape(Mr * x_lsrtm, shape)
img_lsrtm[:, 1:10] .= 0.0f0

figure(figsize=(12,4))
subplot(1,2,1)
imshow(img_rtm_born', cmap="gray", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("RTM (Born data)")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

subplot(1,2,2)
imshow(img_lsrtm', cmap="gray", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("LS-RTM after $niter LSQR iterations")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

tight_layout()

Question 4

Compare the RTM and LS-RTM images.

Comment on: 1. sharpness, 2. amplitude balancing, 3. illumination / artifact reduction.

Write 3 short bullet points or 2–3 sentences here.

9. Residual check

A useful sanity check is whether the LS-RTM solution reduces the data misfit relative to the zero image.

In [25]:

r0 = dp
r_ls = Jp * x_lsrtm - dp

println("||r0||₂      = ", norm(r0))
println("||r_ls||₂    = ", norm(r_ls))
println("Reduction factor ||r_ls|| / ||r0|| = ", norm(r_ls)/norm(r0))

Question 5

Why is lower residual not the same thing as perfect geology?

Give one reason related to physics/modeling assumptions and one reason related to inversion/regularization.

Write 2–3 sentences here.

10. Exploratory questions

Try one or two of the following:

  • change the number of LSQR iterations,
  • reduce or increase the number of sources,
  • replace the synthetic model with your own class velocity model,
  • compare migration with a better and a worse smooth background model,
  • try using d_true instead of d_born inside RTM and discuss the differences.

Optional discussion prompt

If the migration velocity is poor, which part is most affected first: - event positioning, - focusing, - amplitudes, - or all of them?

Explain briefly.

11. Realistic Geology (Bonus question)

In the above problem, the velocity model was made up to keep it simple for the illustration. Use the Compass velocity model (https://www.dropbox.com/scl/fi/2unm2keg1h2rwbp1h5f2i/BGCompass_tti_625m.jld2?rlkey=1cbudz40vb89r23eusynn2chk&st=1fknqf9k&dl=0), resample it to the same size as above problem (or any smaller size) and perform RTM and LS-RTM.

12. Take-away summary

  • Data generation uses a true model.
  • RTM is the adjoint image (J^ op d).
  • LS-RTM solves a least-squares inverse problem iteratively.
  • The quality of both RTM and LS-RTM depends strongly on the background velocity model.
  • Starting from a velocity model and building the workflow from scratch helps separate the roles of: physics, acquisition, imaging, and inversion.