Note
Go to the end to download the full example code.
10. Layered#
Since v1.8.0 there exists the possibility to use locally approximated,
layered models together with empymod to compute responses and the gradient
instead of the 3D model with emg3d.
Here we look at this possibility by reproducing the example 04. Gradient of the misfit function. (Please note, however, that the layered option does not only work for the gradient, but also for pure forward modelling.)
For this example we use the resistivity model created in the example GemPy-II: Perth Basin, together with the survey created in 03. Simulation.
Limitations#
The limitations of the layered modelling are:
Only point and dipole sources.
Only isotropic and VTI models.
Also, using the layered modelling has a few other implications, e.g.:
There are no {e;h}fields, only the fields at receiver locations.
gridding, most ofgridding_opts,solver_opts,receiver_interpolation, andfile_dirhave no effect.The attribute
emg3d.simulations.Simulation.jvecis not implemented.
Some considerations about speed#
If you are concerned about speed take into account that using layered computations changes the things one has to consider.
“Regular” 3D:
Linear with number of source-frequency pairs; parallelized over it.
Depends on gridding options; linear in number of cells.
Extra receivers at almost no extra cost.
“Layered” 1D:
Linear with number of source-receiver pairs; parallelized over sources; serial over receivers.
Only for the gradient: depends linearly on the number of layers in the provided model.
Extra frequencies at little extra cost.
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, data, and model#
First we load the model, survey, and accompanying data as obtained in the examples GemPy-II: Perth Basin and 03. Simulation (check these examples for more information).
isurvey = (
'GemPy-II-survey-A.h5',
'5f2ed0b959a4f80f5378a071e6f729c6b7446898be7689ddc9bbd100a8f5bce7',
'surveys',
)
imodel = (
'GemPy-II.h5',
'ea8c23be80522d3ca8f36742c93758370df89188816f50cb4e1b2a6a3012d659',
'models',
)
# Download model and survey.
for data in [isurvey, imodel]:
path = data[2]+'/'+data[0]
pooch.retrieve(
'https://raw.github.com/emsig/data/2021-05-21/emg3d/'+path,
data[1], fname=data[0], path=data_path,
)
# Load them.
survey = emg3d.load(os.path.join(data_path, isurvey[0]))['survey']
true_model = emg3d.load(os.path.join(data_path, imodel[0]))['model']
grid = true_model.grid
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].
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].
Create an initial model#
To create an initial model we take 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.
# Overwrite all subsurface resistivity values with 1.0
model = true_model.copy()
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.
popts = {'norm': LogNorm(vmin=0.3, vmax=200)}
grid.plot_3d_slicer(
model.property_x, xslice=12000, yslice=7500, pcolor_opts=popts
)
# 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*')
Options for layered computation#
If layered=True, one can provide a layered_opts-dict to a simulation.
For more information about this see towards the end of this example. The
simulation uses some defaults values, if they are not provided. By default a
local 1D model is generated by creating a right elliptic cylinder around
source and receiver position, and average the values for each layer within
this cylinder.
You can read more about the local extraction and the ellipse options in the
functions emg3d.models.Model.extract_1d and
emg3d.maps.ellipse_indices().
Note that using layered computations means that we do not have to provide any gridding options, as there is no gridding happening.
Create the Simulation#
simulation = emg3d.simulations.Simulation(
name="Initial Model",
survey=survey,
model=model,
max_workers=4,
layered=True, # Set layered to True!
# layered_opts=..., # Parameters for the local 1D extractions.
)
# Let's QC our Simulation instance
simulation
The simulation indicates that it uses layered computations with the method
'cylinder' and certain ellipse-parameters.
Compute Gradient#
Let’s compute the gradient (with layered=True!).
Compute empymod 0/3 [00:00]
Compute empymod ███▎ 1/3 [00:00]
Compute empymod ██████▋ 2/3 [00:00]
Compute empymod ██████████ 3/3 [00:00]
Compute empymod 0/3 [00:00]
Compute empymod ███▎ 1/3 [00:08]
Compute empymod ██████▋ 2/3 [00:12]
Compute empymod ██████████ 3/3 [00:12]
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=7500, 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*')
Even though we use layered (1D) models, the gradient looks three dimensional. This is because each source-receiver pair has its own cylinder, from which it generates a 1D model. The resulting gradient is then projected back to the same cylinder. This yields, summing up all source-receiver pairs, a three-dimensional gradient.
Note that this gradient is obtained using the finite-difference method, by
perturbing each layer by a tiny bit. This is in contrast to the gradient
obtained with emg3d, where the adjoint-state method is used.
Compare this gradient to the adjoint-state gradient, obtained by true 3D computations, in example 04. Gradient of the misfit function.
A few words about the 1D model extraction#
There are five possible methods: cylinder, prism, midpoint, source, and receiver. However, in reality they reduce to really two different possibility:
'midpoint': If midpoint is chosen, the selected 1D model is simply the column midway between the two provided points, usually source and receiver.'source'and'receiver'only differ in that both points are set to source or receiver, respectively, meaning that the column at source or receiver location is selected.'cylinder': If cylinder is chosen, a right elliptic cylinder around the two points is selected. How the ellipse is selected depends on a few other parameters, see functionsemg3d.models.Model.extract_1dandemg3d.maps.ellipse_indices(). The only difference for'prism'is that it uses the rectangular, coordinate-system aligned envelope of the ellipse.
The default values are: method='cylinder'. The corresponding default
values of the ellipse-options are factor=1.2, minor=0.8, and
radius is set to one skin depth using the lowest frequency and the
lowest conductivity value in the lowest layer in vertical direction.
Let’s look at the default value the simulation has chosen in this case:
lopts = simulation.layered_opts
lopts
So we see that it chose the default method 'cylinder', the default values
factor=1.2 and minor=0.8, and the radius is set to roughly 712 m.
We select now one source and one receiver, to QC how the selected cylinder looks like.
src = survey.sources['TxED-1']
rec = survey.receivers['RxEP-30']
layered, imat = true_model.extract_1d(
lopts['method'], src.center[:2], rec.center[:2], lopts['ellipse'],
return_imat=True
)
values = true_model.property_x.copy()
values[~np.array(imat, dtype=bool), :] = np.nan
xy = (src.center[:2] + rec.center[:2])/2
grid.plot_3d_slicer(
values, xslice=xy[0], yslice=xy[1], zslice=-4000, pcolor_opts=popts
)
fig = plt.gcf()
fig.suptitle('Right elliptic cylinder for this source and receiver location')
axs = fig.get_children()
axs[1].plot(rec.center[0], rec.center[1], 'bv')
axs[1].plot(src.center[0], src.center[1], 'r*')
axs[2].plot(rec.center[0], rec.center[2], 'bv')
axs[3].plot(rec.center[2], rec.center[1], 'bv')
axs[2].plot(src.center[0], src.center[2], 'r*')
axs[3].plot(src.center[2], src.center[1], 'r*')
We can compare the five different methods, and see what 1D models they extract from the 3D model.
fig, axs = plt.subplots(
1, 5, figsize=(9, 4), sharex=True, sharey=True, constrained_layout=True
)
axs = axs.ravel()
pdepth = np.repeat(true_model.grid.nodes_z[1:-1], 2)
ellipse = lopts['ellipse']
methods = ['source', 'receiver', 'midpoint', 'cylinder', 'prism']
for i, method in enumerate(methods):
axs[i].set_title(method)
p0 = src.center[:2]
p1 = rec.center[:2]
if method == 'source':
p1 = p0
method = 'midpoint'
elif method == 'receiver':
p0 = p1
method = 'midpoint'
oned = true_model.extract_1d(method, p0, p1, ellipse)
res = np.repeat(oned.property_x, 2)[1:-1]
axs[i].plot(res, pdepth)
axs[0].set_xscale('log')
axs[0].set_xlim([0.2, 60])
Ellipse#
As mentioned before, the ellipse is selected using the function
emg3d.maps.ellipse_indices(). The documentation of that function
contains a lot of useful information about it. However, here is a function
that illustrates the ellipse with all its parameters.
def plot_ellipse(ax, p0, p1, radius, factor=1.0, minor=1.0, check_foci=True):
"""Plot annotated explanation of emg3d.maps.ellipse_ind()."""
# Cast input points to numpy arrays.
p0, p1 = np.array(p0), np.array(p1)
# Center coordinates
cx, cy = (p0 + p1) / 2
# Adjacent and opposite sides
dx, dy = (p1 - p0) / 2
# c: linear eccentricity
dxy = np.linalg.norm([dx, dy])
# a: semi-major axis
major = max(dxy * factor, dxy + radius)
# b: semi-minor axis
minor = max(minor * major, radius)
if check_foci:
minor = max(minor, np.sqrt(abs(major**2 - dxy**2)))
# Distance to foci
df = np.sqrt(major**2 - minor**2)
# Take care of division by zero
if dy == 0.0:
cos, sin = 1.0, 0.0
else:
cos, sin = dx/dxy, dy/dxy
# Parametrized ellipse
t = np.linspace(0, 2*np.pi, 101)
cost, sint = np.cos(t), np.sin(t)
x = cx + major*cos*cost - minor*sin*sint
y = cy + major*sin*cost + minor*cos*sint
# Verteces
vr = [cx + major*cos, cy + major*sin]
vl = [cx - major*cos, cy - major*sin]
# Upper co-vertex
vu = [cx - minor*sin, cy + minor*cos]
# Foci
fr = [cx + df*cos, cy + df*sin]
fl = [cx - df*cos, cy - df*sin]
# PLOTTING
# Title with radius and _actual_ minor factor.
ax.set_title(f"r={radius:.1f}; "
f"f={factor:.2f}; "
f"m={1.0 if major == 0.0 else minor/major:.2f}", fontsize=10)
# Draw a, b, c, and r
ax.plot([cx, vr[0]], [cy, vr[1]], 'C4', label=r'a = c + r')
ax.plot([cx, vu[0]], [cy, vu[1]], 'C3', label=r'b = minor · a')
ax.plot([p0[0], cx], [p0[1], cy], 'C0', label=r'c = off / 2')
ax.plot([vl[0], p0[0]], [vl[1], p0[1]], 'C5', label=r'r = radius')
# Draw the two points and their circles with radius
for p in [p0, p1]:
ax.plot(*p, 'o', c='k', ms=8)
ax.plot(p[0] + radius*cost, p[1] + radius*sint, '.7')
# Draw the foci
ax.plot([fr[0], fl[0]], [fr[1], fl[1]], 'C2x', label='Foci')
# Draw the envelopping square
ax.plot([x.min(), x.max(), x.max(), x.min(), x.min()],
[y.min(), y.min(), y.max(), y.max(), y.min()], '.7')
# Draw the ellipse
ax.plot(x, y, 'C1')
# Ensure square axes
ax.axis('equal')
Here the annotation for the above selection. However, you can use this function to test the influence of the different parameters.
fig, ax = plt.subplots(constrained_layout=True)
plot_ellipse(
ax,
p0=src.center[:2],
p1=rec.center[:2],
radius=lopts['ellipse']['radius'],
factor=lopts['ellipse']['factor'],
minor=lopts['ellipse']['minor'],
)
ax.legend(bbox_to_anchor=(1.05, 0.7))
Total running time of the script: (0 minutes 15.331 seconds)
Estimated memory usage: 58 MB