In []:

Introduction to the Python Seismic Imaging Toolbox (PySIT)

MSRI: Introduction to the Mathematics of Seismic Imaging

Russell J. Hewett

Department of Mathematics, Massachusetts Institute of Technology

Python Seismic Imaging Toolbox¶

Full Waveform Inversion (FWI) + Numerical Python

Reproducible FWI experiments in controlled environment

Modular construction for investigating FWI from multiple perspectives

Simple and flexible interface for rapid prototyping of new ideas

(Soon to be) open source with open development model

Python is write once, run everywhere
Construct experiments in Python scripting environment + version control

In []:

################ CARTOON (doesn't run...yet) ###############

# Setup a physical problem and acquisition
M, M0, mesh, domain = marmousi2(scale='mini')

shots = equispaced_acquisition(mesh, domain, nshots=5, nrecs=50)

# Setup solver
solver = ConstantDensityAcousticWave(mesh, trange=(0.0, 3.0))

# Populate synthetic data
generate_seismic_data(shots, solver, base_model)

# Define objective function
objective = TemporalLeastSquares(solver)

# Define optimization algorithm 
invalg = LBFGS(objective)

# Run inversion for 5 steps
result = invalg(shots, initial_value, 5, line_search='backtrack')

# Visualize results
vis.plot(result.C, mesh)

Physical setup and acquisition

Optimization

\[\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Updownarrow\]

Objective Function

\[\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Updownarrow\]

(Modeling)

\[\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Updownarrow\]

Solver

Gradient Descent
Limited Memory BFGS (L-BFGS)
(Gauss-Newton via Krylov)

Backtracking line search

Store full history of optimization process for analysis

Objective Functions¶

Implements

Evaluating objective
Computing gradient
Applying Hessian to perturbation

Time domain least-squares: \(J(m) = \frac{1}{2} \int_0^T || d(t) - \mathbf{S}\mathcal{F}[m](t)||^2\)

Frequency domain least-squares: \(J(m) = \frac{1}{2} \sum_{\omega \in \Omega} || \hat d(\omega) - {\hat{\mathbf{S}\mathcal{F}[m]}}(\omega)||^2\)

Hybrid least-squares: \(J(m) = \frac{1}{2} \sum_{\omega \in \Omega} || \hat d(\omega) - {\hat{\mathbf{S}\mathcal{F}[m]}}(\omega)||^2\)

Modeling

Implementation of \(\mathcal{F}, F,\) and \(F^*\)

Wave: Time-domain Constant Density Acoustics
- Scalar (leap frog) and ODE implementations
- Pure numerical Python and C++ accelerated

Helmholtz: Frequency-domain Constant Density Acoustics
- Scalar implementations
- Utilizes scipy.sparse

Boundary Conditions
- Dirichlet
- PML: Perfectly Matched Layers (advanced absorbing boundary)
- (Neumann)

Generalized to 1D, 2D, and 3D
- Arbitrary order spatial FD

Simple Rapid Prototyping¶

In []:

################ CARTOON (doesn't run...yet) ###############

# Setup a physical problem and acquisition
M, M0, mesh, domain = marmousi2(scale='mini')

shots = equispaced_acquisition(mesh, domain, nshots=5, nrecs=50)

# Setup solver
solver = ConstantDensityAcousticWave(mesh, trange=(0.0, 3.0))

# Populate synthetic data
generate_seismic_data(shots, solver, base_model)

# Define objective function
objective = TemporalLeastSquares(solver)

# Define optimization algorithm 
invalg = LBFGS(objective)

# Run inversion for 5 steps
result = invalg(shots, initial_value, 5, line_search='backtrack')

# Visualize results
vis.plot(result.C, mesh)

Simple Rapid Prototyping¶

In []:

################ CARTOON (doesn't run...yet) ###############

# Setup a physical problem and acquisition
M, M0, mesh, domain = marmousi2(scale='mini')

shots = equispaced_acquisition(mesh, domain, nshots=5, nrecs=50)

# Setup solver
solver = ConstantDensityAcousticWaveLOD(mesh, trange=(0.0, 3.0))

# Populate synthetic data
generate_seismic_data(shots, solver, base_model)

# Define objective function
objective = TemporalLeastSquares(solver)

# Define optimization algorithm 
invalg = LBFGS(objective)

# Run inversion for 5 steps
result = invalg(shots, initial_value, 5, line_search='backtrack')

# Visualize results
vis.plot(result.C, mesh)

2D Horizontal Reflector Demo¶

Start an ipython interpreter: ipython --pylab

In [1]:

import numpy as np
import matplotlib.pyplot as plt 

import pysit # from pysit import *

from pysit.gallery import horizontal_reflector

Define the Physical Domain¶

In [2]:

# Define a boundary condition object
bcx = pysit.PML(0.1, 100)
# bcx = Dirichlet()
bcz = pysit.PML(0.1, 100)

# Configure the domain in PHYSICAL paramters
x_config = (0.1, 1.0, bcx, bcx)
z_config = (0.1, 0.8, bcz, bcz)

In [3]:

# Define the PHYSICAL domain based on the parameters
domain = pysit.RectangularDomain( x_config, z_config )

Define Computational Mesh¶

In [4]:

# Define the COMPUTATIONAL mesh based on the domain
mesh = pysit.CartesianMesh(domain, 90, 70)

Define and Visualize Demo Problem¶

In [5]:

# Create our demo model
C, C0 = horizontal_reflector(mesh, reflector_position=[0.45, 0.65], 
                             reflector_scaling = [1.0, 1.0], base_speed = 1.0,)

(90L, 70L) (6300, 1)

In [6]:

pysit.vis.plot(C0, mesh, axis=plt.subplot(1,2,1), clim=(C.min(),C.max())); 
plt.colorbar();

pysit.vis.plot(C, mesh, axis=plt.subplot(1,2,2)) 
plt.colorbar();

Setup Shots/Data Acquisition¶

In [7]:

# Extract some useful parameters
zmin = domain.z.lbound
zmax = domain.z.rbound
zdepth = zmin + (1./9.)*zmax

In [8]:

shots = pysit.equispaced_acquisition(mesh,
                               pysit.RickerWavelet(10.0),
                               sources=1,
                               source_depth=zdepth,
                               receivers='max',
                               receiver_depth=zdepth) 

print shots

[<pysit.core.shot.Shot object at 0x00000000099095C0>]

In [9]:

solver = pysit.ConstantDensityAcousticWave(mesh,
                                     formulation='scalar',
                                     model_parameters={'C': C},
                                     spatial_accuracy_order=2,
                                     trange=(0.0,3.0),
                                     use_cpp_acceleration=True)

Generate Synthetic Seismic Data¶

In [10]:

wavefields = []
base_model = solver.ModelParameters(mesh,{'C': C})
pysit.generate_seismic_data(shots, solver, base_model, wavefields=wavefields)

In [11]:

pysit.vis.display_seismogram(shots[0], cmap='gray', clim=(-0.01,0.01))

Setup Objective Function¶

In [12]:

objective = pysit.TemporalLeastSquares(solver)

Setup Optimization Method¶

In [13]:

invalg = pysit.LBFGS(objective)

initial_value = solver.ModelParameters(mesh,{'C': C0}) 

sc = {'residual_length_frequency' : 1,
      'objective_frequency' : 1,
      'gradient_frequency' : 1,
      'run_time_frequency' : 1,
                        }

In [15]:

nsteps = 1
result = invalg(shots, initial_value, nsteps, line_search='backtrack', 
                status_configuration=sc, verbose=True);

Iteration 0
  gradnorm 0.848689624313
  objective 0.000328374069781
  residual 0.0256270977592
  Starting:  45590.1551396 0.000328374069781
  Pass 1: a:11.1304089696; 0.000469260081457 ?<= 0.000319450882057
  Pass 2: a:8.90432717571; 0.000314273881233 ?<= 0.000322663229637
  alpha 8.90432717571
  run time 8.98700022697s

In [16]:

nsteps = 15
result = invalg(shots, initial_value, nsteps, line_search='backtrack', 
                status_configuration=sc);

In [17]:

pysit.pysit.vis.plot(result.C, mesh);

In [18]:

obj_vals = np.array([v for k,v in invalg.objective_history.items()])
plt.figure()
plt.semilogy(obj_vals, linewidth=2);

In [19]:

clim = C.min(),C.max()
plt.figure()
plt.subplot(3,1,1)
pysit.vis.plot(C0, mesh, clim=clim, cmap='gray')
plt.title('Initial Model')
plt.xticks([])
plt.subplot(3,1,2)
pysit.vis.plot(C, mesh, clim=clim, cmap='gray')
plt.title('True Model')
plt.xticks([])
plt.subplot(3,1,3)
pysit.vis.plot(result.C, mesh, clim=clim, cmap='gray')
plt.title('Reconstruction');

MSRI: Introduction to the Mathematics of Seismic Imaging

Python Seismic Imaging Toolbox¶

Full Waveform Inversion (FWI) + Numerical Python

Reproducible¶

Modular¶

Optimization¶

Objective Functions¶

Modeling

Solvers¶

Simple Rapid Prototyping¶

Simple Rapid Prototyping¶

2D Horizontal Reflector Demo¶

Define the Physical Domain¶

Define Computational Mesh¶

Define and Visualize Demo Problem¶

Setup Shots/Data Acquisition¶

Define a Solver¶

Generate Synthetic Seismic Data¶

Setup Objective Function¶

Setup Optimization Method¶

Run Optimization¶

Visualize Results¶

Visualize Results¶

Visualize Results¶