{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n2D inversion of Loop-Loop EM Data\n=================================\n\nIn this example, we consider a single line of loop-loop EM data\nat 30kHz with 3 different coil separations [0.32m, 0.71m, 1.18m].\nWe will use only Horizontal co-planar orientations (vertical magnetic dipole),\nand look at the real and imaginary parts of the secondary magnetic field.\n\nWe use the :class:SimPEG.Maps.Surject2Dto3D mapping to invert for a 2D model\nand perform the forward modelling in 3D.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\nimport time\nfrom pymatsolver import Pardiso as Solver\n\nfrom SimPEG import (\n Mesh, Maps, Optimization,\n DataMisfit, Regularization, InvProblem,\n Inversion, Directives, Versions\n)\nfrom SimPEG.EM import FDEM\nfrom SimPEG.Utils import mkvc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Setup\n-----\n\nDefine the survey and model parameters\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sigma_surface = 10e-3\nsigma_deep = 40e-3\nsigma_air = 1e-8\n\ncoil_separations = [0.32, 0.71, 1.18]\nfreq = 30e3\n\nprint(\"skin_depth: {:1.2f}m\".format(500/np.sqrt(sigma_deep*freq)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a dipping interface between the surface layer and the deeper layer\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "z_interface_shallow = -0.25\nz_interface_deep = -1.5\nx_dip = np.r_[0., 8.]\n\n\ndef interface(x):\n interface = np.zeros_like(x)\n\n interface[x < x_dip] = z_interface_shallow\n\n dipping_unit = ((x >= x_dip) & (x <= x_dip))\n x_dipping = (\n -(z_interface_shallow-z_interface_deep)/x_dip\n )*(x[dipping_unit]) + z_interface_shallow\n interface[dipping_unit] = x_dipping\n\n interface[x > x_dip] = z_interface_deep\n\n return interface" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Forward Modelling Mesh\n----------------------\n\nHere, we set up a 3D tensor mesh which we will perform the forward\nsimulations on.\n\n

#### Note

In practice, a smaller horizontal discretization should be used to improve\n accuracy, particularly for the shortest offset (eg. you can try 0.25m).

\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "csx = 0.5 # cell size for the horizontal direction\ncsz = 0.125 # cell size for the vertical direction\npf = 1.3 # expansion factor for the padding cells\n\nnpadx = 7 # number of padding cells in the x-direction\nnpady = 7 # number of padding cells in the y-direction\nnpadz = 11 # number of padding cells in the z-direction\n\ncore_domain_x = np.r_[-11.5, 11.5] # extent of uniform cells in the x-direction\ncore_domain_z = np.r_[-2., 0.] # extent of uniform cells in the z-direction\n\n# number of cells in the core region\nncx = int(np.diff(core_domain_x)/csx)\nncz = int(np.diff(core_domain_z)/csz)\n\n# create a 3D tensor mesh\nmesh = Mesh.TensorMesh(\n [\n [(csx, npadx, -pf), (csx, ncx), (csx, npadx, pf)],\n [(csx, npady, -pf), (csx, 1), (csx, npady, pf)],\n [(csz, npadz, -pf), (csz, ncz), (csz, npadz, pf)]\n ]\n)\n# set the origin\nmesh.x0 = np.r_[\n -mesh.hx.sum()/2., -mesh.hy.sum()/2., -mesh.hz[:npadz+ncz].sum()\n]\n\nprint(\"the mesh has {} cells\".format(mesh.nC))\nmesh.plotGrid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inversion Mesh\n--------------\n\nHere, we set up a 2D tensor mesh which we will represent the inversion model\non\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inversion_mesh = Mesh.TensorMesh([mesh.hx, mesh.hz[mesh.vectorCCz<=0]])\ninversion_mesh.x0 = [-inversion_mesh.hx.sum()/2., -inversion_mesh.hy.sum()]\ninversion_mesh.plotGrid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mappings\n---------\n\nMappings are used to take the inversion model and represent it as electrical\nconductivity on the inversion mesh. We will invert for log-conductivity below\nthe surface, fixing the conductivity of the air cells to 1e-8 S/m\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create a 2D mesh that includes air cells\nmesh2D = Mesh.TensorMesh([mesh.hx, mesh.hz], x0=mesh.x0[[0, 2]])\nactive_inds = mesh2D.gridCC[:, 1] < 0 # active indices are below the surface\n\n\nmapping = (\n Maps.Surject2Dto3D(mesh) * # populates 3D space from a 2D model\n Maps.InjectActiveCells(mesh2D, active_inds, sigma_air) * # adds air cells\n Maps.ExpMap(nP=inversion_mesh.nC) # takes the exponential (log(sigma) --> sigma)\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "True Model\n----------\n\nCreate our true model which we will use to generate synthetic data for\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m_true = np.log(sigma_deep) * np.ones(inversion_mesh.nC)\ninterface_depth = interface(inversion_mesh.gridCC[:, 0])\nm_true[inversion_mesh.gridCC[:, 1] > interface_depth] = np.log(sigma_surface)\n\nfig, ax = plt.subplots(1, 1)\ncb = plt.colorbar(inversion_mesh.plotImage(m_true, ax=ax, grid=True), ax=ax)\ncb.set_label(\"$\\log(\\sigma)$\")\nax.set_title(\"true model\")\nax.set_xlim([-10, 10])\nax.set_ylim([-2, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Survey\n------\n\nCreate our true model which we will use to generate synthetic data for\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "src_locations = np.arange(-11, 11, 0.5)\nsrc_z = 0.25 # src is 0.25m above the surface\norientation = 'z' # z-oriented dipole for horizontal co-planar loops\n\n# reciever offset in 3D space\nrx_offsets = np.vstack([np.r_[sep, 0., 0.] for sep in coil_separations])\n\n# create our source list - one source per location\nsrcList = []\nfor x in src_locations:\n src_loc = np.r_[x, 0., src_z]\n rx_locs = src_loc - rx_offsets\n\n rx_real = FDEM.Rx.Point_bSecondary(\n locs=rx_locs, orientation=orientation, component='real'\n )\n rx_imag = FDEM.Rx.Point_bSecondary(\n locs=rx_locs, orientation=orientation, component='imag'\n )\n\n src = FDEM.Src.MagDipole(\n rxList=[rx_real, rx_imag], loc=src_loc, orientation=orientation,\n freq=freq\n )\n\n srcList.append(src)\n\n# create the survey and problem objects for running the forward simulation\nsurvey = FDEM.Survey(srcList)\nprob = FDEM.Problem3D_b(mesh, sigmaMap=mapping, Solver=Solver)\n\nprob.pair(survey)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Data\n----\n\nGenerate clean, synthetic data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "t = time.time()\ndclean = survey.dpred(m_true)\nprint(\n \"Done forward simulation. Elapsed time = {:1.2f} s\".format(time.time() - t)\n)\n\n\ndef plot_data(data, ax=None, color=\"C0\", label=\"\"):\n if ax is None:\n fig, ax = plt.subplots(1, 3, figsize=(15, 5))\n\n # data is [re, im, re, im, ...]\n data_real = data[0::2]\n data_imag = data[1::2]\n\n for i, offset in enumerate(coil_separations):\n ax[i].plot(\n src_locations, data_real[i::len(coil_separations)],\n color=color, label=\"{} real\".format(label)\n )\n ax[i].plot(\n src_locations, data_imag[i::len(coil_separations)], '--',\n color=color, label=\"{} imag\".format(label)\n )\n\n ax[i].set_title(\"offset = {:1.2f}m\".format(offset))\n ax[i].legend()\n ax[i].grid(which=\"both\")\n ax[i].set_ylim(np.r_[data.min(), data.max()]+1e-11*np.r_[-1, 1])\n\n ax[i].set_xlabel(\"source location x (m)\")\n ax[i].set_ylabel(\"Secondary B-Field (T)\")\n\n plt.tight_layout()\n return ax\n\nax = plot_data(dclean)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up data for inversion\n-------------------------\n\nWe will invert the clean data, and assign a standard deviation of 0.05, and\na floor of 1e-11.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "survey.std = 0.05\nsurvey.eps = 1e-11\nsurvey.dobs = dclean" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up the inversion\n--------------------\n\nWe create the data misfit, simple regularization\n(a Tikhonov-style regularization, :class:SimPEG.Regularization.Simple)\nThe smoothness and smallness contributions can be set by including\nalpha_s, alpha_x, alpha_y as input arguments when the regularization is\ncreated. The default reference model in the regularization is the starting\nmodel. To set something different, you can input an mref into the\nregularization.\n\nWe estimate the trade-off parameter, beta, between the data\nmisfit and regularization by the largest eigenvalue of the data misfit and\nthe regularization. Here, we use a fixed beta, but could alternatively\nemploy a beta-cooling schedule using :class:SimPEG.Directives.BetaSchedule\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dmisfit = DataMisfit.l2_DataMisfit(survey)\nreg = Regularization.Simple(inversion_mesh)\nopt = Optimization.InexactGaussNewton(maxIterCG=10, remember=\"xc\")\ninvProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)\n\nbetaest = Directives.BetaEstimate_ByEig(beta0_ratio=0.25)\ntarget = Directives.TargetMisfit()\n\ndirectiveList = [betaest, target]\ninv = Inversion.BaseInversion(invProb, directiveList=directiveList)\n\nprint(\"The target misfit is {:1.2f}\".format(target.target))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the inversion\n------------------\n\nWe start from a half-space equal to the deep conductivity.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m0 = np.log(sigma_deep) * np.ones(inversion_mesh.nC)\n\nt = time.time()\nmrec = inv.run(m0)\nprint(\"\\n Inversion Complete. Elapsed Time = {:1.2f} s\".format(time.time() - t))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the predicted and observed data\n------------------------------------\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 3, figsize=(15, 5))\nplot_data(dclean, ax=ax, color=\"C0\", label=\"true\")\nplot_data(invProb.dpred, ax=ax, color=\"C1\", label=\"predicted\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the recovered model\n------------------------\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 2, figsize = (12, 5))\n\n# put both plots on the same colorbar\nclim = np.r_[np.log(sigma_surface), np.log(sigma_deep)]\n\n# recovered model\ncb = plt.colorbar(\n inversion_mesh.plotImage(mrec, ax=ax, clim=clim), ax=ax,\n)\nax.set_title(\"recovered model\")\ncb.set_label(\"$\\log(\\sigma)$\")\n\n# true model\ncb = plt.colorbar(\n inversion_mesh.plotImage(m_true, ax=ax, clim=clim), ax=ax,\n)\nax.set_title(\"true model\")\ncb.set_label(\"$\\log(\\sigma)$\")\n\n# # uncomment to plot the true interface\n# x = np.linspace(-10, 10, 50)\n# [a.plot(x, interface(x), 'k') for a in ax]\n\n[a.set_xlim([-10, 10]) for a in ax]\n[a.set_ylim([-2, 0]) for a in ax]\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print the version of SimPEG and dependencies\n--------------------------------------------\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Versions()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Moving Forward\n--------------\n\nIf you have suggestions for improving this example, please create a pull request on the example in SimPEG _\n\nYou might try:\n - improving the discretization\n - changing beta\n - changing the noise model\n - playing with the regulariztion parameters\n - ...\n\n" ] } ], "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": 0 }