From fa979e0d0af460077b733b9e19932314a79fe1ed Mon Sep 17 00:00:00 2001 From: dccowan Date: Thu, 25 Jan 2024 13:17:05 -0800 Subject: [PATCH] copy inversion exercise to live folder --- live/02-inversion-dc-resistivity-2d.ipynb | 1163 +++++++++++++++++++++ 1 file changed, 1163 insertions(+) create mode 100644 live/02-inversion-dc-resistivity-2d.ipynb diff --git a/live/02-inversion-dc-resistivity-2d.ipynb b/live/02-inversion-dc-resistivity-2d.ipynb new file mode 100644 index 0000000..5b3cb63 --- /dev/null +++ b/live/02-inversion-dc-resistivity-2d.ipynb @@ -0,0 +1,1163 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2.5D DC Resistivity Inversion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This tutorial demonstrates DC resistivity inversion with SimPEG. We load normalized voltage data and 2D topography for a synthetic pole-dipole survey. Once loaded, we demonstrate a standard set of steps for recovering a model that characterizing subsurface electrical resistivity using weighted least-squares inversion.\n", + "\n", + "The following items are emphasized in this tutorial:\n", + "\n", + "- Assigning uncertainties to DC resistivity data.\n", + "- Designing a suitable mesh based on survey geometry.\n", + "- What model parameters we should invert for.\n", + "- Defining the inverse problem (data misfit, regularization, optimization).\n", + "- Stopping criteria for the inversion.\n", + "- Assessing inversion outputs.\n", + "\n", + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 0: Importing Modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# SimPEG functionality\n", + "from SimPEG.electromagnetics.static import resistivity as dc\n", + "from SimPEG.electromagnetics.static.utils.static_utils import (\n", + " plot_pseudosection,\n", + " apparent_resistivity_from_voltage,\n", + " pseudo_locations,\n", + ")\n", + "from SimPEG.utils.io_utils.io_utils_electromagnetics import read_dcip2d_ubc\n", + "from SimPEG.utils import model_builder, Counter\n", + "from SimPEG import (\n", + " maps,\n", + " data_misfit,\n", + " regularization,\n", + " optimization,\n", + " inverse_problem,\n", + " inversion,\n", + " directives,\n", + ")\n", + "\n", + "# discretize functionality\n", + "from discretize import TreeMesh\n", + "from discretize.utils import active_from_xyz\n", + "\n", + "# Basic Python functionality\n", + "import numpy as np\n", + "import matplotlib as mpl\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import LogNorm\n", + "\n", + "# Pooch for downloading files from the web\n", + "import pooch\n", + "\n", + "mpl.rcParams.update({\"font.size\": 14}) # default font size\n", + "cmap = mpl.cm.RdYlBu_r # default colormap" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Below is a class specifically used for this tutorial. It is not important to spent time reading through the code in this cell." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "class SaveInversionProgress(directives.InversionDirective, Counter):\n", + " \"\"\"\n", + " A custom directive to save items of interest during the course of an inversion\n", + " \"\"\"\n", + "\n", + " def initialize(self):\n", + " \"\"\"\n", + " This is called when we first start running an inversion\n", + " \"\"\"\n", + " # initialize an empty dictionary for storing results\n", + " self.inversion_results = {\n", + " \"iteration\": [],\n", + " \"beta\": [],\n", + " \"phi_d\": [],\n", + " \"phi_m\": [],\n", + " \"phi_m_small\": [],\n", + " \"phi_m_smooth_x\": [],\n", + " \"phi_m_smooth_z\": [],\n", + " \"dpred\": [],\n", + " \"model\": [],\n", + " }\n", + "\n", + " def endIter(self):\n", + " \"\"\"\n", + " This is run at the end of every iteration. So here, we just append\n", + " the new values to our dictionary\n", + " \"\"\"\n", + "\n", + " self.inversion_results[\"iteration\"].append(self.opt.iter)\n", + " self.inversion_results[\"beta\"].append(self.invProb.beta)\n", + " self.inversion_results[\"phi_d\"].append(2 * self.invProb.phi_d)\n", + " self.inversion_results[\"phi_m\"].append(2 * self.invProb.phi_m)\n", + " self.inversion_results[\"dpred\"].append(self.invProb.dpred)\n", + " self.inversion_results[\"model\"].append(self.invProb.model)\n", + "\n", + " reg = self.reg.objfcts[0]\n", + " phi_s = reg.objfcts[0](self.invProb.model) * reg.multipliers[0]\n", + " phi_x = reg.objfcts[1](self.invProb.model) * reg.multipliers[1]\n", + " phi_z = reg.objfcts[2](self.invProb.model) * reg.multipliers[2]\n", + "\n", + " self.inversion_results[\"phi_m_small\"].append(phi_s)\n", + " self.inversion_results[\"phi_m_smooth_x\"].append(phi_x)\n", + " self.inversion_results[\"phi_m_smooth_z\"].append(phi_z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step1: Load Data and Plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Download the Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "topo_url = \"https://github.com/simpeg/agrogeo24/raw/main/data/topo_2d.txt\"\n", + "topo_hash = \"1B2DAC5C46E48AB1C5606BB14621279C0A86DFD7790674A9FAAD83F3595E7030\"\n", + "\n", + "topo_filename = pooch.retrieve(topo_url, known_hash=topo_hash)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_url = \"https://github.com/simpeg/agrogeo24/raw/main/data/dc_data.obs\"\n", + "data_hash = \"EBA4F058B5189DB974E261F4C0B4C18B7841F29E6304F3138E94FACB436A2E58\"\n", + "\n", + "data_filename = pooch.retrieve(data_url, known_hash=data_hash)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load Topography" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "topo_2d = np.loadtxt(str(topo_filename))\n", + "print(topo_2d)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the Observed Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dc_data = read_dcip2d_ubc(data_filename, \"volt\", \"general\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot the Topography and Electrode Locations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "unique_locations = dc_data.survey.unique_electrode_locations\n", + "\n", + "fig = plt.figure(figsize=(10, 2))\n", + "ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])\n", + "ax.plot(topo_2d[:, 0], topo_2d[:, -1], color=\"b\", linewidth=1)\n", + "ax.plot(unique_locations[:, 0], unique_locations[:, -1], \"ro\", markersize=4)\n", + "ax.set_xlim([topo_2d[:, 0].min(), topo_2d[:, 0].max()])\n", + "ax.set_xlabel(\"x (m)\", labelpad=5)\n", + "ax.set_ylabel(\"y (m)\", labelpad=5)\n", + "ax.grid(True)\n", + "ax.set_title(\"Topography and Electrode Locations\", fontsize=16, pad=10)\n", + "plt.show(fig)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot Observed Data in Pseudo-Section" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot voltages pseudo-section\n", + "fig = plt.figure(figsize=(8, 2.75))\n", + "ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78])\n", + "plot_pseudosection(\n", + " dc_data,\n", + " plot_type=\"scatter\",\n", + " ax=ax1,\n", + " scale=\"log\",\n", + " cbar_label=\"V/A\",\n", + " scatter_opts={\"cmap\": mpl.cm.viridis},\n", + ")\n", + "ax1.set_title(\"Normalized Voltages\")\n", + "plt.show()\n", + "\n", + "# Get apparent conductivities from volts and survey geometry\n", + "apparent_resistivities = apparent_resistivity_from_voltage(dc_data.survey, dc_data.dobs)\n", + "\n", + "# Plot apparent resistivity pseudo-section\n", + "fig = plt.figure(figsize=(8, 2.75))\n", + "ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78])\n", + "plot_pseudosection(\n", + " dc_data.survey,\n", + " apparent_resistivities,\n", + " plot_type=\"contourf\",\n", + " data_locations=True,\n", + " ax=ax1,\n", + " scale=\"log\",\n", + " cbar_label=\"$\\Omega m$\",\n", + " mask_topography=True,\n", + " contourf_opts={\"levels\": 20, \"cmap\": mpl.cm.RdYlBu_r},\n", + ")\n", + "ax1.set_title(\"Apparent Resistivity\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 2: Assign Uncertainties \n", + "\n", + "\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Uncertainties are estimates of the standard deviations of the noise on our data. Assign uncertainties of 1e-7 V/A + 8 % to all data. We can do this by setting the *standard_deviation* property of our *data object*." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "inds_sort = np.argsort(np.abs(dc_data.dobs))\n", + "sorted_data = np.abs(dc_data.dobs[inds_sort])\n", + "sorted_std = dc_data.standard_deviation[inds_sort]\n", + "\n", + "fig = plt.figure(figsize=(7, 5))\n", + "ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])\n", + "\n", + "x_plot = np.arange(0, len(sorted_data))\n", + "ax.semilogy(x_plot, sorted_data, \"k.\", markersize=1)\n", + "\n", + "ax.set_title(\"Sorted Data and Uncertainties\")\n", + "ax.fill_between(x_plot, sorted_data - sorted_std, sorted_data + sorted_std, alpha=0.5)\n", + "ax.grid(alpha=0.5)\n", + "ax.set_ylabel(\"Observed Data (V/A)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 3: Designing a (Tree) Mesh\n", + "\n", + "\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Define and provide reasonable values for the following variables:\n", + "\n", + "* **dh** (minimum cell width)\n", + "* **dom_width_x** (width of the domain along x)\n", + "* **dom_width_y** (width of the domain along y)\n", + "\n", + "To do this, we access information from our data object. Here we:\n", + "\n", + "1. Obtain the unique electrode locations from the *dc_data* object\n", + "2. Use this to find the minimum and maximum electrode spacing\n", + "3. Use this to determine good values for *dh*, *dom_width_x* and *dom_width_y*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nbcx = 2 ** int(np.round(np.log(dom_width_x / dh) / np.log(2.0))) # num. base cells x\n", + "nbcy = 2 ** int(np.round(np.log(dom_width_y / dh) / np.log(2.0))) # num. base cells y\n", + "\n", + "# Define the base mesh with top at z = 0 m\n", + "hx = [(dh, nbcx)]\n", + "hy = [(dh, nbcy)]\n", + "mesh = TreeMesh([hx, hy], origin=\"CN\")\n", + "\n", + "# Shift top to maximum topography and center of survey line\n", + "y_topo_max = np.max(topo_2d[:, -1])\n", + "mesh.origin = mesh.origin + np.r_[np.median(topo_2d[:, 0]), y_topo_max]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Mesh refinement based on topography\n", + "mesh.refine_surface(\n", + " topo_2d[np.abs(topo_2d[:, 0]) < 150.0, :],\n", + " padding_cells_by_level=[0, 0, 4, 4],\n", + " finalize=False,\n", + ")\n", + "\n", + "# Extract unique electrode locations.\n", + "unique_locations = dc_data.survey.unique_electrode_locations\n", + "\n", + "# Mesh refinement near electrodes.\n", + "mesh.refine_points(\n", + " unique_locations, padding_cells_by_level=[10, 8, 8, 8, 4, 4], finalize=False\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Finalize the mesh\n", + "mesh.finalize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(8, 4))\n", + "\n", + "ax1 = fig.add_axes([0.14, 0.17, 0.8, 0.7])\n", + "mesh.plot_grid(ax=ax1, linewidth=1)\n", + "ax1.grid(False)\n", + "ax1.set_xlim(-200, 200)\n", + "ax1.set_ylim(y_topo_max - 200, y_topo_max)\n", + "ax1.set_title(\"Mesh\")\n", + "ax1.set_xlabel(\"x (m)\")\n", + "ax1.set_ylabel(\"y (m)\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 4: Define the Active Cells\n", + "\n", + "Use the [active_from_xyz](https://discretize.simpeg.xyz/en/main/api/generated/discretize.utils.active_from_xyz.html) utility function to find the indices of the active cells using the mesh and surface topography. The output quantity is a ``bool`` array." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "active_cells = active_from_xyz(mesh, topo_2d)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 5: Project Survey to Discretized Topography\n", + "\n", + "\n", + "\n", + "Use the [drape_electrodes_on_topography](https://docs.simpeg.xyz/content/api/generated/SimPEG.electromagnetics.static.resistivity.Survey.drape_electrodes_on_topography.html#SimPEG.electromagnetics.static.resistivity.Survey.drape_electrodes_on_topography) method to project electrodes to the discrete surface topography." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dc_data.survey.drape_electrodes_on_topography(mesh, active_cells, option=\"top\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 6: Mapping from the Model to the Mesh\n", + "\n", + "\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "1. Determine the number of active cells (i.e. the number of model parameters)\n", + "2. Generate an [SimPEG.maps.ExpMap](https://docs.simpeg.xyz/content/api/generated/SimPEG.maps.ExpMap.html#SimPEG.maps.ExpMap) mapping object to map log-resistivities to resistivities.\n", + "3. Generate a [SimPEG.maps.InjectActiveCells](https://docs.simpeg.xyz/content/api/generated/SimPEG.maps.InjectActiveCells.html#SimPEG.maps.InjectActiveCells) mapping object to project values on active cells to the entire mesh. Make sure to fix the value of inactive (air) cells to 1e8 $\\Omega m$.\n", + "4. Use the \\* operator to combine the separate mapping objects into a single mapping. Your final mapping should be called **log_resistivity_map**." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 7: Define the Forward Simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dc_simulation = dc.simulation_2d.Simulation2DNodal(\n", + " mesh, survey=dc_data.survey, rhoMap=log_resistivity_map, storeJ=True\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 8: Starting/Reference Models\n", + "\n", + "The **starting model** defines a reasonable starting point for the inversion.\n", + "\n", + "The **reference model** is used to include a-priori geological information in the recovered model.\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Use the median apparent resistivity value to define a homogeneous starting model on all active cells. Define this using the variable name **starting_model**. The number of elements must equal the number of model parameters, so we must first determine the number of active cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 9: Define the Inverse Problem\n", + "\n", + "Geophysical inversion is frequently formulated as an optimization problem. The solution to the inverse problem is the model $\\mathbf{m}$ which minimizes an objective (penalty) function $\\phi (\\mathbf{m})$ of the form:\n", + "\n", + "$$\n", + "\\min_{\\mathbf{m}} \\; \\phi (\\mathbf{m}) = \\phi_d (\\mathbf{m}) + \\beta \\, \\phi_m (\\mathbf{m})\\\\\n", + "$$\n", + "\n", + "where\n", + "\n", + "* $\\phi_d (\\mathbf{m})$ is the data misfit. It is a measure of how well data predicted with a model $\\mathbf{m}$ fits the data.\n", + "* $\\phi_m (\\mathbf{m})$ is the regularization. It is where we impose structural constraints on the recovered model.\n", + "* $\\beta$ is the trade-off parameter that balances the relative emphasis of fitting the data and imposing structure." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 9a: Define the Data Misfit\n", + "\n", + "The data misfit $\\phi_d (\\mathbf{m})$ is a measure of how well data predicted with a given model $\\mathbf{m}$ fits the observed data. Here, the data misfit is define as the L2-norm of a weighted residual:\n", + "\n", + "$$\n", + "\\phi_d (\\mathbf{m})\n", + "= \\sum_{i=1}^{nD} \\, \\Bigg ( \\frac{d_i^{pred} - d_i^{obs}}{\\varepsilon_i} \\Bigg )^2\n", + "= \\big \\| \\mathbf{W_d} \\big ( F[\\mathbf{m}] - \\mathbf{d^{obs}} \\big ) \\big \\|^2\n", + "$$\n", + "\n", + "We define $F[\\mathbf{m}]$ as the predicted data for a given model. And $\\mathbf{W_d}$ is a diagonal matrix that weights the residual by the data uncertainties.\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Generate a [data_misfit.L2DataMisfit](https://docs.simpeg.xyz/content/api/generated/SimPEG.data_misfit.L2DataMisfit.html) object and call it **dmis**. To instantiate this objected, we must set the following keyword arguments:\n", + "\n", + "* **simulation=dc_simulation** (our simulation object)\n", + "* **data=dc_data** (our data object)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 9b: Regularization\n", + "\n", + "The regularization function $\\phi_m (\\mathbf{m})$ is where we impose structural constraints on the recovered model. Here, we use a weighted least-squares regularization:\n", + "\n", + "$$\n", + "\\phi_m (\\mathbf{m})\n", + "= \\alpha_s \\, \\big \\| \\mathbf{W_s} \\big ( \\mathbf{m - m}_{ref} \\big ) \\big \\|^2\n", + "+ \\alpha_x \\, \\big \\| \\mathbf{W_x \\, G_x \\, m} \\big \\|^2\n", + "+ \\alpha_y \\, \\big \\| \\mathbf{W_z \\, G_y \\, m} \\big \\|^2\n", + "$$\n", + "\n", + "where \n", + "\n", + "* $\\alpha_s, \\alpha_x, \\alpha_y$ balance the relative emphasis on each term,\n", + "* $\\mathbf{W}_s, \\mathbf{W}_x, \\mathbf{W}_y$ are matrices that apply cell weights,\n", + "* $\\mathbf{G_x}$ and $\\mathbf{G_y}$ are partial gradient operators,\n", + "* and $\\mathbf{m}_{ref}$ is a user-defined reference model\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Generate a [regularization.WeightedLeastSquares](myst:SimPEG#SimPEG.regularization.WeightedLeastSquares) object and call it **reg**. Set to invert for the smoothest model by instantiating the object with:\n", + "\n", + "* **mesh** (the mesh as an input argument)\n", + "* **active_cells=active_cells** (regularization only acts on active cells)\n", + "* **reference_model=starting_model** (set starting model as the reference model)\n", + "* **alpha_s=1e-5** (set impact of the smallness term to be negligible)\n", + "\n", + "### Exercise (advanced):\n", + "\n", + "Set the regularization to enforce smoothness roughly 10x more than smallness. Do this by setting **length_scale_x** and **length_scale_y** to **10** as keyword arguments instead of setting **alpha_s**. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 9c: Optimization\n", + "\n", + "\n", + "\n", + "For a given trade-off parameter $\\beta$, we must solve the following optimization problem:\n", + "\n", + "$$\n", + "\\min_{\\mathbf{m}} \\; \\phi (\\mathbf{m}) = \\phi_d (\\mathbf{m}) + \\beta \\phi_m (\\mathbf{m}) \\\\\n", + "$$\n", + "\n", + "Because our data misfit and regularization functions are parabolic, the solution $\\mathbf{m^*}$ occurs when the gradient is equal to zero, i.e.:\n", + "\n", + "$$\n", + "\\nabla \\phi(\\mathbf{m^*}) = \\mathbf{0}\n", + "$$\n", + "\n", + "In practice, this problem is solved iteratively using quasi-Newton methods (see figure). Where $\\mathbf{H}$ represents an approximation of the Hessian, a model update $\\delta \\mathbf{m}$ direction is obtained by solving the following linear system:\n", + "\n", + "$$\n", + "\\mathbf{H}_k\\, \\delta \\mathbf{m}_k = - \\nabla \\! \\phi_k\n", + "$$\n", + "\n", + "We then determine an appropriate step length $\\gamma$ and update the model as follows:\n", + "\n", + "$$\n", + "\\mathbf{m}_{k+1} = \\mathbf{m}_k + \\gamma \\, \\delta \\mathbf{m}_k\n", + "$$\n", + "\n", + "This process is repeated until stopping criteria for the algorithm is reached.\n", + "\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Generate an [optimization.InexactGaussNewton](https://docs.simpeg.xyz/content/api/generated/SimPEG.optimization.InexactGaussNewton.html) object to solve the invert problem and call it **opt**. Set the following properties for the algorithm via keyword arguments:\n", + "\n", + "* **maxIter=40** (maximum number of iterations; beta and Newton)\n", + "* **maxIterLS=20** (maximum number of line searches for $\\gamma$)\n", + "* **maxIterCG=20** (maximum number of conjugate gradient iterations)\n", + "* **tolCG=1e-3** (minimum relative error threshold)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 9d: Inverse Problem\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Use the [BaseInvProblem](https://docs.simpeg.xyz/content/api/generated/SimPEG.inverse_problem.BaseInvProblem.html) class to fully define the inverse problem that is solved at each beta (trade-off parameter) iteration. Use the variable name **inv_prob**. Instantiation requires the following as input arguments:\n", + "\n", + "* **dmis** (the data misfit object)\n", + "* **reg** (the regularization object)\n", + "* **opt** (the optimization object)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 10: Defining Inversion Directives\n", + "\n", + "Directives provide specific instructions to the inversion while it is running. Here, we define SimPEG directives commonly used for least-squares DC resistivitiy inversion." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 10a: Beta Cooling Directives\n", + "\n", + "\n", + "\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "Define the starting beta, schedule and target misfit directives.\n", + "\n", + "* Generate a [directives.BetaEstimate_ByEig](myst:SimPEG#SimPEG.directives.BetaEstimate_ByEig) object and call it **starting_beta**. Set keyword argument **beta0_ratio=100**\n", + "* Generate a [directives.BetaSchedule](myst:SimPEG#SimPEG.directives.BetaSchedule) and call it **beta_schedule**. Set the following keyword arguments:\n", + "\n", + " * **coolingFactor=2** (reduction factor for $\\beta$)\n", + " * **coolingRate=2** (number of Newton iterations at each $\\beta$)\n", + "\n", + "* Generate a [directives.TargetMisfit](myst:SimPEG#SimPEG.directives.TargetMisfit) to set stopping criteria for the inversion and call it **target_misfit**. Set keyword argument **chifact=1** to so the inversion terminates when data misfit equals number of data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Other Directives\n", + "\n", + "Other common directives for weighted least-squares inversion of DC resistivity data are:\n", + "\n", + "- [UpdateSensitivityWeights](myst:SimPEG#SimPEG.directives.UpdateSensitivityWeights): apply sensitivity weighting to counteract the natural tendancy of DC resistivity inversion to place materials near the electrodes.\n", + "\n", + "- [UpdatePreconditioner](myst:SimPEG#SimPEG.directives.UpdatePreconditioner): Apply Jacobi preconditioner when solving the optimization problem to reduce the number of conjugate gradient iterations.\n", + "\n", + "- **SaveInversionProgress:** a directive we defined earlier in the notebook to allow us to better Q.C. the inversion result." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sensitivity_weights = directives.UpdateSensitivityWeights(\n", + " every_iteration=True, threshold_value=1e-2\n", + ")\n", + "update_jacobi = directives.UpdatePreconditioner(update_every_iteration=True)\n", + "save_inversion = SaveInversionProgress()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The directive objects are organized in a ``list``. The order of the list matters so the directives have been organized for you." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "directives_list = [\n", + " sensitivity_weights,\n", + " update_jacobi,\n", + " starting_beta,\n", + " beta_schedule,\n", + " save_inversion,\n", + " target_misfit,\n", + "]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 10: Define and Run the Inversion\n", + "\n", + "### Exercise (beginner):\n", + "\n", + "1. Define the inversion by instantiating an [inversion.BaseInversion](https://docs.simpeg.xyz/content/api/generated/SimPEG.inversion.BaseInversion.html) object. Call this object **inv**. This object requires the **inv_prob** and **directives_list** as input arguments.\n", + "2. Use the **run** method to run the inversion. This method requires the **starting_model** as in input argument." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 11: Inversion Outputs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 11a: Extract Some Inversion Outputs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iteration = save_inversion.inversion_results[\"iteration\"]\n", + "beta = save_inversion.inversion_results[\"beta\"]\n", + "phi_d = save_inversion.inversion_results[\"phi_d\"]\n", + "phi_m = save_inversion.inversion_results[\"phi_m\"]\n", + "dpred = save_inversion.inversion_results[\"dpred\"]\n", + "models = save_inversion.inversion_results[\"model\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 11b: Plot Tikhonov Curve" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15, 6))\n", + "ax1a = fig.add_axes([0.05, 0.05, 0.35, 0.9])\n", + "ax1b = ax1a.twinx()\n", + "ax2a = fig.add_axes([0.55, 0.05, 0.35, 0.9])\n", + "ax2b = ax2a.twinx()\n", + "\n", + "ax1a.plot(iteration, phi_d, \"b-o\")\n", + "ax1a.plot(\n", + " iteration,\n", + " target_misfit.chifact * dc_data.survey.nD * np.ones_like(iteration),\n", + " \"k--\",\n", + ")\n", + "ax1b.plot(iteration, phi_m, \"r-o\")\n", + "\n", + "ax1a.set_xlabel(\"Iteration\", fontsize=20)\n", + "ax1a.set_ylabel(\"Data Misfit ($\\phi_d$)\", color=\"b\", fontsize=18)\n", + "ax1a.tick_params(axis=\"y\", colors=\"b\")\n", + "ax1a.set_title(\"Tikhonov Curve\", fontsize=20)\n", + "ax1a.set_xlim([np.min(iteration), np.max(iteration)])\n", + "ax1b.tick_params(axis=\"y\", colors=\"r\")\n", + "ax1b.set_ylabel(\"Regularization ($\\phi_m$)\", color=\"r\", fontsize=18)\n", + "\n", + "ax2a.semilogy(iteration, phi_d, \"b-o\")\n", + "ax2a.semilogy(\n", + " iteration,\n", + " target_misfit.chifact * dc_data.survey.nD * np.ones_like(iteration),\n", + " \"k--\",\n", + ")\n", + "ax2b.semilogy(iteration, phi_m, \"r-o\")\n", + "ax2a.set_xlabel(\"Iteration\", fontsize=20)\n", + "ax2a.set_ylabel(\"Data Misfit ($\\phi_d$)\", color=\"b\", fontsize=18)\n", + "ax2a.tick_params(axis=\"y\", which=\"both\", colors=\"b\")\n", + "ax2a.set_title(\"Tikhonov Curve\", fontsize=20)\n", + "ax2a.set_xlim([np.min(iteration), np.max(iteration)])\n", + "ax2b.tick_params(axis=\"y\", which=\"both\", colors=\"r\")\n", + "ax2b.set_ylabel(\"Regularization ($\\phi_m$)\", color=\"r\", fontsize=18)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 11c: Choose an Iteration to Further Examine" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model_iteration = -1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 11d: Plot Observed Data, Prediced Data and Normalized Misfit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dobs = dc_data.dobs\n", + "std = dc_data.standard_deviation\n", + "\n", + "fig = plt.figure(figsize=(7, 9))\n", + "data_array = [\n", + " np.abs(dobs),\n", + " np.abs(dpred[model_iteration]),\n", + " (dobs - dpred[model_iteration]) / std,\n", + "]\n", + "plot_title = [\"Observed Data\", \"Predicted Data\", \"Normalized Misfit\"]\n", + "plot_units = [\"V/A\", \"V/A\", \"\"]\n", + "scale = [\"log\", \"log\", \"linear\"]\n", + "cmap_list = [mpl.cm.viridis, mpl.cm.viridis, mpl.cm.RdYlBu]\n", + "\n", + "ax1 = 3 * [None]\n", + "cax1 = 3 * [None]\n", + "cbar = 3 * [None]\n", + "cplot = 3 * [None]\n", + "\n", + "for ii in range(0, 3):\n", + " ax1[ii] = fig.add_axes([0.15, 0.72 - 0.33 * ii, 0.65, 0.21])\n", + " cax1[ii] = fig.add_axes([0.81, 0.72 - 0.33 * ii, 0.03, 0.21])\n", + " cplot[ii] = plot_pseudosection(\n", + " dc_data.survey,\n", + " data_array[ii],\n", + " \"contourf\",\n", + " data_locations=True,\n", + " ax=ax1[ii],\n", + " cax=cax1[ii],\n", + " scale=scale[ii],\n", + " cbar_label=plot_units[ii],\n", + " mask_topography=True,\n", + " contourf_opts={\"levels\": 25, \"cmap\": cmap_list[ii]},\n", + " )\n", + " ax1[ii].set_title(plot_title[ii])\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 11e: Plot Observed and Predicted Apparent Resistivities" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot apparent conductivity pseudo-section\n", + "fig = plt.figure(figsize=(8, 2.75))\n", + "ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78])\n", + "plot_pseudosection(\n", + " dc_data.survey,\n", + " dobs=apparent_resistivities,\n", + " plot_type=\"contourf\",\n", + " data_locations=True,\n", + " ax=ax1,\n", + " scale=\"log\",\n", + " cbar_label=\"$\\Omega m$\",\n", + " mask_topography=True,\n", + " contourf_opts={\"levels\": 40, \"cmap\": mpl.cm.RdYlBu_r},\n", + ")\n", + "ax1.set_title(\"Apparent Resistivity (Observed)\")\n", + "plt.show()\n", + "\n", + "# Plot apparent conductivity pseudo-section\n", + "apparent_resistivities_dpred = apparent_resistivity_from_voltage(\n", + " dc_data.survey, dpred[model_iteration]\n", + ")\n", + "\n", + "fig = plt.figure(figsize=(8, 2.75))\n", + "ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78])\n", + "plot_pseudosection(\n", + " dc_data.survey,\n", + " dobs=apparent_resistivities_dpred,\n", + " plot_type=\"contourf\",\n", + " data_locations=True,\n", + " ax=ax1,\n", + " scale=\"log\",\n", + " cbar_label=\"$\\Omega m$\",\n", + " mask_topography=True,\n", + " contourf_opts={\"levels\": 40, \"cmap\": mpl.cm.RdYlBu_r},\n", + ")\n", + "ax1.set_title(\"Apparent Resistivity (Predicted)\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 11f: Plot Recovered Model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "air_resistivity = 1e8\n", + "background_resistivity = 1e2\n", + "block_resistivity = 1e3\n", + "basement_resistivity = 1e1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define resistivity model\n", + "true_resistivity = background_resistivity * np.ones(mesh.n_cells)\n", + "\n", + "ind_basement = mesh.cell_centers[:, -1] < -16.0\n", + "true_resistivity[ind_basement] = basement_resistivity\n", + "\n", + "ind_block = model_builder.get_indices_block(\n", + " np.r_[5.0, -10.0], np.r_[15.0, 0.0], mesh.cell_centers\n", + ")\n", + "true_resistivity[ind_block] = block_resistivity\n", + "\n", + "true_resistivity = true_resistivity[active_cells]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Spatial plotting limits\n", + "x_lim = [-60, 60]\n", + "y_lim = [y_topo_max - 50, y_topo_max]\n", + "\n", + "# Define a mapping to plot models and ignore inactive cells\n", + "plotting_map = maps.InjectActiveCells(mesh, active_cells, np.nan)\n", + "plotting_model = [true_resistivity, np.exp(models[model_iteration])]\n", + "norm = LogNorm(vmin=1e1, vmax=1e3)\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\n", + "ax1 = 2 * [None]\n", + "ax2 = 2 * [None]\n", + "title_str = [\n", + " \"True Resistivity\",\n", + " \"Recovered Resistivity\",\n", + "]\n", + "\n", + "for ii in range(0, 2):\n", + " ax1[ii] = fig.add_axes([0.14, 0.55 - 0.5 * ii, 0.68, 0.35])\n", + " mesh.plot_image(\n", + " plotting_map * plotting_model[ii],\n", + " ax=ax1[ii],\n", + " grid=False,\n", + " pcolor_opts={\"norm\": norm, \"cmap\": mpl.cm.RdYlBu_r},\n", + " )\n", + "\n", + " ax1[ii].set_xlim(x_lim)\n", + " ax1[ii].set_ylim(y_lim)\n", + " ax1[ii].set_title(title_str[ii])\n", + " ax1[ii].set_xlabel(\"x (m)\")\n", + " ax1[ii].set_ylabel(\"y (m)\")\n", + "\n", + " ax2[ii] = fig.add_axes([0.84, 0.55 - 0.5 * ii, 0.03, 0.35])\n", + " cbar = mpl.colorbar.ColorbarBase(\n", + " ax2[ii], norm=norm, orientation=\"vertical\", cmap=mpl.cm.RdYlBu_r\n", + " )\n", + " cbar.set_label(r\"$\\rho (\\Omega m)$\", rotation=270, labelpad=15, size=12)\n", + "\n", + "pseudo_locations_xy = pseudo_locations(dc_data.survey)\n", + "ax1[1].scatter(pseudo_locations_xy[:, 0], pseudo_locations_xy[:, -1], 1.5, \"k\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}