{ "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 blocky model\n================================\n\nThis tutorial shows how an built-in forward operator is used for inversion.\nA DC 1D (VES) modelling is used to generate data, noisify and invert them.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We import numpy, matplotlib and the 1D plotting function\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\n\nimport pygimli as pg\nfrom pygimli.viewer.mpl import drawModel1D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "some definitions before (model, data and error)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nlay = 4 # number of layers\nlam = 200. # (initial) regularization parameter\nerrPerc = 3. # relative error of 3 percent\nab2 = np.logspace(-1, 2, 50) # AB/2 distance (current electrodes)\nmn2 = ab2 / 3. # MN/2 distance (potential electrodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "initialize the forward modelling operator\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "f = pg.core.DC1dModelling(nlay, ab2, mn2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "other ways are by specifying a Data Container or am/an/bm/bn distances\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "synres = [100., 500., 20., 800.] # synthetic resistivity\nsynthk = [0.5, 3.5, 6.] # synthetic thickness (nlay-th layer is infinite)" ] }, { "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": [ "rhoa = f(synthk+synres)\nrhoa = rhoa * (pg.randn(len(rhoa), seed=0) * errPerc / 100. + 1.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "create some transformations used for inversion\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "transThk = pg.trans.TransLog() # log-transform ensures thk>0\ntransRho = pg.trans.TransLogLU(1, 1000) # lower and upper bound\ntransRhoa = pg.trans.TransLog() # log transformation for data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "set model transformation for thickness and resistivity\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "f.region(0).setTransModel(transThk) # 0=thickness\nf.region(1).setTransModel(transRho) # 1=resistivity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "generate start model values from median app. resistivity & spread\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "paraDepth = max(ab2) / 3. # rule-of-thumb for Wenner/Schlumberger\nf.region(0).setStartValue(paraDepth / nlay / 2)\nf.region(1).setStartValue(np.median(rhoa))" ] }, { "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, True) # data vector, fop, verbose\n# could also be set by inv.setTransData(transRhoa)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "set error model, regularization strength and Marquardt scheme\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRelativeError(errPerc / 100.0) # alternative: setAbsoluteError in Ohmm\ninv.setLambda(lam) # (initial) regularization parameter\ninv.setMarquardtScheme(0.9) # decrease lambda by factor 0.9\nmodel = f.createStartVector() # creates from region start value\nmodel[nlay] *= 1.5 # change default model by changing 2nd layer resistivity\ninv.setModel(model) #" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "run actual inversion and extract resistivity and thickness\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model = inv.run() # result is a pg.Vector, but compatible to numpy array\nres, thk = model[nlay-1:nlay*2-1], model[0:nlay-1]\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "show estimated&synthetic models and data with model response in 2 subplots\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, plot='semilogx', color='r')\ndrawModel1D(ax[0], thk, res, color='b')\nax[0].grid(True, which='both')\nax[0].set_ylabel('z (m)')\nax[0].set_xlabel(r'$\\rho$ ($\\Omega$m)')\nax[1].loglog(rhoa, ab2, 'rx-', label='data') # sounding curve\nax[1].loglog(inv.response(), ab2, 'b-', label='response')\nax[1].set_ylim((max(ab2), min(ab2))) # downwards according to penetration\nax[1].grid(True, which='both')\nax[1].set_xlabel(r'$\\rho_a$ ($\\Omega$m)')\nax[1].set_ylabel('AB/2 (m)')\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 }