4. Gradient of the misfit function

A basic example how to use the emg3d.optimize.gradient() routine to compute the adjoint-state gradient of the misfit function. Here we just show its usage.

For this example we use the survey and data as obtained in the example 3. Simulation.

import os
import pooch
import emg3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm, SymLogNorm


# Adjust this path to a folder of your choice.
data_path = os.path.join('..', 'download', '')

Load survey and data

First we load the survey and accompanying data as obtained in the example 3. Simulation.

fname = 'GemPy-II-survey-A.h5'
pooch.retrieve(
    'https://raw.github.com/emsig/data/2021-05-21/emg3d/surveys/'+fname,
    '5f2ed0b959a4f80f5378a071e6f729c6b7446898be7689ddc9bbd100a8f5bce7',
    fname=fname,
    path=data_path,
)
survey = emg3d.load(data_path + fname)['survey']

# Let's have a look
survey

Out:

Data loaded from «/home/dtr/Codes/emsig/emg3d-gallery/examples/download/GemPy-II-survey-A.h5»
[emg3d v0.17.1.dev22+ge23c468 (format 1.0) on 2021-04-14T22:24:03.229285].

Survey «GemPy-II Survey A»

<xarray.Dataset>
Dimensions:    (freq: 2, rec: 45, src: 3)
Coordinates:
  * src        (src) <U6 'TxED-1' 'TxED-2' 'TxED-3'
  * rec        (rec) <U7 'RxEP-01' 'RxEP-02' 'RxEP-03' ... 'RxEP-44' 'RxEP-45'
  * freq       (freq) <U3 'f-1' 'f-2'
Data variables:
    observed   (src, rec, freq) complex128 (-5.965512862515066e-14-5.52219071...
    synthetic  (src, rec, freq) complex128 (-1.821164983623268e-14-5.42793941...
Attributes:
    noise_floor:     1e-15
    relative_error:  0.05


We can see that the survey consists of three sources, 45 receivers, and two frequencies.

Create an initial model

To create an initial model we load the true model, but set all subsurface resistivities to 1 Ohm.m. So we are left with a homogeneous three-layer model air-seawater-subsurface, which includes the topography of the seafloor.

# Load true model
fname = "GemPy-II.h5"
pooch.retrieve(
    'https://raw.github.com/emsig/data/2021-05-21/emg3d/models/'+fname,
    'ea8c23be80522d3ca8f36742c93758370df89188816f50cb4e1b2a6a3012d659',
    fname=fname,
    path=data_path,
)
model = emg3d.load(data_path + fname)['model']
grid = model.grid

Out:

Data loaded from «/home/dtr/Codes/emsig/emg3d-gallery/examples/download/GemPy-II.h5»
[emg3d v1.0.0rc3.dev5+g0cd9e09 (format 1.0) on 2021-05-21T18:40:16.721968].
# Overwrite all subsurface resistivity values with 1.0
res = model.property_x
subsurface = (res > 0.5) & (res < 1000)
res[subsurface] = 1.0
model.property_x = res

# QC the initial model and the survey.
grid.plot_3d_slicer(model.property_x, xslice=12000, yslice=7000,
                    pcolor_opts={'norm': LogNorm(vmin=0.3, vmax=200)})

# Plot survey in figure above
fig = plt.gcf()
fig.suptitle('Initial resistivity model (Ohm.m)')
axs = fig.get_children()
rec_coords = survey.receiver_coordinates()
src_coords = survey.source_coordinates()
axs[1].plot(rec_coords[0], rec_coords[1], 'bv')
axs[2].plot(rec_coords[0], rec_coords[2], 'bv')
axs[3].plot(rec_coords[2], rec_coords[1], 'bv')
axs[1].plot(src_coords[0], src_coords[1], 'r*')
axs[2].plot(src_coords[0], src_coords[2], 'r*')
axs[3].plot(src_coords[2], src_coords[1], 'r*')
plt.show()
Initial resistivity model (Ohm.m)

Options for automatic gridding

gridding_opts = {
    'center': (src_coords[0][1], src_coords[1][1], -2200),
    'properties': [0.3, 10, 1, 0.3],
    'domain': (
        [rec_coords[0].min()-100, rec_coords[0].max()+100],
        [rec_coords[1].min()-100, rec_coords[1].max()+100],
        [-5500, -2000]
    ),
    'min_width_limits': (100, 100, 50),
    'stretching': (None, None, [1.05, 1.5]),
}

Create the Simulation

simulation = emg3d.simulations.Simulation(
    name="Initial Model",    # A name for this simulation
    survey=survey,           # Our survey instance
    model=model,             # The model
    gridding='both',         # Src- and freq-dependent grids
    max_workers=4,           # How many parallel jobs
    # solver_opts=...,       # Any parameter to pass to emg3d.solve
    gridding_opts=gridding_opts,
    receiver_interpolation='linear',  # For proper adjoint-state gradient
)

# Let's QC our Simulation instance
simulation

Simulation «Initial Model»

  • Survey «GemPy-II Survey A»: 3 sources; 45 receivers; 2 frequencies
  • Model: resistivity; isotropic; 100 x 100 x 102 (1,020,000)
  • Gridding: Frequency- and source-dependent grids; 160 x 40 x 48 (307,200) - 192 x 48 x 64 (589,824)


Compute Gradient

Out:

Compute efields:           0/6  [00:00]
Compute efields: #6        1/6  [01:32]
Compute efields: ########3 5/6  [01:42]
Compute efields: ##########6/6  [01:42]

Back-propagate :           0/6  [00:00]
Back-propagate : #6        1/6  [01:56]
Back-propagate : ########3 5/6  [02:05]
Back-propagate : ##########6/6  [02:05]

QC Gradient

# Set the gradient of air and water to NaN.
# This will eventually move directly into emgd3 (active and inactive cells).
grad[~subsurface] = np.nan


# Plot the gradient
grid.plot_3d_slicer(
        grad.ravel('F'), xslice=12000, yslice=7000, zslice=-4000,
        pcolor_opts={'cmap': 'RdBu_r',
                     'norm': SymLogNorm(
                         linthresh=1e-2, base=10, vmin=-1e1, vmax=1e1)}
        )

# Add survey
fig = plt.gcf()
fig.suptitle('Gradient of the misfit function')
axs = fig.get_children()
axs[1].plot(rec_coords[0], rec_coords[1], 'bv')
axs[2].plot(rec_coords[0], rec_coords[2], 'bv')
axs[3].plot(rec_coords[2], rec_coords[1], 'bv')
axs[1].plot(src_coords[0], src_coords[1], 'r*')
axs[2].plot(src_coords[0], src_coords[2], 'r*')
axs[3].plot(src_coords[2], src_coords[1], 'r*')
plt.show()
Gradient of the misfit function
Sat Jul 17 17:56:28 2021 CEST
OS Linux CPU(s) 4 Machine x86_64
Architecture 64bit RAM 15.5 GiB Environment Python
Python 3.9.4 | packaged by conda-forge | (default, May 10 2021, 22:13:33) [GCC 9.3.0]
numpy 1.21.0 scipy 1.7.0 numba 0.53.1
emg3d 1.1.0 empymod 2.1.2 xarray 0.18.2
discretize 0.7.0 h5py 3.3.0 matplotlib 3.4.2
tqdm 4.61.2
Intel(R) oneAPI Math Kernel Library Version 2021.2-Product Build 20210312 for Intel(R) 64 architecture applications


Total running time of the script: ( 4 minutes 12.030 seconds)

Estimated memory usage: 473 MB

Gallery generated by Sphinx-Gallery