{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Checkout www.pygimli.org for more examples\n%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nVES inversion for a smooth model\n================================\n\nThis tutorial shows how an built-in forward operator is used for an Occam type\n(smoothness-constrained) inversion with fixed regularization (most natural).\nA direct current (DC) one-dimensional (1D) VES (vertical electric sounding)\nmodelling operator is used to generate data, add noise and inversion.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We import numpy numerics, mpl plotting, pygimli and the 1D plotting function\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nnp.random.seed(1337)\n\nimport matplotlib.pyplot as plt\n\nimport pygimli as pg\nfrom pygimli.viewer.mpl import drawModel1D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##############################################################################\n initialize the forward modelling operator and compute synthetic noisy data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "synres = [100., 500., 20., 800.] # synthetic resistivity\nsynthk = [4, 6, 10] # synthetic thickness (lay layer is infinite)\nab2 = np.logspace(-1, 2, 25) # 0.1 to 100 in 25 steps (8 points per decade)\nfBlock = pg.core.DC1dModelling(len(synres), ab2, ab2/3)\nrhoa = fBlock(synthk+synres)\n# The data are noisified using a\nerrPerc = 3. # relative error of 3 percent\nrhoa = rhoa * (np.random.randn(len(rhoa)) * errPerc / 100. + 1.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The forward operator can be called by f.response(model) or simply f(model)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "thk = np.logspace(-0.5, 0.5, 30)\nf = pg.core.DC1dRhoModelling(thk, ab2, ab2/3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create some transformations used for inversion\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "transRho = pg.trans.TransLogLU(1, 1000) # lower and upper bound\ntransRhoa = pg.trans.TransLog() # log transformation also for data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up inversion\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv = pg.core.Inversion(rhoa, f, transRhoa, transRho, False) # data vector, f, ...\n# The transformations can also be omitted and set individually by\n# inv.setTransData(transRhoa)\n# inv.setTransModel(transRho)\ninv.setRelativeError(errPerc / 100.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a homogeneous starting model\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model = pg.Vector(len(thk)+1, np.median(rhoa)) # uniform values\ninv.setModel(model) #" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set pretty large regularization strength and run inversion\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"inversion with lam=200\")\ninv.setLambda(100)\nres100 = inv.run() # result is a pg.Vector, but compatible to numpy array\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))\n# Decrease the regularization (smoothness) and start (from old result)\nprint(\"inversion with lam=20\")\ninv.setLambda(10)\nres10 = inv.run() # result is a pg.Vector, but compatible to numpy array\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))\n# We now optimize lambda such that data are fitted within noise (chi^2=1)\nprint(\"chi^2=1 optimized inversion\")\nresChi = inv.runChi1() # ends up in a lambda of about 3\nprint(\"optimized lambda value:\", inv.getLambda())\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show model (inverted and synthetic) as well as model response and data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(ncols=2, figsize=(8, 6)) # two-column figure\ndrawModel1D(ax[0], synthk, synres, color='b', label='synthetic',\n plot='semilogx')\ndrawModel1D(ax[0], thk, res100, color='g', label=r'$\\lambda$=100')\ndrawModel1D(ax[0], thk, res10, color='c', label=r'$\\lambda$=10')\ndrawModel1D(ax[0], thk, resChi, color='r',\n label=r'$\\chi$={:.2f}'.format(inv.getLambda()))\nax[0].grid(True, which='both')\nax[0].legend(loc='best')\nax[1].loglog(rhoa, ab2, 'rx-', label='measured')\nax[1].loglog(inv.response(), ab2, 'b-', label='fitted')\nax[1].set_ylim((max(ab2), min(ab2)))\nax[1].grid(True, which='both')\nax[1].set_xlabel(r'$\\rho_a$ [$\\Omega$m]')\nax[1].set_ylabel('AB/2 [m]')\nax[1].yaxis.set_label_position('right')\nax[1].yaxis.set_ticks_position('right')\nax[1].legend(loc='best')\nplt.show()" ] } ], "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 }