Sensitivities

In this notebook we show the gradient with respect to the squared slowness computed with various physics. Gradeitns with respect to additional parameters will be added in the future to each propagator according to their physical parameters.

import os
os.environ['OMP_NUM_THREADS'] = '8'
os.environ['DEVITO_LOGGING'] = 'INFO'

from recipes import recipes_registry, layered_model
import matplotlib.pyplot as plt
import numpy as np

Model

For this illustrative purpose, we use a simple layered model and compute the gradient with a smooth background. We expect to obtain RTM-looking images of the layers

so = 16
model = layered_model('visco-vti', nlayers=3, shape=(101, 21, 101), spacing=(10, 10, 10),
                      density=True, nbl=40, space_order=so, fs=True)
model0 = layered_model('visco-vti', nlayers=3, shape=(101, 21, 101), spacing=(10, 10, 10), density=True,
                       space_order=so, dt=model.critical_dt, nbl=40, fs=True)
model0.smooth(('vp', 'b', 'epsilon', 'delta'), sigma=3)
iy = model.shape[1] // 2 + model.nbl
Operator `initdamp` ran in 0.03 s
Operator `initdamp` ran in 0.01 s
physics = ('iso-acoustic', 'tti-zhang')
options = {'nt': 1750, 'f0': 0.015, 'space_order': so}
z0 = 0 if model.fs else model.nbl
slicevp = model.vp.data[40:-40, iy, z0:-40]
plt.figure(figsize=(10, 10))
plt.subplot(2,2,1)
plt.imshow(slicevp.T, aspect="auto")
plt.title("P-wave velocity")
plt.colorbar()
plt.subplot(2,2,2)
plt.imshow(model.b.data[40:-40, iy, z0:-40].T, aspect="auto")
plt.title("Buoyancy")
plt.colorbar()
plt.subplot(2,2,3)
plt.imshow(model.epsilon.data[40:-40, iy, z0:-40].T, aspect="auto")
plt.title(r"$\epsilon$")
plt.colorbar()
plt.subplot(2,2,4)
plt.imshow(model.delta.data[40:-40, iy, z0:-40].T, aspect="auto")
plt.title(r"$\delta$")
plt.colorbar()
plt.tight_layout()

plt.figure(figsize=(10, 5))
plt.subplot(1,2,1)
plt.imshow(model0.vp.data[40:-40, iy, z0:-40].T, aspect="auto")
plt.title("P-wave background velocity")
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(model0.vp.data[40:-40, iy, z0:-40].T - model.vp.data[40:-40, iy, z0:-40].T, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
plt.title("dm")
plt.colorbar()
plt.tight_layout()

Inversion

for p in physics:
    plt.figure(figsize=(10, 10))
    solver = recipes_registry[p](model, options)
    solver0 = recipes_registry[p](model0, options)
    # Forward in true model
    solver.forward()
    # Forward in background model
    solver0.forward(save=True)
    plt.subplot(2,2,1)
    plt.imshow(solver.rec.data, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
    plt.title("True data")
    plt.subplot(2,2,2)
    plt.imshow(solver0.rec.data, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
    plt.title("Background data")
    plt.subplot(2,2,3)
    plt.imshow(solver.rec.data-solver0.rec.data, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
    plt.title("Residual")
    # Imaging
    solver0.rec.data[:] -= solver.rec.data[:]
    solver0.jacobian_adjoint()
    slicedm = np.array(solver0.perturbation(param='m').data[40:-40, iy, z0:-40])
    clip = np.quantile(np.abs(slicedm.reshape(-1)), 0.9)
    plt.subplot(2,2,4)
    plt.imshow(slicedm.T, aspect="auto", cmap="Greys", vmin=-clip, vmax=clip)
    plt.imshow(slicevp.T, aspect="auto", cmap="terrain", alpha=.5)
    plt.title("Gradient")
    plt.tight_layout()
    plt.show()
Operator `ForwardAcousticIsotropic` ran in 76.44 s
Operator `ForwardAcousticIsotropic` ran in 76.75 s
Operator `JacobianAdjAcousticIsotropic` ran in 135.44 s

Operator `ForwardZhangTTI` ran in 129.30 s
Operator `ForwardZhangTTI` ran in 131.43 s
Operator `JacobianAdjZhangTTI` ran in 194.76 s

Let’s run zhang considering it self-adjoint

plt.figure(figsize=(10, 10))
solver = recipes_registry['tti-zhang'](model, options)
solver0 = recipes_registry['tti-zhang'](model0, options)
# Forward in true model
solver.forward()
# Forward in background model
solver0.forward(save=True)
plt.subplot(2,2,1)
plt.imshow(solver.rec.data, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
plt.title("True data")
plt.subplot(2,2,2)
plt.imshow(solver0.rec.data, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
plt.title("Background data")
plt.subplot(2,2,3)
plt.imshow(solver.rec.data-solver0.rec.data, aspect="auto", vmin=-.1, vmax=.1, cmap="seismic")
plt.title("Residual")
# Imaging
solver0.rec.data[:] -= solver.rec.data[:]
solver0.jacobian_adjoint(true_adj=False)
slicedm = np.array(solver0.perturbation(param='m').data[40:-40, iy, z0:-40])
clip = np.quantile(np.abs(slicedm.reshape(-1)), 0.9)
plt.subplot(2,2,4)
plt.imshow(slicedm.T, aspect="auto", cmap="Greys", vmin=-clip, vmax=clip)
plt.imshow(slicevp.T, aspect="auto", cmap="terrain", alpha=.5)
plt.title("Gradient")
plt.tight_layout()
plt.show()
Operator `ForwardZhangTTI` ran in 128.25 s
Operator `ForwardZhangTTI` ran in 129.83 s
Operator `JacobianAdjZhangTTI` ran in 124.28 s

Back to top