Inversion of the Caber VTEM max survey for a thin plate target

If you are running this notebook locally, make sure you’ve followed steps here to set up the environment. (This environment.yml file specifies a list of packages required to run the notebooks)

# -------------------------------------------------------- #
#                                                          #
#     Uncomment below to set up environment on "colab"     #
#                                                          #
# -------------------------------------------------------- #

# !pip install -U cofi
# !pip install git+https://github.com/JuergHauser/PyP223.git
# !git clone https://github.com/inlab-geo/cofi-examples.git
# %cd cofi-examples/examples/vtem_max

# If this notebook is run locally PyP223 needs to be installed separately by uncommenting the following line,
# that is by removing the # and the white space between it and the exclamation mark.
# !pip install git+https://github.com/JuergHauser/PyP223.git

import numpy
import pandas
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import cofi
import json

from vtem_max_forward_lib import (
    ForwardWrapper,
    problem_setup,
    plot_field_data,
    plot_transient,
    plot_predicted_profile,
    plot_observed_profile,
    get_subset_of_survey,
    plot_plate_faces,
    plot_plate_faces_single,
    plot_survey_map
)

numpy.random.seed(42)

Background

(Prikhodko et. al,2010)

The Caber orebody is a volcanogenic massive sulfide (VMS) deposit in western Quebec, Canada, that hosts significant Zn-Cu-Ag mineralisation and VTEM max surveys have been collected over the target and previously been inverted recovering subvertical body dipping at about 80 degrees to the southwest. Here we invert the vertical compnent for a starting model with a dip azimuth of \(180 \degree\) and expect the inversion to change the dip azimuth to a southwesterly direction.

Further Reading

Prikhodko, A., Morrison, E., Bagrianski, A., Kuzmin, P., Tishin, P., & Legault, J. (2010). Evolution of VTEM? technical solutions for effective exploration. ASEG Extended Abstracts, 2010(1), 1-4.

McMillan, M. S., Schwarzbach, C., Haber, E., & Oldenburg, D. W. (2015). 3D parametric hybrid inversion of time-domain airborne electromagnetic data. Geophysics, 80(6), K25-K36.

Decimated survey data

# Load the data files
with open('caber_survey.npy', 'rb') as f:
    survey_setup = numpy.load(f,allow_pickle=True)
    system_spec = numpy.load(f,allow_pickle=True)
    data_obs = numpy.load(f,allow_pickle=True)

survey_setup=survey_setup.tolist()
system_spec=system_spec.tolist()

plot_survey_map(survey_setup)

_, axes = plt.subplots(5, 1, figsize=(7, 10))

for i in range(5):
    plot_observed_profile(survey_setup, data_obs, label="observed data", gate_idx=numpy.arange(0, 34, 2),
                              line_id=[i], ax=axes[i], color="pink")
plt.tight_layout()

# initial guess
initial_model = {
    "res": numpy.array([1000, 1000]),
    "thk": numpy.array([25]),
    "peast": numpy.array([710200]),
    "pnorth": numpy.array([5513450]),
    "ptop": numpy.array([150]),
    "pres": numpy.array([0.008]),
    "plngth1": numpy.array([150]),
    "plngth2": numpy.array([150]),
    "pwdth1": numpy.array([0.1]),
    "pwdth2": numpy.array([500]),
    "pdzm": numpy.array([225]),
    "pdip": numpy.array([80])
}

problem_setup['cellw']=100

forward = ForwardWrapper(initial_model, problem_setup, system_spec, survey_setup,
                         ["pdip","pdzm", "pwdth2","ptop","plngth1","plngth2"], data_returned=["vertical"])

['pdip', 'pdzm', 'pwdth2', 'ptop', 'plngth1', 'plngth2']

forward.params_to_invert

['pdip', 'pdzm', 'plngth1', 'plngth2', 'ptop', 'pwdth2']

init_param_value=[90,200,150,150,150,500]
init_param_value=[90,180,150,150,150,500]

Parameter estimation

def my_objective(model):
    dpred = forward(model)
    residual = dpred - data_obs
    return residual.T @ residual

def my_gradient(model):
    dpred = forward(model)
    jacobian = forward.jacobian(model)
    residual = dpred - data_obs
    return 2 * jacobian.T @ residual

def my_hessian(model):
    jacobian = forward.jacobian(model)
    return 2 * jacobian.T @ jacobian

class CallbackFunction:
    def __init__(self):
        self.x = None
        self.i = 0

    def __call__(self, xk):
        print(f"Iteration #{self.i+1}")
        print(f"  objective value: {my_problem.objective(xk)}")
        self.x = xk
        self.i += 1

my_problem = cofi.BaseProblem()
my_problem.set_objective(my_objective)
my_problem.set_gradient(my_gradient)
my_problem.set_hessian(my_hessian)
my_problem.set_initial_model(init_param_value)

my_options = cofi.InversionOptions()
my_options.set_tool("scipy.optimize.minimize")
#my_options.set_params(method="newton-cg", options={"maxiter": 5}, callback=CallbackFunction())
my_options.set_params(method="trust-ncg", options={"maxiter": 10}, callback=CallbackFunction())

my_inversion = cofi.Inversion(my_problem, my_options)
my_result = my_inversion.run()
print(my_result.model)

Iteration #1
  objective value: 3682.494724069316
Iteration #2
  objective value: 3652.5047048741376
Iteration #3
  objective value: 3598.26112692279
Iteration #4
  objective value: 3510.5582707397016
Iteration #5
  objective value: 3393.85818986716
Iteration #6
  objective value: 3337.42629970419
Iteration #7
  objective value: 3334.662725587539
Iteration #8
  objective value: 3334.662722409527
Iteration #9
  objective value: 3334.6627235490664
Iteration #10
  objective value: 3334.6627247255874
[ 89.08492569 206.93166544 135.4257979  142.5526461  199.57746436
 505.35124188]

my_result.message

'Maximum number of iterations has been exceeded.'

_, ax = plt.subplots()
plot_transient(init_param_value, forward, "data from init model", ax, color="green", linestyle=":")
plot_transient(my_result.model, forward, "data from inverted model", ax, color="red", linestyle="-.")
ax.legend(loc="upper center")
ax.set_title("vertical")
plt.tight_layout()

idx_to_plot = numpy.arange(8, 30)       # (0, 44)

_, axes = plt.subplots(5, 1, figsize=(7,10))

for i in range(5):
    plot_predicted_profile(init_param_value, forward, "data from init model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="green", linestyle=":")
    plot_predicted_profile(my_result.model, forward, "data from inverted model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="red", linestyle="-.")
ax.legend(loc="upper center")
plt.tight_layout()

Despite only a small improvement in fit to the data the recovered thing plate shows the expected southwesterly dip direction.

_, axes = plt.subplots(2, 2)
axes[1,1].axis("off")
plot_plate_faces(
    "plate_init", forward, init_param_value,
    axes[0,0], axes[0,1], axes[1,0], color="green", label="Starting model",
    surface_elevation=0
)
plot_plate_faces(
    "plate_inverted", forward, my_result.model,
    axes[0,0], axes[0,1], axes[1,0], color="red", label="MAP model", linestyle="dotted",
    surface_elevation=0
)

plt.tight_layout()
point = Line2D([0], [0], label='Fiducial', marker='o', markersize=5,
         markeredgecolor='orange', markerfacecolor='orange', linestyle='')

handles, labels = axes[1,0].get_legend_handles_labels()
handles.extend([point])

axes[1,0].legend(handles=handles,bbox_to_anchor=(1.04, 0), loc="lower left")

<matplotlib.legend.Legend object at 0x7f02c8bdd310>

---

Watermark

watermark_list = ["cofi", "numpy", "scipy", "matplotlib"]
for pkg in watermark_list:
    pkg_var = __import__(pkg)
    print(pkg, getattr(pkg_var, "__version__"))

cofi 0.2.11
numpy 2.3.5
scipy 1.17.0
matplotlib 3.10.8

sphinx_gallery_thumbnail_number = -1