Reflection Seismology: CMP Processing

Reflection Seismology: CMP Processing#

Colab note: This notebook is designed to run on Google Colab. The first code cell installs dependencies. Open In Colab

Learning Objectives:

  • Sort seismic data into Common Midpoint (CMP) gathers

  • Perform Normal Moveout (NMO) correction

  • Stack traces to improve signal-to-noise ratio

  • Estimate stacking velocities

Prerequisites: Basic signal processing

Reference: Shearer, Chapter 6 (Reflection Seismology)

Notebook Outline:

We will:

  1. Generate a synthetic CMP gather (one reflector).

Hide code cell source

 
# Install dependencies (for Google Colab or missing packages)
import sys

# Check if running in Colab
try:
    import google.colab
    IN_COLAB = True
    print("Running in Google Colab")
except:
    IN_COLAB = False
    print("Running in local environment")

# Install required packages if needed
required_packages = {
    'numpy': 'numpy',
    'matplotlib': 'matplotlib',
    'scipy': 'scipy',
    'obspy': 'obspy'
}

missing_packages = []
for package, pip_name in required_packages.items():
    try:
        __import__(package)
        print(f"✓ {package} is already installed")
    except ImportError:
        missing_packages.append(pip_name)
        print(f"✗ {package} not found")

if missing_packages:
    print(f"\nInstalling missing packages: {', '.join(missing_packages)}")
    import subprocess
    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q"] + missing_packages)
    print("✓ Installation complete!")
else:
    print("\n✓ All required packages are installed!")

import numpy as np
import matplotlib.pyplot as plt

# ---------- utilities ----------
def ricker(t, f0):
    """Ricker wavelet (Mexican hat) with peak frequency f0 (Hz)."""
    a = (np.pi*f0*t)**2
    return (1 - 2*a) * np.exp(-a)

def add_event(trace, t, t0, amp=1.0, f0=20.0):
    """Add a Ricker wavelet centered at t0 to a trace."""
    trace += amp * ricker(t - t0, f0)
    return trace

def nmo_time(t0, x, v):
    """Hyperbolic moveout: t(x) = sqrt(t0^2 + (x/v)^2)."""
    return np.sqrt(t0**2 + (x/v)**2)

def apply_nmo(data, t, offsets, v, t0min=0.0):
    """NMO correction: map each trace from t(x) -> t0 using assumed v.
    Simple interpolation; good enough for a toy notebook."""
    dt = t[1]-t[0]
    out = np.zeros_like(data)
    for i, x in enumerate(offsets):
        # For each zero-offset time t0, find moveout time tx and sample original
        t0 = t.copy()
        tx = nmo_time(t0, x, v)
        # avoid extrapolating beyond record length
        tx = np.clip(tx, t[0], t[-1])
        out[:, i] = np.interp(t0, tx, data[:, i])
    return out

def stack_traces(data):
    """Simple average stack across traces."""
    return data.mean(axis=1)

np.random.seed(7)

1) Synthetic CMP gather: one reflector, hyperbolic moveout#

We model a single horizontal reflector at two-way zero-offset time \(t_0\) in a constant-velocity overburden. For each offset \(x\) the reflection arrives at:

\[ t(x)=\sqrt{t_0^2 + \left(\frac{x}{v}\right)^2} \]

Prediction (write it down)#

If I double the velocity v, does the curvature increase or decrease? Why?

# time axis
dt = 0.002
tmax = 3.0
t = np.arange(0, tmax+dt, dt)

# offsets (m) and model parameters
offsets = np.linspace(0, 2000, 41)   # 0–2 km
v_true = 2000.0                      # m/s (toy)
t0_true = 1.2                        # s two-way zero-offset time
f0 = 25.0                            # Hz wavelet peak

# build CMP gather: [nt, nx]
data = np.zeros((t.size, offsets.size))
for i, x in enumerate(offsets):
    tx = nmo_time(t0_true, x, v_true)
    add_event(data[:, i], t, tx, amp=1.0, f0=f0)

# add random noise
noise_level = 0.25
data_noisy = data + noise_level*np.random.normal(size=data.shape)

# plot gather
plt.figure()
plt.imshow(data_noisy, aspect='auto', origin='upper',
           extent=[offsets[0]/1000, offsets[-1]/1000, t[-1], t[0]])
plt.xlabel('Offset x (km)')
plt.ylabel('Time (s)')
plt.title('Synthetic CMP gather (noisy): one reflector')
plt.colorbar(label='amplitude')
plt.show()

2) Velocity estimate from a \(t^2\)–\(x^2\) fit (hyperbola linearization)#

For a single constant-velocity layer above the reflector:

\[ t(x)^2 = t_0^2 + \frac{x^2}{v^2} \]

So if you pick arrival times \(t(x)\), a line fit in \((x^2, t^2)\) gives:

  • intercept: \(t_0^2\)

  • slope: \(1/v^2\)

Prediction#

Noise in picks: does it bias \(v\) high or low, or mainly increase uncertainty?

# crude 'picker': take time of max absolute amplitude within a window
# (in real life: do semblance / coherence / interactive picking)
win = (t > 0.8) & (t < 1.8)

t_picks = np.zeros(offsets.size)
for i in range(offsets.size):
    idx = np.argmax(np.abs(data_noisy[win, i]))
    t_picks[i] = t[win][idx]

# linear fit in x^2, t^2
x2 = offsets**2
t2 = t_picks**2

coef = np.polyfit(x2, t2, 1)   # t^2 = a*x^2 + b
a, b = coef[0], coef[1]
v_est = 1/np.sqrt(a)
t0_est = np.sqrt(max(b, 0))

print(f"True v = {v_true:.1f} m/s  |  Estimated v = {v_est:.1f} m/s")
print(f"True t0 = {t0_true:.3f} s  |  Estimated t0 = {t0_est:.3f} s")

# plot t(x) picks and fit
plt.figure()
plt.plot(offsets/1000, t_picks, 'o', label='picks (toy)')
plt.plot(offsets/1000, nmo_time(t0_est, offsets, v_est), '-', label='fit hyperbola')
plt.plot(offsets/1000, nmo_time(t0_true, offsets, v_true), '--', label='true')
plt.xlabel('Offset x (km)')
plt.ylabel('Time (s)')
plt.legend()
plt.title('Moveout picks and hyperbola fit')
plt.gca().invert_yaxis()
plt.show()

# show linearization
plt.figure()
plt.plot(x2/1e6, t2, 'o', label='(x^2, t^2)')
plt.plot(x2/1e6, a*x2 + b, '-', label='line fit')
plt.xlabel('x^2 (km^2)')
plt.ylabel('t^2 (s^2)')
plt.legend()
plt.title('Linearized moveout: t^2 = t0^2 + x^2/v^2')
plt.show()

3) NMO correction + stack#

Apply NMO using: $\( t_0 = \sqrt{t(x)^2 - (x/v)^2} \)$ (implemented by interpolation in the notebook).

Then stack across offsets:

  • Signal that aligns at \(t_0\) adds coherently.

  • Random noise adds incoherently (SNR improves ~ \(\sqrt{N}\) for uncorrelated noise).

Prediction#

Try a deliberately wrong velocity:

  • If \(v\) is too fast, does the corrected gather “smile” upward or downward?

# NMO with estimated velocity
data_nmo = apply_nmo(data_noisy, t, offsets, v_est)
stk = stack_traces(data_nmo)

# Compare with wrong velocities
v_slow = 0.8*v_true
v_fast = 1.2*v_true
data_nmo_slow = apply_nmo(data_noisy, t, offsets, v_slow)
data_nmo_fast = apply_nmo(data_noisy, t, offsets, v_fast)

# plots
def plot_gather(d, title):
    plt.figure()
    plt.imshow(d, aspect='auto', origin='upper',
               extent=[offsets[0]/1000, offsets[-1]/1000, t[-1], t[0]])
    plt.xlabel('Offset x (km)')
    plt.ylabel('Time (s)')
    plt.title(title)
    plt.colorbar(label='amplitude')
    plt.show()

plot_gather(data_nmo, f'NMO corrected gather (v_est={v_est:.0f} m/s)')
plot_gather(data_nmo_slow, f'NMO corrected gather (too slow: v={v_slow:.0f} m/s)')
plot_gather(data_nmo_fast, f'NMO corrected gather (too fast: v={v_fast:.0f} m/s)')

# stack vs a single trace
plt.figure()
plt.plot(t, data_noisy[:, offsets.size//2], label='single mid-offset trace', alpha=0.8)
plt.plot(t, stk, label='NMO stack', linewidth=2)
plt.xlabel('Time (s)')
plt.ylabel('amplitude')
plt.title('Stack boosts coherent reflection')
plt.legend()
plt.gca().invert_xaxis()  # optional: many reflection plots use down-going time; keep normal here
plt.show()

4) Diffraction hyperbola and a toy “migration” (collapse to the apex)#

A point diffractor (or an edge/corner) in a zero-offset section produces a diffraction hyperbola:

\[ t(x)=\frac{2}{v}\sqrt{h^2 + x^2} \]

Idea of migration (diffraction summation):

  • Assume each image point is the apex of a hyperbola.

  • Sum along that hyperbola.

  • Diffractions collapse back to their true subsurface location if velocity is correct.

We’ll implement the concept in 2D (x–t) with a single diffractor.

Prediction#

If migration velocity is too slow, do diffractions collapse too much (overmigrate) or too little (undermigrate)?

# Build a toy zero-offset section: one point diffractor
x = np.linspace(-2000, 2000, 81)  # m
t = np.arange(0, 3.0+dt, dt)

h = 1000.0          # diffractor depth (m)
v = 2000.0          # m/s
f0 = 25.0

section = np.zeros((t.size, x.size))
for ix, xx in enumerate(x):
    tx = (2/v)*np.sqrt(h**2 + xx**2)
    add_event(section[:, ix], t, tx, amp=1.0, f0=f0)

section += 0.15*np.random.normal(size=section.shape)

plt.figure()
plt.imshow(section, aspect='auto', origin='upper',
           extent=[x[0]/1000, x[-1]/1000, t[-1], t[0]])
plt.xlabel('Distance x (km)')
plt.ylabel('Time (s)')
plt.title('Zero-offset section with a point diffractor (diffraction hyperbola)')
plt.colorbar(label='amplitude')
plt.show()

# --- toy migration: for each image point (x0, t0) we sum along its predicted hyperbola ---
def migrate_diffraction_summation(section, t, x, v_mig, t0_grid):
    """Return a migrated image on (t0_grid, x). For each (x0, t0),
    sum data along the corresponding diffraction hyperbola in the input section."""
    out = np.zeros((t0_grid.size, x.size))
    for ix0, x0 in enumerate(x):
        # for each output time t0, sample along hyperbola: t(x) = sqrt(t0^2 + (2*(x-x0)/v)^2)? 
        # For exploding reflector / zero-offset, a convenient mapping is:
        # t_in(x) = sqrt(t0^2 + (2*(x-x0)/v)^2). (captures curvature)
        for it0, t0 in enumerate(t0_grid):
            tin = np.sqrt(t0**2 + (2*(x - x0)/v_mig)**2)
            tin = np.clip(tin, t[0], t[-1])
            vals = np.array([np.interp(tin[j], t, section[:, j]) for j in range(x.size)])
            out[it0, ix0] = vals.mean()
    return out

t0_grid = np.arange(0.3, 2.5, 0.01)

mig_true = migrate_diffraction_summation(section, t, x, v_mig=v, t0_grid=t0_grid)
mig_slow = migrate_diffraction_summation(section, t, x, v_mig=0.8*v, t0_grid=t0_grid)
mig_fast = migrate_diffraction_summation(section, t, x, v_mig=1.2*v, t0_grid=t0_grid)

def plot_mig(img, title):
    plt.figure()
    plt.imshow(img, aspect='auto', origin='upper',
               extent=[x[0]/1000, x[-1]/1000, t0_grid[-1], t0_grid[0]])
    plt.xlabel('Distance x (km)')
    plt.ylabel('t0 (s)')
    plt.title(title)
    plt.colorbar(label='amplitude')
    plt.show()

plot_mig(mig_true, 'Toy migrated image (v correct)')
plot_mig(mig_slow, 'Toy migrated image (v too slow)')
plot_mig(mig_fast, 'Toy migrated image (v too fast)')

Wrap-up prompts (discussion)#

  1. What did NMO need? What information about the Earth does it assume?

  2. Velocity errors: In which steps do they hurt you most: NMO, stack, migration?

  3. Noise types: Which arrivals don’t stack coherently (ground roll, direct wave, multiples)?

  4. Bridge to “real”: What changes in a layered velocity gradient v(z)? What breaks the simple hyperbola?

Extension (homework / advanced):

  • Add a second reflector and a multiple; see what stacking does.

  • Replace constant v with v(z) and ray-trace a moveout curve.

  • Add random time shifts and solve for shifts that maximize stack power.

Contents