{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Checkout www.pygimli.org for more examples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Structural coupled cooperative inversion (SCCI)\n\nJoint inversion is an important method to improve resolution properties by\ncombining different methods. G\u00fcnther & R\u00fccker (2006) and G\u00fcnther et al.\n(2010) introduced a scheme later refined by Hellmann et al. (2017), Ronczka\net al. (2017), and Skibbe et al. (2018, 2021).\n\nWe use the model already used in the pyGIMLi paper (R\u00fccker et al., 2017) for\npetrophysical joint inversion.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We import the necessary libraries and the SCCI class\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport numpy as np\nimport pygimli as pg\nfrom pygimli import meshtools as mt\nfrom pygimli.physics import ert\nfrom pygimli.physics import traveltime as tt\nfrom pygimli.frameworks import SCCI\nfrom pygimli.viewer.mpl.meshview import drawCWeight\nfrom pygimli.physics.petro import transFwdArchieS as ArchieTrans\nfrom pygimli.physics.petro import transFwdWyllieS as WyllieTrans" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def createSynthModel(nSeg=32):\n \"\"\"Return the modelling mesh, the porosity distribution and the\n parametric mesh for inversion.\n \"\"\"\n # Create the synthetic model\n world = mt.createCircle(boundaryMarker=-1, nSegments=nSeg*2)\n tri = mt.createPolygon([[-0.8, -0], [-0.5, -0.7], [0.7, 0.5]],\n isClosed=True, area=0.0015)\n c1 = mt.createCircle(radius=0.2, pos=[-0.2, 0.5], segments=16,\n area=0.0025, marker=3)\n c2 = mt.createCircle(radius=0.2, pos=[0.32, -0.3], segments=16,\n area=0.0025, marker=3)\n\n poly = mt.mergePLC([world, tri, c1, c2])\n\n poly.addRegionMarker([0.0, 0, 0], 1, area=0.0015)\n poly.addRegionMarker([-0.9, 0, 0], 2, area=0.0015)\n\n c = mt.createCircle(radius=0.99, nSegments=nSeg, start=np.pi, end=np.pi*3)\n [poly.createNode(p.pos(), -99) for p in c.nodes()]\n poly.save(\"poly.bms\")\n mesh = pg.meshtools.createMesh(poly, q=34.4, smooth=[1, 10])\n mesh.scale(1.0/5.0)\n mesh.rotate([0., 0., 3.1415/3])\n mesh.rotate([0., 0., 3.1415])\n\n petro = pg.solver.parseArgToArray([[1, 0.9], [2, 0.6], [3, 0.3]],\n mesh.cellCount(), mesh)\n\n # Create the parametric mesh that only reflect the domain geometry\n world = mt.createCircle(boundaryMarker=-1, nSegments=nSeg*2, area=0.0051)\n paraMesh = pg.meshtools.createMesh(world, q=34.0, smooth=[1, 10])\n paraMesh.scale(1.0/5.0)\n\n return mesh, paraMesh, petro" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mMesh, pMesh, saturation = createSynthModel()\nrKW = dict(logScale=True, cMin=250, cMax=2500, cMap=\"Spectral_r\")\nvKW = dict(logScale=True, cMin=1000, cMax=2500, cMap=\"Spectral_r\")\nertTrans = ArchieTrans(rFluid=20, phi=0.3)\nres = ertTrans(saturation)\nttTrans = WyllieTrans(vm=4000, phi=0.3)\nvel = 1./ttTrans(saturation)\nsensors = mMesh.positions()[mMesh.findNodesIdxByMarker(-99)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pg.info(\"Simulate ERT\")\nERT = ert.ERTManager(verbose=False, sr=False)\nertScheme = ert.createERTData(sensors, schemeName='dd', closed=1)\nertData = ert.simulate(mMesh, scheme=ertScheme, res=res, noiseLevel=0.01)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pg.info(\"Simulate Traveltime\")\nTT = tt.TravelTimeManager(verbose=False)\nttScheme = tt.createRAData(sensors)\nttData = tt.simulate(mMesh, scheme=ttScheme, vel=vel, secNodes=5,\n noiseLevel=0.001, noiseAbs=2e-6)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ERT = ert.ERTManager(ertData, verbose=True, sr=False)\nERT.setMesh(pMesh)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "TT = tt.TravelTimeManager(ttData, verbose=True)\nTT.errIsAbsolute = True\nTT.setMesh(pMesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ERT.setRegularization(cType=1)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ERT.invert(zWeight=1, lam=100, maxIter=3)\n# ERT.showResult(**rKW)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "TT.invert(zWeight=1, lam=400, secNodes=5, startModel=1/2000, maxIter=3)\n# ax, cb = TT.showResult(**vKW)\n# ax.plot(pg.x(ttData), pg.y(ttData), \"ko\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "scci = SCCI([ERT, TT], names=[\"ERT\", \"TT\"])\nscci.a = 0.1\nscci.b = 0.02\nscci.c = 1\nscci.cmax = 100" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cw = scci.singleCWeights()\nprint(min(cw[0]), min(cw[1]), np.mean(cw[0]), np.mean(cw[1]))\nfig, ax = pg.plt.subplots(ncols=2)\naa, cb = ERT.showResult(ax=ax[0], **rKW)\ndrawCWeight(aa, ERT.paraDomain, cw[0])\naa, cb = TT.showResult(ax=ax[1], **vKW)\ndrawCWeight(aa, TT.paraDomain, cw[1])\nvPre = pg.Vector(TT.inv.model)\nrPre = pg.Vector(ERT.inv.model)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "TT.inv.lam = 1000\nERT.inv.lam = 200\nscci.run(maxIter=5) # save=True)\nprint(ERT.inv.chi2(), TT.inv.chi2())\nprint(min(ERT.inv.inv.cWeight()), min(TT.inv.inv.cWeight()))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ERT.inv.model = ERT.inv.inv.model()\nTT.inv.model = 1./TT.inv.inv.model()\nvPost = TT.inv.model\nrPost = ERT.inv.model\nprint(rPre[0], rPost[0], pg.math.rrms(rPre, rPost))\nprint(vPre[0], vPost[0], pg.math.rms(vPre, vPost))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = pg.plt.subplots(ncols=2)\naa, cb = ERT.showResult(ax=ax[0], **rKW)\naa.plot(pg.x(ttData), pg.y(ttData), \"ko\")\n# drawCWeight(aa, ERT.paraDomain, ERT.inv.inv.cWeight())\naa, cb = TT.showResult(ax=ax[1], **vKW)\naa.plot(pg.x(ttData), pg.y(ttData), \"ko\")\n# drawCWeight(aa, TT.paraDomain, TT.inv.inv.cWeight())" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# References\n#\n# G\u00fcnther, T., Dlugosch, R., Holland, R. & Yaramanci, U. (2010): Aquifer characterization using coupled inversion of DC/IP and MRS data on a hydrogeophysical test-site. - SAGEEP 23, 39 (2010); Keystone, CO.\n# G\u00fcnther, T. & R\u00fccker, C. (2006): A new joint inversion approach applied to the combined tomography of dc resistivity and seismic refraction data. - Ext. abstract, 19. EEGS annual meeting (SAGEEP), 02.-06.04.2006; Seattle, USA.\n# Hellman, K., Ronczka, M., G\u00fcnther, T., Wennermark, M., R\u00fccker, C. & Dahlin, T. (2017): Structurally coupled inversion of ERT and refraction seismic data combined with cluster-based model integration. Journal of Applied Geophysics 143, 169-181, doi:10.1016/j.jappgeo.2017.06.008.\n# Ronczka, M., Hellman, K., G\u00fcnther, T., Wisen, R., Dahlin, T. (2017): Electric resistivity and seismic refraction tomography, a challenging joint underwater survey at Asp\u00f6 hard rock laboratory. Solid Earth 8, 671-682. doi:10.5194/se-8-671-2017.\n# R\u00fccker, C., G\u00fcnther, T., Wagner, F.M. (2017): pyGIMLi: An open-source library for modelling and inversion in geophysics, Computers & Geosciences 109, 106-123, doi:10.1016/j.cageo.2017.07.011.\n# Skibbe, N., G\u00fcnther, T. & M\u00fcller-Petke, M. (2018): Structurally coupled cooperative inversion of magnetic resonance with resistivity soundings. Geophysics 83(6), JM51-JM63, doi:10.1190/geo2018- 0046.1.\n# Skibbe, N., G\u00fcnther, T. & M\u00fcller-Petke, M. (2021): Improved hydrogeophysical imaging by structural coupling of two-dimensional magnetic resonance and electrical resistivity tomography. Geophysics 86 (5), WB135-WB146, doi:10.1190/geo2020-0593.1." ] } ], "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.11.2" } }, "nbformat": 4, "nbformat_minor": 0 }