Invert three surveys line for a thin plate

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 matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import cofi

from vtem_max_forward_lib import (
    problem_setup,
    system_spec,
    survey_setup,
    ForwardWrapper,
    plot_predicted_profile,
    plot_transient,
    plot_plate_faces,
    plot_plate_faces_single
)

numpy.random.seed(42)

Background

This example inverts a three survey line of VTEM max data using the vertical component for a thin plate target. It thus becomes possible to invert for the easting, depth of the plate reference point, the plate dip and plate azimuth and the plate length. The forward problem and model setup is described here.

Problem definition

tx_min = 115
tx_max = 281
tx_interval = 15
ty_min = 25
ty_max = 176
ty_interval = 75
tx_points = numpy.arange(tx_min, tx_max, tx_interval)
ty_points = numpy.arange(ty_min, ty_max, ty_interval)
n_transmitters = len(tx_points) * len(ty_points)
tx, ty = numpy.meshgrid(tx_points, ty_points)
tx = tx.flatten()
ty = ty.flatten()

fiducial_id = numpy.arange(len(tx))
line_id = numpy.zeros(len(tx), dtype=int)
line_id[ty==ty_points[0]] = 0
line_id[ty==ty_points[1]] = 1
line_id[ty==ty_points[2]] = 2

survey_setup = {
    "tx": tx,                                                   # transmitter easting/x-position
    "ty": ty,                                                   # transmitter northing/y-position
    "tz": numpy.array([50]*n_transmitters),                     # transmitter height/z-position
    "tazi": numpy.deg2rad(numpy.array([90]*n_transmitters)),    # transmitter azimuth
    "tincl": numpy.deg2rad(numpy.array([6]*n_transmitters)),    # transmitter inclination
    "rx": tx,                                                   # receiver easting/x-position
    "ry": numpy.array([100]*n_transmitters),                    # receiver northing/y-position
    "rz": numpy.array([50]*n_transmitters),                     # receiver height/z-position
    "trdx": numpy.array([0]*n_transmitters),                    # transmitter receiver separation inline
    "trdy": numpy.array([0]*n_transmitters),                    # transmitter receiver separation crossline
    "trdz": numpy.array([0]*n_transmitters),                    # transmitter receiver separation vertical
    "fiducial_id": fiducial_id,                                 # unique id for each transmitter
    "line_id": line_id                  # id for each line
}

true_model = {
    "res": numpy.array([300, 1000]),
    "thk": numpy.array([20]),
    "peast": numpy.array([175]),
    "pnorth": numpy.array([100]),
    "ptop": numpy.array([30]),
    "pres": numpy.array([0.1]),
    "plngth1": numpy.array([100]),
    "plngth2": numpy.array([100]),
    "pwdth1": numpy.array([0.1]),
    "pwdth2": numpy.array([90]),
    "pdzm": numpy.array([75]),
    "pdip": numpy.array([60])
}

forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup,
                         ["pdip","pdzm", "peast", "ptop", "pwdth2"], data_returned=["vertical"])

['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2']

# check the order of parameters in a model vector
forward.params_to_invert

['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2']

true_param_value = numpy.array([60,65, 175, 30, 90])

True model

_, axes = plt.subplots(2, 2)
axes[1,1].axis("off")
plot_plate_faces(
    "plate_true", forward, true_param_value,
    axes[0,0], axes[0,1], axes[1,0], color="purple", label="True model"
)
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 0x7f02ca7f60d0>

Generate synthetic data

# The data
absolute_noise= 0.05

# create data and ad a realisation of the noise
data_pred_true = forward(true_param_value)
data_obs = data_pred_true + numpy.random.randn(len(data_pred_true))*absolute_noise

# define data covariance matrix
sigma=absolute_noise
Cdinv=numpy.identity(len(data_obs))*(1.0/(sigma*sigma))

Parameter estimation

Initialise a model for inversion

init_param_value = numpy.array([45, 90, 150, 20, 80])

Define helper functions for CoFI

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

def my_gradient(model):
    dpred = forward(model)
    jacobian = forward.jacobian(model, relative_step=0.1)
    residual = dpred - data_obs
    return jacobian.T @ Cdinv @ residual

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

class PerIterationCallbackFunction:
    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

Define CoFI problem

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)

Define CoFI options

my_options = cofi.InversionOptions()
my_options.set_tool("scipy.optimize.minimize")
my_options.set_params(method="Newton-CG",callback=PerIterationCallbackFunction())

Run CoFI inversion

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

Iteration #1
  objective value: 25216.587572321012
Iteration #2
  objective value: 7290.678246750972
Iteration #3
  objective value: 3101.2129903115647
Iteration #4
  objective value: 1982.6413889061635
Iteration #5
  objective value: 1684.5275371363728
Iteration #6
  objective value: 1585.5627896068984
Iteration #7
  objective value: 1569.5279365139831
Iteration #8
  objective value: 1553.5849324815508
Iteration #9
  objective value: 1552.8316200317665
Iteration #10
  objective value: 1551.64208272022
Iteration #11
  objective value: 1550.9006994110896
[ 59.89497772  65.24627406 174.87958828  29.9595983   90.03333672]

Plotting

Data - Fiducials

_, ax = plt.subplots()
plot_transient(true_param_value, forward, "data from true model", ax, color="purple")
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()

Data - Profiles

# Select gates to plot
idx_to_plot = numpy.arange(8, 30)

_, axes = plt.subplots(3, 1, figsize=(12,12))

for i in range(3):
    plot_predicted_profile(true_param_value, forward, "Data from true model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="purple")
    plot_predicted_profile(init_param_value, forward, "Data from starting model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="green", linestyle=":")
    plot_predicted_profile(my_result.model, forward, "Data from MAP model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="red", linestyle="-.")
    axes[i].set_title("Crossline {} (m)".format(ty_points[i]))
    axes[i].legend(bbox_to_anchor=(1.04, 0), loc="lower left")
plt.tight_layout()

Model

As a general rule we can only constrain parateters if there are fiducials in the survey that are ssnsitvit tothem and also fiducials that are not sensitive to them. In order words the anomly needs to be closed with respect to the model paratmeter in question. So to be able to constraint plate lenghy we would need survey lines to the north and south that are not overflying the thin plate.

_, axes = plt.subplots(2, 2)
axes[1,1].axis("off")
plot_plate_faces(
    "plate_true", forward, true_param_value,
    axes[0,0], axes[0,1], axes[1,0], color="purple", label="True model"
)
plot_plate_faces(
    "plate_init", forward, init_param_value,
    axes[0,0], axes[0,1], axes[1,0], color="green", label="Starting model"
)
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"
)

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 0x7f02ca7f5f90>

---

Watermark

watermark_list = ["cofi", "numpy", "scipy", "matplotlib","smt"]
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
smt 2.10.1

sphinx_gallery_thumbnail_number = -1