{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Effect of noise and height on the inversion\n", "This study has been demonstrated in the EMagPy paper [McLachlan et al. (2021)](https://doi.org/10.1016/j.cageo.2020.104561).\n", "\n", "All inversions are performed with the ROPE solver on a two-layer model with a varying depth. (a) Inversion with 0% noise with device on the ground. (b) Inversion with 5% noise on the ground. (c) Inversion with 0% noise at 1 m above the ground (d) Inversion with 5% noise at 1 m above the ground. The red line represents the true interface between the two layers." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import sys\n", "sys.path.append('../src') # add path where emagpy is\n", "from emagpy import Problem\n", "\n", "datadir = '../src/examples/'\n", "\n", "letters = ['a','b','c','d','e','f','g','h','i','j']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# parameters for the synthetic model\n", "nlayer = 2 # number of layer\n", "npos = 20 # number of sampling positions\n", "conds = np.ones((npos, nlayer))*[20, 100]\n", "x = np.linspace(0.1, 2, npos)[:,None]\n", "depths = 0.65 + 0.2 * np.sin(x*np.pi*2) # wave\n", "coils0 = ['VCP1.48f10000h0', 'VCP2.82f10000h0', 'VCP4.49f10000h0',\n", " 'HCP1.48f10000h0', 'HCP2.82f10000h0', 'HCP4.49f10000h0']\n", "coils1 = ['VCP1.48f10000h1', 'VCP2.82f10000h1', 'VCP4.49f10000h1',\n", " 'HCP1.48f10000h1', 'HCP2.82f10000h1', 'HCP4.49f10000h1']\n", "coils = [coils0, coils0, coils1, coils1]\n", "noises = [0, 0.05, 0, 0.05]\n", "ks = []\n", "# generate ECa using forward model\n", "for i in range(4):\n", " k = Problem()\n", " k.setModels([depths], [conds])\n", " _ = k.forward(forwardModel='FSlin', coils=coils[i], noise=noises[i])\n", " ks.append(k)\n", "\n", "# invert\n", "for k in ks:\n", " k.setInit(depths0=np.array([0.5]), fixedDepths=[False])\n", " k.invert(forwardModel='FSlin', method='ROPE', regularization='l1',\n", " bnds=[(0.05, 2.5), (5, 150), (5, 150)], rep=1000, njobs=-1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# graph of inverted values\n", "titles = ['(a) at 0 m (no noise)', '(b) at 0 m (with 5% noise)',\n", " '(c) at 1 m (no noise)', '(d) at 1 m (with 5% noise)']\n", "fig, axs = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(14,3))\n", "for i in range(4):\n", " ax = axs[i]\n", " ks[i].showResults(ax=ax, vmin=0, vmax=120, maxDepth=2, errorbar=True)\n", " rmseDepths = np.sqrt(np.sum((ks[i].depths[0][:,0] - depths[:,0])**2)/len(depths[:,0]))\n", " ax.set_title('{:s} \\nRMSE\$_{{depth}}\$={:.2f}'.format(titles[i], rmseDepths))\n", " ax.step(np.arange(depths.shape[0])+0.5, -depths, 'r', where='post') # true depth\n", " if i < 3:\n", " fig.axes[-1].remove()\n", " if i > 0:\n", " ax.set_ylabel('')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# graph of apparent values\n", "titles = ['at 0 m (no noise)', 'at 0 m (with 5% noise)',\n", " 'at 1 m (no noise)', 'at 1 m (with 5% noise)']\n", "fig, axs = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(14,3))\n", "for i in range(4):\n", " ax = axs[i]\n", " ks[i].show(ax=ax, vmin=10, vmax=100)\n", " ax.set_title('({:s}) {:s}'.format(letters[i], titles[i]))\n", " if i > 0:\n", " ax.get_legend().remove()\n", " ax.set_ylabel('')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }