3D Forward Simulation of Magnetic Gradiometry Data for Magnetic Vector Models

Author: Devin C. Cowan

Keywords: gradiometry, magnetic vector model, forward simulation, integral formulation, tree mesh.

Summary: Here we use the module SimPEG.potential_fields.magnetics to simulate magnetic gradiometry data for a magnetic vector model. We use the 3D integral formulation and carry out the forward simulation on a tree mesh.

Learning Objectives:

How to simulate magnetic data for 3D structures with SimPEG.
How to create magnetic gradiometry surveys; i.e. managing multiple data components.
How to design tensor meshes for magnetic simulations using the integral formulation.
How to construct a magnetic vector model.
How to predict magnetic gradiometry data for a magnetic vector model.
How to include surface topography in the forward simulation.
What are the units of the magnetic vector model and resulting data.

Import Modules¶

Here, we import all of the functionality required to run the notebook for the tutorial exercise. All of the functionality specific to simulating magnetic data are imported from SimPEG.potential_fields.magnetics. We also import some useful utility functions from SimPEG.utils. To simulate magnetic data, we need to define our problem geometry on a numerical grid (or mesh). To generate the mesh, we used the discretize package.

# SimPEG functionality
from SimPEG.potential_fields import magnetics
from SimPEG.utils import plot2Ddata, model_builder, mat_utils
from SimPEG import maps

# discretize functionality
from discretize import TreeMesh
from discretize.utils import mkvc, active_from_xyz

# Common Python functionality
import numpy as np
from scipy.interpolate import LinearNDInterpolator
from scipy.constants import mu_0
import matplotlib as mpl

mpl.rcParams.update({"font.size": 14})
import matplotlib.pyplot as plt
import os

save_output = False  # Optional

Define the Topography¶

Surface topography is defined as an (N, 3) numpy.ndarray for 3D simulations. Here, we create basic topography for the forward simulation. For user-specific simulations, you may load topography from an XYZ file.

[x_topo, y_topo] = np.meshgrid(np.linspace(-200, 200, 41), np.linspace(-200, 200, 41))
z_topo = (
    -15 * np.exp(-(x_topo**2 + y_topo**2) / 80**2)
    + 100.0
    + 0.5 * np.random.rand(*x_topo.shape)
)

fig = plt.figure(figsize=(6, 6))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection="3d")
ax.set_zlim([z_topo.max() - 40, z_topo.max()])
ax.plot_surface(x_topo, y_topo, z_topo, color="r", edgecolor="k", linewidth=0.5)
ax.set_box_aspect(aspect=None, zoom=0.85)
ax.set_xlabel("X (m)", labelpad=10)
ax.set_ylabel("Y (m)", labelpad=10)
ax.set_zlabel("Z (m)", labelpad=10)
ax.set_title("Topography (Exaggerated z-axis)", fontsize=16, pad=-20)
ax.view_init(elev=20.0)

x_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)
topo_xyz = np.c_[x_topo, y_topo, z_topo]

Define the Survey¶

Surveys within SimPEG generally require the user to create and connect three types of objects:

receivers: which define the locations of field measurements and type of data being measured.
sources: the passive or active sources responsible for generating geophysical responses, and their associated receivers.
survey: the object which stores and organizes all of the sources and receivers.

Here, we define the survey that will be used for the forward simulation. Magnetic surveys are simple to create. The user only needs an (N, 3) numpy.ndarray to define the xyz locations of the observation locations, the field components being measured, and the Earth’s inducing field. For the tutorial simulation, the receivers are located 10 m above the surface topography and spaced 10 m apart.

# Define the observation locations as an (N, 3) numpy array or load them.
x = np.linspace(-80.0, 80.0, 17)
y = np.linspace(-80.0, 80.0, 17)
x, y = np.meshgrid(x, y)
x, y = mkvc(x.T), mkvc(y.T)
fun_interp = LinearNDInterpolator(np.c_[x_topo, y_topo], z_topo)
z = fun_interp(np.c_[x, y]) + 10  # Flight height 10 m above surface.
receiver_locations = np.c_[x, y, z]

# Define the component(s) of the field we want to simulate as strings within
# a list. Here we measure the x, y and z derivatives of the Bz anomaly at
# each observation location.
components = ["bxz", "byz", "bzz"]

# Use the observation locations and components to define the receivers. To
# simulate data, the receivers must be defined as a list.
receiver_list = magnetics.receivers.Point(receiver_locations, components=components)
receiver_list = [receiver_list]

# Define the inducing field
field_inclination = 90   # inclination [deg]
field_declination = 0    # declination [deg]
field_amplitude = 50000  # amplitude [nT]

source_field = magnetics.sources.UniformBackgroundField(
    receiver_list=receiver_list,
    amplitude=field_amplitude,
    inclination=field_inclination,
    declincation=field_declination
)

# Define the survey
survey = magnetics.survey.Survey(source_field)

If desired, we can extract various objects and properties from the objects used to generate the survey. E.g.

print("# of locations: {}".format(survey.nRx))  # number of receiver locations
print("# of data: {}".format(survey.nD))  # number of data that will be simulated
print(survey.source_field.inclination)  # inclination of the source field
print(survey.source_field.receiver_list[0])  # the receiver object
print(receiver_list[0].locations[:5, :])  # the first 5 receiver locations

# of locations: 289
# of data: 867
90.0
<SimPEG.potential_fields.magnetics.receivers.Point object at 0x7fd15b077160>
[[-80.         -80.         108.10265364]
 [-70.         -80.         107.59535864]
 [-60.         -80.         107.19181796]
 [-50.         -80.         106.60620485]
 [-40.         -80.         106.00709073]]

Design a (Tree) Mesh¶

Meshes are designed using the discretize package. See the discretize user tutorials to learn more about creating meshes. Here, the forward simulation is computed for a tree mesh. If you wanted to generate a tensor mesh instead, you can use the code snippet from the 3D Forward Simulation of Total Magnetic Intensity Data tutorial.

The integral formulation for magnetics essentially sums the independent contribution for every magnetized voxel cell in the mesh. Since the kernel function that computes the contribution for a single cell is an analytic solution, small cells are not required to accurately compute the contributions from coarse structures with constant magnetization (e.g. a rectangular prism). For complex structures however, or to define surface topography more accurately, finer cells may be needed. Furthermore, cells do not need to be cubic. Since the analytic solution is only valid outside the magnetized region, please do no place receivers within the Earth.

Here, a core cell width of 5 m is used within our survey region. Padding is used to extend the mesh outside the immediate survey area. The mesh will be plotted after we define our density contrast model.

dx = 5  # minimum cell width (base mesh cell width) in x
dy = 5  # minimum cell width (base mesh cell width) in y
dz = 5  # minimum cell width (base mesh cell width) in z

x_length = 240.0  # domain width in x
y_length = 240.0  # domain width in y
z_length = 120.0  # domain width in z

# Compute number of base mesh cells required in x and y
nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))
nbcz = 2 ** int(np.round(np.log(z_length / dz) / np.log(2.0)))

# Define the base mesh. Top defined at z = 0 m.
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
hz = [(dz, nbcz)]
mesh = TreeMesh([hx, hy, hz], x0="CCN")

# Shift vertically to top same as maximum topography
mesh.origin += np.r_[0.0, 0.0, z_topo.max()]

# Refine based on surface topography
mesh.refine_surface(topo_xyz, padding_cells_by_level=[2, 2], finalize=False)

# Refine box based on region of interest
wsb_corner = np.c_[-100, -100, 20]
ent_corner = np.c_[100, 100, 100]
# Note -1 is a flag for smallest cell size
mesh.refine_box(wsb_corner, ent_corner, levels=[-1], finalize=False)

mesh.finalize()

If desired, we can extract various properties of the mesh. E.g.

print("# of cells: {}".format(mesh.n_cells))  # Number of cells
print("# of x-faces: {}".format(mesh.n_faces_x))  # Number of x-faces
print("Origin: {}".format(mesh.origin))  # bottom-southewest corner
print("Max cell volume: {}".format(mesh.cell_volumes.max()))  # Largest cell size
print(mesh.cell_centers[0:5, :])  # Cell center locations

# of cells: 42032
# of x-faces: 41884
Origin: [-160.         -160.          -59.50188146]
Max cell volume: 64000.0
[[-140.         -140.          -39.50188146]
 [-100.         -140.          -39.50188146]
 [-140.         -100.          -39.50188146]
 [-100.         -100.          -39.50188146]
 [-140.         -140.            0.49811854]]

Define the Active Cells¶

Whereas cells below the Earth’s surface contribute towards the simulated magnetic anomaly, air cells do not. The set of mesh cells used in the forward simulation are referred to as ‘active cells’. Unused cells (air cells) are ‘inactive cells’. Here, the discretize active_from_xyz utility function is used to find the indices of the active cells using the mesh and surface topography. The output quantity is a bool array.

# Indices of the active mesh cells from topography (e.g. cells below surface)
active_cells = active_from_xyz(mesh, topo_xyz)

Mapping from the Model to Active Cells¶

In SimPEG, the term ‘model’ is not synonymous with the physical property values defined on the mesh. For example, the model may be defined as the logarithms of the physical property values, or be parameters defining a layered Earth. When simulating magnetic data using the integral formulation, we must define a mapping from the set of model parameters to the active cells in the mesh. Mappings are created using the SimPEG.maps module. For the tutorial exercise, the model is the density contrast values for all active cells. As such, our mapping is an identity mapping, whose dimensions are equal to the number of active cells.

# Define mapping from model to active cells. The model consists of a density
# contrast value for each cell below the Earth's surface.
n_active = int(active_cells.sum())
model_map = maps.IdentityMap(nP=3 * n_active)

Define the Magnetic Vector Model¶

Magnetic vector models are defined by three-component effective susceptibilities. Where $\vec{M}_{ind}$ represents the induced magnetization, $\vec{M}_{rem}$ represents the remanent magnetization, and $\vec{H}_0$ is the inducing magnetic field intensity, the magnetic vector model defines the components of:

\vec{\chi}_{eff} = \frac{\vec{M}_{ind} + \vec{M}_{rem}}{\left | \vec{H}_0 \right |}

(1)

As a result, magnetic vector models in SimPEG are defined using SI units. To create a magnetic vector model, we must

Compute the induced magnetization vector for each cell by multiplying the susceptibility χ by the inducing magnetic field intensity $\vec{H}_0$ . The resulting quantity has units A/m.
Define the remanent magnetization vector for each cell, once again in A/m.
Sum the induced and remanent contributions.
Normalize by the amplitude of inducing magnetic field intensity; i.e. $\left | \vec{H}_0 \right |$ in units A/m. The resulting quantity is unitless (SI).
Re-organize the quantity as a 1D numpy.ndarray organized by vector component; i.e. np.r_[chi_1, chi_2, chi_3]

Here, the model consists of a susceptible and remanently magnetized sphere within a negligibly susceptible host. The horizontal components of the induced and remanent magnetization has been balance so that they cancel out, and the net magnetization is downward. We plot the magnetic vector models using the plot_slice method.

# Define susceptibility values for each unit in SI
background_susceptibility = 0.0001
sphere_susceptibility = 0.01

# Compute the induced magnetization vector (A/m) for every active cell
susceptibility_model = background_susceptibility * np.ones(n_active)
ind_sphere = model_builder.get_indices_sphere(
    np.r_[0.0, 0.0, 55.0], 16.0, mesh.cell_centers
)
ind_sphere = ind_sphere[active_cells]
susceptibility_model[ind_sphere] = sphere_susceptibility

# Compute the unit direction of the inducing field in Cartesian coordinates
field_direction = mat_utils.dip_azimuth2cartesian(field_inclination, field_declination)

# Inducing magnetic field intensity (A/m)
H0 = 1e-9 * field_amplitude * field_direction / mu_0

# Compute induced magnetization
induced_magnetization = np.outer(
    susceptibility_model, H0
)  # (n_active, 3) numpy.ndarray

# Define the remanent magnetization vector (A/m) for every active cell.
# inclination (deg), declination (deg), amplitude (A/m)
remanence_inclination = 45.0
remanence_declination = 210.0
remanence_amplitude = 0.39788735751313814

remanent_magnetization = np.zeros_like(induced_magnetization)
remanent_magnetization_sphere = remanence_amplitude * mat_utils.dip_azimuth2cartesian(
    remanence_inclination, remanence_declination
)
remanent_magnetization[ind_sphere, :] = remanent_magnetization_sphere

# Compute total magnetization (A/m) for all active cells
total_magnetization = induced_magnetization + remanent_magnetization

# Define effective susceptibility model as a vector np.r_[chi_x, chi_y, chi_z]
model = mkvc(total_magnetization) / np.linalg.norm(H0)

# Mapping to ignore inactive cells when plotting MVI amplitude
plotting_map = maps.InjectActiveCells(mesh, active_cells, np.nan)

# MVI model organized into (N, 3) numpy.ndarray
plotting_model = total_magnetization / np.linalg.norm(H0)

# Amplitude of MVI
model_amplitude = np.sqrt(np.sum(plotting_model**2, axis=1))

fig = plt.figure(figsize=(14, 4))

norm = mpl.colors.Normalize(vmin=0, vmax=np.max(model_amplitude))

ax1 = fig.add_axes([0.05, 0.12, 0.5, 0.78])
mesh.plot_slice(
    plotting_map * model_amplitude,
    normal="Y",
    ax=ax1,
    ind=int(mesh.h[1].size / 2),
    grid=True,
    pcolor_opts={"cmap": mpl.cm.plasma, "norm": norm},
)
ax1.set_title("MVI Model at y = 0 m (amplitude)")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")

ax2 = fig.add_axes([0.62, 0.12, 0.25, 0.78])
mesh.plot_slice(
    plotting_map * plotting_model,
    v_type="CCv",
    view="vec",
    normal="Y",
    ax=ax2,
    ind=int(mesh.h[1].size / 2),
    grid=True,
    pcolor_opts={"cmap": mpl.cm.plasma, "norm": norm},
    quiver_opts={
        "pivot": "mid",
        "width": 0.01,
        "headwidth": 3.0,
        "headlength": 3.0,
        "headaxislength": 3.0,
        "scale": 0.25,
    },
)
ax2.set_title("MVI Model at y = 0 m")
ax2.set_xlim([-25, 25])
ax2.set_ylim([30, 80])
ax2.set_xlabel("x (m)")
ax2.set_ylabel("z (m)")

cx = fig.add_axes([0.89, 0.12, 0.02, 0.78])
cbar = mpl.colorbar.ColorbarBase(
    cx, norm=norm, orientation="vertical", cmap=mpl.cm.plasma
)
cbar.set_label("SI", rotation=270, labelpad=15)

Define the Forward Simulation¶

In SimPEG, the physics of the forward simulation is defined by creating an instance of an appropriate simulation class. In this case, we use the simulation class for the 3D integral formulation. To fully define the forward simulation, we need to connect the simulation object to:

the survey
the mesh
the indices of the active cells
the mapping from the model to the active cells
the model type: “scalar” for susceptibility model, “vector” for magnetic vector model

This is accomplished by setting each one of the aforementioned items as a property of the simulation object. Additional keyword arguments can also be set which impact the forward simulation. Because we are only simulating data for a single model, there is no benefit to storing the sensitivities for the forward simulation. store_sensitivities property to ‘forward_only’

simulation = magnetics.simulation.Simulation3DIntegral(
    survey=survey,
    mesh=mesh,
    model_type="vector",
    chiMap=model_map,
    ind_active=active_cells,
    store_sensitivities="forward_only",
)

Simulate Magnetic Gradiometry Data¶

Once any simulation within SimPEG has been properly constructed, simulated data for a given model vector can be computed using the dpred method. In SimPEG, magnetic gradiometry data values are in nT/m.

dpred = simulation.dpred(model)

Data are ordered by component, then by location. Here, we reshape the predicted data vector into an array for easier plotting.

n_loc = survey.nRx  # Number of receiver locations
n_comp = len(components)  # Number of data components
dpred_plotting = np.reshape(dpred, (n_loc, n_comp))

fig = plt.figure(figsize=(10, 3))
v_max = np.max(np.abs(dpred))

ax = 3 * [None]
cplot = 3 * [None]
comp_list = ["x", "y", "z"]

norm = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)

for ii in range(0, 3):
    ax[ii] = fig.add_axes([0.1 + ii * 0.26, 0.15, 0.25, 0.78])
    cplot[ii] = plot2Ddata(
        receiver_locations,
        dpred_plotting[:, ii],
        ax=ax[ii],
        ncontour=60,
        contourOpts={"cmap": "bwr", "norm": norm},
    )
    ax[ii].set_title(r"$\partial B_z /\partial {}$".format(comp_list[ii]))
    ax[ii].set_xlabel("x (m)")
    if ii == 0:
        ax[ii].set_ylabel("y (m)")
    else:
        ax[ii].set_yticks([])

cx = fig.add_axes([0.89, 0.13, 0.02, 0.79])
cbar = mpl.colorbar.ColorbarBase(cx, norm=norm, orientation="vertical", cmap=mpl.cm.bwr)
cbar.set_label("$nT/m$", rotation=270, labelpad=10, size=12)

plt.show()

SimPEG User Tutorials

3D Forward Simulation of TMI Data

SimPEG User Tutorials

3D Inversion of TMI Data to Recover a Susceptibility Model