Seismic Interpretation Workflow

End-to-end pipeline for seismic data analysis, from loading SEG-Y files through signal processing, rock physics modelling, and 3D visualization.

Skill Chain

segyio          obspy           bruges          disba           pyvista
[SEG-Y I/O] --> [Signal Proc] --> [Rock Physics] --> [Dispersion] --> [3D Viz]
  |               |                |                 |                |
  Load traces     Filter/FFT       AVO modelling     Surface waves    Volume render
  Read headers    Instrument resp   Fluid sub         Phase velocity   Slice display
  3D geometry     Spectral anal    Synthetics        Inversion        Horizon pick

Decision Points

Step-by-Step Orchestration

Stage 1: Data Loading (segyio)

import segyio
import numpy as np

# Load and inspect SEG-Y
with segyio.open('survey.sgy', 'r', iline=189, xline=193) as f:
    n_traces = f.tracecount
    n_samples = len(f.samples)
    dt = f.samples[1] - f.samples[0]  # ms
    cube = segyio.tools.cube(f)        # (ilines, xlines, samples)

    # Extract header info for coordinates
    cdp_x = f.attributes(segyio.TraceField.CDP_X)[:]
    cdp_y = f.attributes(segyio.TraceField.CDP_Y)[:]

Stage 2: Signal Processing (obspy)

from obspy import Trace, Stream
from obspy.signal.filter import bandpass

# Convert numpy trace to obspy for processing
tr = Trace(data=cube[50, 100, :].astype(np.float64))
tr.stats.sampling_rate = 1000.0 / dt  # Hz
tr.stats.delta = dt / 1000.0          # seconds

# Bandpass filter
tr.filter('bandpass', freqmin=5, freqmax=80, corners=4, zerophase=True)

# Spectral analysis
from obspy.signal.tf_misfit import cwt
scalogram = cwt(tr.data, dt=tr.stats.delta, w0=6, fmin=1, fmax=100, nf=50)

Stage 3: Rock Physics (bruges)

from bruges.reflection import zoeppritz, shuey
from bruges.filters import ricker
from bruges.rockphysics import gassmann
from scipy.signal import convolve

# AVO modelling
theta = np.arange(0, 45, 1)
Rpp = zoeppritz(vp1, vs1, rho1, vp2, vs2, rho2, theta)

# Fluid substitution
vp_new, vs_new, rho_new = gassmann(
    vp_sat, vs_sat, rho_sat,
    k_mineral, rho_mineral,
    k_brine, rho_brine,
    k_oil, rho_oil,
    porosity
)

# Synthetic seismogram
impedance = vp * rho
rc = np.diff(impedance) / (impedance[:-1] + impedance[1:])
_, wavelet = ricker(duration=0.128, dt=0.001, f=25)
synthetic = convolve(rc, wavelet, mode='same')

Stage 4: Surface Wave Analysis (disba, optional)

from disba import PhaseDispersion

# Define velocity model (thickness, Vp, Vs, density)
velocity_model = disba.PhaseDispersion(
    thickness=[5, 10, 20, 0],       # m (0 = half-space)
    velocity_p=[300, 600, 1200, 2500],
    velocity_s=[150, 300, 600, 1200],
    density=[1.8, 1.9, 2.1, 2.4]
)
# Compute Rayleigh wave dispersion
pd = velocity_model(wave="rayleigh", mode=0, period=np.linspace(0.01, 1.0, 50))

Stage 5: Visualization (pyvista)

import pyvista as pv

# Create structured grid from seismic cube
grid = pv.ImageData(dimensions=cube.shape)
grid.point_data['amplitude'] = cube.flatten(order='F')

# Volume rendering
plotter = pv.Plotter()
plotter.add_volume(grid, scalars='amplitude', cmap='seismic', opacity='sigmoid')
plotter.show()

# Inline slice
inline_slice = grid.slice(normal='x', origin=grid.center)
plotter = pv.Plotter()
plotter.add_mesh(inline_slice, scalars='amplitude', cmap='seismic')
plotter.show()

Common Pipelines

Basic Seismic QC

- [ ] Load SEG-Y with `segyio.open()`, check trace count and geometry
- [ ] Inspect text and binary headers for acquisition parameters
- [ ] Read trace headers to verify coordinates and fold
- [ ] Compute amplitude statistics (min, max, RMS) per trace
- [ ] Display representative inlines/crosslines with matplotlib
- [ ] Check for dead traces (zero or constant amplitude)
- [ ] Verify sample rate and recording length

AVO Analysis Pipeline

- [ ] Load near and far angle stacks from SEG-Y (segyio)
- [ ] Extract well log Vp, Vs, density at target zone (lasio/welly)
- [ ] Compute AVO response with `bruges.reflection.zoeppritz()`
- [ ] Classify AVO anomaly (Class I-IV) from intercept and gradient
- [ ] Perform Gassmann fluid substitution for scenario modelling
- [ ] Generate synthetic seismogram and compare to seismic
- [ ] Create AVO crossplot (intercept vs gradient)

Surface Wave Analysis (MASW)

- [ ] Load shot gathers from SEG-Y (segyio)
- [ ] Apply f-k or phase-shift transform to extract dispersion image
- [ ] Pick fundamental mode dispersion curve
- [ ] Invert dispersion curve for 1D Vs profile using disba
- [ ] Repeat for multiple shot points
- [ ] Interpolate Vs profiles into 2D section (verde)
- [ ] Visualize with pyvista or matplotlib

Seismic-to-Well Tie

- [ ] Load seismic section (segyio) and well logs (lasio)
- [ ] Compute acoustic impedance from Vp and density logs
- [ ] Generate reflectivity series from impedance
- [ ] Create Ricker wavelet at dominant frequency (bruges)
- [ ] Convolve reflectivity with wavelet for synthetic
- [ ] Extract seismic trace at well location
- [ ] Cross-correlate synthetic with seismic to find time shift
- [ ] Apply bulk shift and stretch/squeeze for best tie

When to Use

Use the seismic interpretation workflow when:

• Processing seismic data from SEG-Y files through to interpretation
• Building AVO analysis pipelines that combine well logs and seismic
• Performing seismic QC, attribute extraction, or volume visualization
• Running surface wave analysis (MASW/SASW) workflows
• Creating synthetic seismograms for well ties

Use individual domain skills when:

• Only reading/writing SEG-Y files (use segyio alone)
• Only processing waveforms without seismic context (use obspy alone)
• Only computing rock physics equations (use bruges alone)

Common Issues