Quick Start

Quick Start#

This page walks through a minimal end-to-end workflow: generate a geological model, compute elastic properties, simulate seismic data, and run a deterministic inversion.

1. Generate a Geological Model#

import torch
from geobrain.geomodel import Simulator, SimulationConfig

config = SimulationConfig(
    shape=(100, 200),
    lh=15, lv=5,
    mean=0.25, std=0.05,
    seed=2025,
)
sim = Simulator.create("fft_ma")
porosity = sim.simulate(config)
print(f"Porosity: shape={porosity.shape}, range=[{porosity.min():.3f}, {porosity.max():.3f}]")
Porosity field generated by FFT-MA

Fig. 7 Porosity field generated by FFT-MA.#

2. Rock Physics Forward Model#

Convert porosity to elastic velocities using the rock physics module:

from geobrain.physics.rock import SoftSand, Gassmann, DensityModel, v_from_moduli

soft_sand = SoftSand()
gassmann = Gassmann()

K_m, G_m = 36.6, 44.0        # Mineral moduli (GPa)
K_fl, rho_fl = 3.06, 1.08    # Fluid properties
rho_m = 2.65                  # Mineral density (g/cm3)

pressure = torch.ones_like(porosity) * 20
k_dry, g_dry = soft_sand(K_m, G_m, porosity, 0.4, 7, pressure)
k_sat, g_sat = gassmann(k_dry, g_dry, K_m, K_fl, porosity)
rho = DensityModel()(porosity, rho_m, rho_fl)
vp, vs = v_from_moduli(k_sat, g_sat, rho)
Rock physics: porosity to elastic velocities

Fig. 8 Rock physics: porosity to elastic velocities.#

3. Seismic Forward Modeling#

Generate a synthetic seismic section using AVO reflectivity and wavelet convolution:

from geobrain.physics.wave import (
    RickerWavelet, compute_reflectivity, create_conv_matrix,
)

# Reflectivity from layer contrasts
vp1, vp2 = vp[:-1], vp[1:]
vs1, vs2 = vs[:-1], vs[1:]
rho1, rho2 = rho[:-1], rho[1:]

refl = compute_reflectivity(vp1, vs1, rho1, vp2, vs2, rho2,
                            theta=[0], method='shuey')

# Wavelet convolution
ricker = RickerWavelet()
wavelet, _ = ricker(f0=25.0, dt=0.001)
W = create_conv_matrix(wavelet, refl.shape[-1])
nsW = len(wavelet) // 2
W = W[nsW:-nsW].T
seismic = refl @ W.T
Synthetic seismic angle gathers

Fig. 9 Synthetic seismic angle gathers.#

4. Deterministic Inversion#

Recover the velocity model from the synthetic seismic data:

from geobrain import InverseProblem

problem = InverseProblem(
    forward_fn=my_forward,      # Your forward operator
    observed=seismic,
    noise_std=0.01,
)

inverter = problem.create_inverter(initial_model=vp_init)
result = inverter.run(problem.observed, max_epochs=200, lr=0.01)
print(f"Final loss: {result.loss_history[-1]:.6f}")

5. Bayesian Uncertainty Quantification#

Quantify uncertainty using SVGD:

from geobrain import SVGD

posterior = problem.as_posterior(log_prior=my_prior)
svgd = SVGD(target=posterior, lr=0.005)
result = svgd.run(n_samples=50, n_steps=500)
print(f"Ensemble mean shape: {result.samples.mean(dim=0).shape}")
print(f"Ensemble std shape:  {result.samples.std(dim=0).shape}")
Bayesian inversion posterior comparison

Fig. 10 Bayesian inversion posterior comparison.#

Next Steps#

Contents