Exercise 7: From velocity model to subsurface offset gathers and extended LS-RTM in JUDI/ImageGather

This notebook is a start-from-scratch imaging and inversion exercise for a seismic imaging class.

You will go through this workflow:

define or load a true velocity model,
build a smooth background / migration model,
generate synthetic seismic data,
compute subsurface common offset gathers (SCOGs / CIGs),
run extended LS-RTM (eLSRTM) with an extended Born operator.

This replaces the RTM / LS-RTM exercise with an extended imaging point of view: instead of producing only one migrated image, we build an image volume that also depends on subsurface offset.

Learning goals

After this exercise, you should be able to explain:

why we distinguish between a true model and a background model,
how subsurface offset gathers are computed from the adjoint of an extended Born operator,
why the zero-offset image is only one slice of a larger extended image volume,
how extended LS-RTM differs from simply applying the adjoint once,
how event focusing / defocusing in the offset dimension is connected to background-model quality.

1. Background: imaging vs extended imaging vs inversion

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

Synthetic data generation

For a true model (m_{}), synthetic data are generated as [ d_{} = F(m_{}),q. ]

Background-model prediction

For a smooth background model (m_0), we can also generate [ d_0 = F(m_0),q. ]

The data residual [ d = d_{} - d_0 ] contains the scattered information that the background model fails to explain.

Extended linearized modeling

In extended imaging, we replace the standard Born operator by an extended Jacobian [ J_{}(m_0), ] whose adjoint produces a subsurface offset gather volume [ I(x,z,h) = J_{}(m_0)^d. ]

Here (h) is the subsurface horizontal offset.

Extended LS-RTM

Instead of using only the adjoint image, we solve [ x |J{} x - d|_2^2, ] where (x) is now an extended image volume. This is often called extended LS-RTM.

2. Packages

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

In [5]:

using Pkg
# Pkg.add(["PyPlot", "IterativeSolvers", "LinearAlgebra", "Statistics", "Random"])
# Pkg.add(name="JUDI", version="3.3.8")
# Pkg.add(name="ImageGather", version="0.2.6")

using JUDI
using ImageGather
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 model so the notebook is self-contained.
For 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

In [7]:

# Grid
shape   = (301, 201)              # (nx, nz)
n = shape
spacing = (10.0f0, 10.0f0)        # meters
origin  = (0.0f0, 0.0f0)

# True velocity model (replace 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

# 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

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

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 unknown in a real field-data problem?

Write 1–2 sentences here.

4. Build a smooth background model

For migration / linearized inversion we usually work with a smooth background model. This model is meant to capture the main kinematics without containing every reflector explicitly.

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

In [10]:

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_bg = 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_bg', cmap="jet", aspect="auto",
       extent=[0, (shape[1]-1)*spacing[1], (shape[2]-1)*spacing[2], 0])
title("Smooth background velocity")
xlabel("x (m)")
ylabel("z (m)")
colorbar()

tight_layout()

We now convert velocity to squared slowness and define the JUDI models.

In [12]:

m_true = 1.0f0 ./ vp_true.^2
m_bg   = 1.0f0 ./ vp_bg.^2
dm     = m_true .- m_bg

model_true = Model(shape, spacing, origin, m_true; nb=80)
model_bg   = Model(shape, spacing, origin, m_bg;   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_bg")
tight_layout()

Question 2

Why do we smooth the velocity model before constructing the imaging / inversion operator?
How is this related to kinematics versus reflectivity?

Write 2–3 sentences here.

5. Acquisition geometry and source wavelet

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

This follows the same spirit as the earlier acoustic imaging exercise: students explicitly build the grid, geometry, and wavelet before imaging.

In [15]:

# 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)

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

figure(figsize=(6,3))
plot((0:length(wavelet)-1) .* sampling_rate .* 1f-3, wavelet)
xlabel("Time (s)")
ylabel("Amplitude")
title("Source wavelet")
grid(true)
tight_layout()

6. Generate synthetic data

We compute:

observed data from the true model,
predicted data from the background model,
the data residual used in extended imaging: [ d = d_{} - d_0. ]

Optionally, we can add band-limited random noise.

In [17]:

Random.seed!(1234)

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

F_true = judiModeling(model_true, q.geometry, rec_geom; options=opt)
F_bg   = judiModeling(model_bg,   q.geometry, rec_geom; options=opt)

@time d_obs = F_true * q
@time d0    = F_bg * q

# Residual data
d_res = d_obs - d0

# Optional: add band-limited noise to the observed data
add_noise = false
snr = 15f0

if add_noise
    noise_ = deepcopy(d_obs)
    for l = 1:nsrc
        noise_.data[l] = randn(Float32, size(d_obs.data[l]))
        noise_.data[l] = real.(ifft(fft(noise_.data[l]) .* fft(q.data[1])))
    end
    noise_ = noise_ / norm(noise_) * norm(d_obs) * 10f0^(-snr/20f0)
    d_obs = d_obs + noise_
    d_res = d_obs - d0
end

println("Number of shots: ", nsrc)
println("Residual norm  : ", norm(d_res))

In [18]:

data_extent = [xrec[1], xrec[end], 1f-3*record_time, 0]

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

    subplot(3,3,3+i)
    imshow(d0.data[i], cmap="gray", aspect="auto", extent=data_extent)
    title("Background, shot $i")
    xlabel("x (m)")
    ylabel("t (s)")

    subplot(3,3,6+i)
    imshow(d_res.data[i], cmap="gray", aspect="auto", extent=data_extent)
    title("Residual, shot $i")
    xlabel("x (m)")
    ylabel("t (s)")
end
tight_layout()

Question 3

Compare d_obs, d0, and d_res.

Which events are already explained by the background model, and which events remain in the residual?
Why is the residual the natural input for linearized imaging?

Write 2–3 sentences here.

7. Subsurface common offset gathers (extended adjoint imaging)

We now build an extended Jacobian and apply its adjoint to the residual data.

The result is a 3D volume: [ I(x,z,h) = J_{}^ op d, ] where (h) is the subsurface horizontal offset.

The offsetrange below defines the offsets we want in the extended image.

In [21]:

# Offsets for the extended image volume (in meters)
n_offsets = 41
offset_start = -200
offset_end = 200
offsetrange = range(offset_start, stop=offset_end, length=n_offsets)
nh = length(offsetrange)

J_ext = judiExtendedJacobian(F_bg(m_bg), q, offsetrange)

@time cig = adjoint(J_ext) * d_res

# Mute shallow cells to reduce source/receiver imprint
cig[:, :, 1:20] .= 0.0f0

println("Extended image size: ", size(cig))
println("Offsets (m): ", offsetrange)

In [66]:

# Butterfly plot
function plot_cig(cig)
    d = (10.0f0, 10.0f0)
    # font size and label size
    fs = 10
    ls = 10
    cig_fs = 32
    
    # start plotting
    y = reshape(permutedims(cig, [2, 3, 1]), n[1], n[2], n_offsets, 1);
    PyPlot.rc("figure", titlesize=cig_fs)
    PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=cig_fs); PyPlot.rc("ytick", labelsize=cig_fs)
    PyPlot.rc("axes", labelsize=cig_fs)     # Default fontsize for x and y labels
    PyPlot.rc("axes", titlesize=cig_fs)     # Default fontsize for titles
    ### X, Z position in km
    
    xpos = 1.2f3
    zpos = 1.3f3
    # xpos = 2.4f3
    # zpos = 1.0f3
    xgrid = Int(round(xpos / d[1]))
    zgrid = Int(round(zpos / d[2]))
    # Create a figure and a 2x2 grid of subplots
    fig, axs = subplots(2, 2, figsize=(20,12), gridspec_kw = Dict("width_ratios" => [3, 1], "height_ratios" => [1, 3]))
    # Adjust the spacing between the plots
    subplots_adjust(hspace=0.0, wspace=0.0)
    
    # vmin1, vmax1 = (-1, 1) .* 23716
    # vmin2, vmax2 = (-1, 1) .* 22608
    # vmin3, vmax3 = (-1, 1) .* 42511
    vmin1, vmax1 = (-1, 1) .* quantile(abs.(vec(y[:,zgrid,:,1])), 0.98)
    vmin2, vmax2 = (-1, 1) .* quantile(abs.(vec(y[:,:,div(n_offsets,2)+1,1])), 0.88)
    vmin3, vmax3 = (-1, 1) .* quantile(abs.(vec(y[xgrid,:,:,1])), 0.94)
    sca(axs[1, 1])
    # Top left subplot
    axs[1, 1].imshow(y[:,zgrid,:,1]', aspect="auto", cmap="gray", interpolation="none",vmin=vmin1, vmax=vmax1, extent=(0f0, (n[1]-1)*d[1], offset_start, offset_end))
    # axs[1, 1].imshow(y[:,zgrid,:,1]', aspect="auto", cmap="gray", interpolation="none",vmin=-0.2, vmax=0.2, extent=(0f0, (n[1]-1)*d[1], offset_start, offset_end))
    # axs[1, 1].set_ylabel("Offset [m]", fontsize=cig_fs)
    axs[1, 1].set_xticklabels([])
    axs[1, 1].set_xlabel("")
    hlines(y=0, colors=:b, xmin=0, xmax=(n[1]-1)*d[1], linewidth=3)
    vlines(x=xpos, colors=:b, ymin=offset_start, ymax=offset_end, linewidth=3)
    # Bottom left subplot
    sca(axs[2, 1])
    axs[2, 1].imshow(y[:,:,div(n_offsets,2)+1,1]', aspect="auto", cmap="gray", interpolation="none",vmin=vmin2, vmax=vmax2, extent=(0f0, (n[1]-1)*d[1], (n[2]-1)*d[2], 0f0))
    # axs[2, 1].imshow(y[:,:,div(n_offsets,2)+1,1]', aspect="auto", cmap="gray", interpolation="none",vmin=-0.2, vmax=0.2,
    # extent=(0f0, (n[1]-1)*d[1], (n[2]-1)*d[2], 0f0))
    axs[2, 1].set_xlabel("X [m]", fontsize=cig_fs)
    axs[2, 1].set_ylabel("Z [m]", fontsize=cig_fs)
    axs[2, 1].set_xticks([0, 500, 1000, 1500, 2000, 2500])
    axs[2, 1].set_xticklabels(["0", "500", "1000", "1500", "2000", "2500"])
    axs[2, 1].set_yticks([500, 1000, 1500])
    axs[2, 1].set_yticklabels(["500", "1000", "1500"])
    # axs[2, 2].get_shared_x_axes().join(axs[1, 1], axs[2, 1])
    vlines(x=xpos, colors=:b, ymin=0, ymax=(n[2]-1)*d[2], linewidth=3)
    hlines(y=zpos, colors=:b, xmin=0, xmax=(n[1]-1)*d[1], linewidth=3)
    # Top right subplot
    axs[1, 2].set_visible(false)
    # Bottom right subplot
    sca(axs[2, 2])
    axs[2, 2].imshow(y[xgrid,:,:,1], aspect="auto", cmap="gray", interpolation="none",vmin=vmin3, vmax=vmax3,
    extent=(offset_start, offset_end, (n[2]-1)*d[2], 0f0))
    # axs[2, 2].imshow(y[xgrid,:,:,1], aspect="auto", cmap="gray", interpolation="none",vmin=-0.2, vmax=0.2,
    # extent=(offset_start, offset_end, (n[2]-1)*d[2], 0f0))
    axs[2, 2].set_xlabel("Offset [m]", fontsize=cig_fs)
    # Share y-axis with bottom left
    # axs[2, 2].get_shared_y_axes().join(axs[2, 2], axs[2, 1])
    axs[2, 2].set_yticklabels([])
    axs[2, 2].set_ylabel("")
    vlines(x=0, colors=:b, ymin=0, ymax=(n[2]-1)*d[2], linewidth=3)
    hlines(y=zpos, colors=:b, xmin=offset_end, xmax=offset_start, linewidth=3)
    # Remove the space between subplots and hide the spines
    for ax in reshape(axs, :)
        for spine in ["top", "right", "bottom", "left"]
            ax.spines[spine].set_visible(false)
        end
    end
    
    PyPlot.show()
end

plot_cig (generic function with 2 methods)

In [22]:

plot_cig(cig)

9. Extended LS-RTM

We now solve the least-squares problem [ x rac12 |J{} x - d|_2^2, ] where (x) is an extended image volume.

This is the extended analogue of LS-RTM: instead of solving for one image (x(x,z)), we solve for (x(x,z,h)).

In [24]:

niter = 10

# extended image in the same shape returned by the adjoint
x = zeros(Float32, size(cig))
fval = zeros(Float32, niter)

for k = 1:niter
    println("Iteration ", k)

    # forward
    r = J_ext * x - d_res

    # gradient
    g = J_ext' * r

    # basic step length
    α = norm(r)^2 / norm(g)^2

    x .-= α .* g
    fval[k] = 0.5f0 * norm(r)^2
end

img_ext_ls = copy(x)          # (nh, nx, nz) in your case
img_ext_ls[:, :, 1:20] .= 0f0

In [25]:

plot_cig(img_ext_ls)

Question 5

Compare the adjoint extended image and the extended LS-RTM result.

Comment on: 1. sharpness, 2. amplitude balancing, 3. artifact suppression.

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

10. Residual check

A useful sanity check is whether the eLSRTM solution reduces the data misfit compared with the zero extended image.

In [28]:

r0 = d_res
r_ls = J_ext * x - d_res

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

Question 6

Why does a lower least-squares residual not automatically mean that the recovered image is geologically perfect?

Give one reason related to physics / modeling and one reason related to inversion / nonuniqueness.

Write 2–3 sentences here.

11. Optional exploration

Try one or two of the following:

change the offsetrange,
increase or decrease the number of LS iterations,
add noise by setting add_noise = true,
compare a better and a worse background model,
choose several ix_gather values and compare local gathers,
replace the synthetic model with your own class velocity model.

Optional discussion prompt

If the background velocity is poor, where do you expect the error to appear first: - in the data residual, - in the focusing across subsurface offset, - in the stacked image?

12. Take-away summary

Observed data come from the true model.
Residual data are formed by subtracting the background-model prediction.
The adjoint of the extended Jacobian produces a subsurface offset gather volume.
The zero-offset slice is only one part of that extended image.
Extended LS-RTM solves a least-squares inverse problem for the full extended image volume.
The quality of the gathers and the inversion depends strongly on the background velocity model.

Common Image Gathers and extended LS-RTM

Learning goals

1. Background: imaging vs extended imaging vs inversion

Synthetic data generation

Background-model prediction

Extended linearized modeling

Extended LS-RTM

2. Packages

3. Define or load a true velocity model

Question 1

4. Build a smooth background model

Question 2

5. Acquisition geometry and source wavelet

6. Generate synthetic data

Question 3

7. Subsurface common offset gathers (extended adjoint imaging)

9. Extended LS-RTM

Question 5

10. Residual check

Question 6

11. Optional exploration

Optional discussion prompt

12. Take-away summary