{ "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# Regularization - concepts explained\n\nIn geophysical inversion, we minimize the data objective functional as\nthe L2 norm of the misfit between data $d$ and the forward\nresponse $f$ of the model $m$, weighted by the data error\n$\\epsilon$:\n\n\\begin{align}\\Phi_d = \\sum\\limits_i^N \\left(\\frac{d_i-f_i(m)}{\\epsilon_i}\\right)^2=\\|W_d(d-f(m))\\|^2\\end{align}\n\nAs this minimization problem is non-unique and ill-posed, we introduce a\nregularization term $\\Phi$, weighted by a regularization parameter\n$\\lambda$:\n\n\\begin{align}\\Phi = \\Phi_d + \\lambda \\Phi_m\\end{align}\n\nThe regularization strength $\\lambda$ should be chosen so that the\ndata are fitted within noise, i.e. $\\chi^2=\\Phi_d/N=1$.\n\nIn the term $\\Phi-m$ we put our expectations to the model, e.g. to\nbe close to any prior model. In many cases we do not have much\ninformation and aim for the smoothest model that is able to fit our\ndata. We decribe it by the operator $W_m$:\n\n\\begin{align}\\Phi_m=\\|W_m (m-m_{ref})\\|^2\\end{align}\n\nThe regularization operator is defined by some constraint operator\n$C$ weighted by some weighting function $w$ so that\n$W_m=\\mbox{diag}(w) C$. The operator $C$ can be a discrete\nsmoothness operator, or the identity to keep the model close to the\nreference model $m_{ref}$.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with importing the numpy, matplotlib and pygimli libraries\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\nimport pygimli as pg\nimport pygimli.meshtools as mt\nfrom pygimli.math.matrix import GeostatisticConstraintsMatrix\nfrom pygimli.core.math import symlog" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Regularization drives the model where the data are too weak to constrain\nthe model. In order to explain different kinds of regularization (also\ncalled constraints), we use a very simple mapping forward operator: The\nvalues at certain positions are picked.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pygimli.frameworks import PriorModelling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implementation 1. determine the indices where the cells are\n\n::\n\n ind = [mesh.findCell(po).id() for po in pos]\n\n2. forward response: take the model at indices\n\n::\n\n response = model[ind]\n\n3. Jacobian matrix\n\n::\n\n J = pg.SparseMapMatrix()\n J.resize(len(ind), mesh.cellCount())\n for i, n in enumerate(self.ind):\n self.J.setVal(i, n, 1.0)\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We exemplify this on behalf of a simple triangular mesh in a rectangular\ndomain.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rect = mt.createRectangle(start=[0, -10], end=[10, 0])\nmesh = mt.createMesh(rect, quality=34.5, area=0.3)\nprint(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define two positions where we associate two arbitrary values.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pos = [[3, -3], [7, -7]]\nvals = np.array([20., 15.])\nfop = PriorModelling(mesh, pos)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We set up an inversion instance with the forward operator and prepare\nthe keywords for running the inversion always the same way: - the data\nvector - the error vector (as relative error) - a starting model value\n(could also be vector)\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv = pg.Inversion(fop=fop, verbose=False)\ninvkw = dict(dataVals=vals, errorVals=np.ones_like(vals)*0.03, startModel=10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classical smoothness constraints\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(cType=1) # the default\nresult = inv.run(**invkw)\npg.show(mesh, result);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will have a closer look at the regularization matrix $C$.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "C = fop.constraints()\nprint(C.rows(), C.cols(), mesh)\nax, _ = pg.show(fop.constraints(), markersize=1)\n\nrow = C.row(111)\nnz = np.nonzero(row)\nprint(nz, row[nz])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does that change the regularization matrix $C$?\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(cType=1, zWeight=0.2) # the default\nresult = inv.run(**invkw)\npg.show(mesh, result)\n\nRM = fop.regionManager()\ncw = RM.constraintWeights()\nprint(min(cw), max(cw))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we try some other regularization options.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(cType=0) # damping difference to starting model\nresult = inv.run(**invkw)\nax, _ = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obviously, the damping keeps the model small ($\\log 1=0$) as the\nstarting model is NOT a reference model by default. We will enable this\nby specifying the isReference switch.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "invkw[\"isReference\"] = True\nresult = inv.run(**invkw)\nax, cb = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "cType=10 means a mix between 1st order smoothness (1) and damping (0)\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(cType=10) # mix of 1st order smoothing and damping\nresult = inv.run(**invkw)\nax, _ = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the matrix both contributions are under each other\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "C = fop.constraints()\nprint(C.rows(), C.cols())\nprint(mesh)\nax, _ = pg.show(fop.constraints(), markersize=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that we have the first order smoothness and the identity matrix\nbelow each other. We can also use a second-order (-1 2 -1) smoothness\noperator by cType=2.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(cType=2) # 2nd order smoothing\nresult = inv.run(**invkw)\nax, _ = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have a closer look at the constraints matrix\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "C = fop.constraints()\nprint(C.rows(), C.cols(), mesh)\nax, _ = pg.show(C, markersize=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like a Laplace operator and seems to have a wider range\ncompared to first-order smoothness.\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geostatistical regularization\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The idea is that not only neighbors are correlated to each other but to\nhave a wider correlation by using an operator\n\n\\begin{align}\\textbf{C}_{\\text{M},ij}=\\sigma^{2}\\exp{\\left(\n -\\sqrt{\n \\left(\\frac{\\textbf{H}^x_{ij}}{I_{x}}\\right)^{2}+\n \\left(\\frac{\\textbf{H}^y_{ij}}{I_{y}}\\right)^{2}\n }\\right)}.\\end{align}\n\nMore details can be found in\nhttps://www.pygimli.org/_tutorials_auto/3_inversion/plot_6-geostatConstraints.html\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We generate such a matrix and multiply it with a zero vector of just one 1.\nFor displaying the wide range of magnitudes we use the symlog function\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "C = GeostatisticConstraintsMatrix(mesh=mesh, I=[8, 4], dip=-20)\nprint(C)\n\nvec = pg.Vector(mesh.cellCount())\nvec[mesh.findCell([5, -5]).id()] = 1.0\nax, _ = pg.show(mesh, symlog(C*vec, 1e-2), cMin=-2, cMax=2, cMap=\"bwr\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For comparison, we use a much finer mesh and compute the same matrix\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fineMesh = mt.createMesh(rect, area=0.03)\nCfine = GeostatisticConstraintsMatrix(mesh=fineMesh, I=[8, 4], dip=-20)\nvec = pg.Vector(fineMesh.cellCount())\nvec[fineMesh.findCell([5, -5]).id()] = 1.0\nax, _ = pg.show(fineMesh, symlog(Cfine*vec, 1e-2), cMin=-1, cMax=1, cMap=\"bwr\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Application\nWe can pass the correlation length directly to the inversion instance\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(correlationLengths=[2, 2, 2])\nresult = inv.run(**invkw)\nax, cb = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This look structurally similar to the second-order smoothness, but can\ndrive values outside the expected range in regions of no data coverage.\nWe change the correlation lengths and the dip to be inclining\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inv.setRegularization(correlationLengths=[2, 0.5, 2], dip=-20)\nresult = inv.run(**invkw)\nax, cb = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now add many more points.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N = 30\nx = np.random.rand(N) * 10\ny = -np.random.rand(N) * 10\nv = np.random.rand(N) * 10 + 10\nplt.plot(x, y, \"*\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and repeat the above computations\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fop = PriorModelling(mesh, zip(x, y))\ninv = pg.Inversion(fop=fop, verbose=True)\ninv.setRegularization(correlationLengths=[4, 4])\nresult = inv.run(v, np.ones_like(v)*0.03, startModel=10)\nax, cb = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comparing the data with the model response is always a good idea.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.plot(v, inv.response, \"*\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Individual regularization operators\n\nSay you want to combine geostatistic operators with a damping, you can\ncreate a block matrix pasting the matric vertically.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "C = pg.matrix.BlockMatrix()\nG = pg.matrix.GeostatisticConstraintsMatrix(mesh=mesh, I=[2, 0.5], dip=-20)\nI = pg.matrix.IdentityMatrix(mesh.cellCount(), val=0.1)\nC.addMatrix(G, 0, 0)\nC.addMatrix(I, mesh.cellCount(), 0)\nax, _ = pg.show(C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that in pg.matrix you find a lot of matrices and matrix generators.\n\nWe set this matrix directly and do the inversion.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fop.setConstraints(C)\nresult = inv.run(v, np.ones_like(v)*0.03, startModel=10, isReference=True)\nax, cb = pg.show(mesh, result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you are using a method manager, you access the inversion instance by\nmgr.inv and the forward operator by mgr.fop.\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

#### Note

Take-away messages\n\n - regularization drives the model where data are weak\n - think and play with your assumptions to the model\n - there are several predefined options\n - geostatistical regularization can be superior, because:\n - it is mesh-independent\n - it better fills the data gaps (e.g.\u00a03D inversion of 2D profiles)

\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.8.10" } }, "nbformat": 4, "nbformat_minor": 0 }