302 lines
7.2 KiB
Plaintext
302 lines
7.2 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Checkout www.pygimli.org for more examples"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"# VES inversion for a blocky model\n",
|
|
"\n",
|
|
"This tutorial shows how an built-in forward operator is used for inversion.\n",
|
|
"A 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,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import pygimli as pg\n",
|
|
"from pygimli.physics import VESManager"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"some definitions before (model, data and error)\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"#ab2 = np.logspace(-0.5, 2.5, 40) # AB/2 distance (current electrodes)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"define a synthetic model and do a forward simulatin including noise\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"#synres = [100., 500., 30.] # synthetic resistivity\n",
|
|
"#synthk = [0.5, 3.5] # 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,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"#synthModel = synthk + synres # concatenate thickness and resistivity\n",
|
|
"#ves = VESManager()\n",
|
|
"#rhoa, err = ves.simulate(synthModel, ab2=ab2, mn2=0.5,\n",
|
|
"# noiseLevel=0.03, seed=1337)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# use data from Bausenberg\n",
|
|
"ab2 = [1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.5, 10.0, 12.5, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 75.0, 100.]\n",
|
|
"rhoa = [64.0, 70., 77., 89., 103., 115., 152., 179., 206., 244., 305., 360., 410., 489., 540., 588., 672., 605., 535., 433., 283.]\n",
|
|
"err = np.ones(len(rhoa))*0.02\n",
|
|
"err[-4:] = 0.05\n",
|
|
"mn2 = np.ones(len(rhoa))*0.5"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"print(ab2)\n",
|
|
"print(rhoa)\n",
|
|
"print(err)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"nlay = 3\n",
|
|
"ves = VESManager()\n",
|
|
"ves.invert(rhoa, err, ab2=ab2, mn2=mn2,\n",
|
|
" nLayers=nlay, lam=1000, lambdaFactor=0.8, cTpye=2, verbose=False)"
|
|
]
|
|
},
|
|
{
|
|
"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,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, ax = plt.subplots(ncols=2, figsize=(8, 6)) # two-column figure\n",
|
|
"#ves.showModel(synthModel, ax=ax[0], label=\"synth\", plot=\"semilogy\", zmax=20)\n",
|
|
"ves.showModel(ves.model, ax=ax[0], label=\"model\", zmax=50, color=\"C1\")\n",
|
|
"ves.showData(rhoa, ax=ax[1], label=\"data\", color=\"C0\", marker=\"o\")\n",
|
|
"out = ves.showData(ves.inv.response, ax=ax[1], label=\"model\", color=\"C1\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We are interested in the model uncertaincies and through model covariance\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from pygimli.frameworks.resolution import modelCovariance\n",
|
|
"var, MCM = modelCovariance(ves.inv)\n",
|
|
"pg.info(var)\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"im = ax.imshow(MCM, vmin=-1, vmax=1, cmap=\"bwr\")\n",
|
|
"plt.colorbar(im)\n",
|
|
"labels = [rf'$d_{i+1}$' for i in range(nlay-1)] + \\\n",
|
|
" [rf'$\\rho_{i+1}$' for i in range(nlay)]\n",
|
|
"plt.xticks(np.arange(nlay*2-1), labels)\n",
|
|
"_ = plt.yticks(np.arange(nlay*2-1), labels)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The model covariance matrix delivers variances and a scaled (dimensionless)\n",
|
|
"correlation matrix. The latter show the interdependency of the parameters\n",
|
|
"among each other. The first and last resistivity is best resolved, also the\n",
|
|
"first layer thickness. The remaining resistivities and thicknesses are highly\n",
|
|
"correlated. The variances can be used as error bars in the model plot.\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"thk = ves.model[:nlay-1]\n",
|
|
"res = ves.model[nlay-1:]\n",
|
|
"z = np.cumsum(thk)\n",
|
|
"mid = np.hstack([z - thk/2, z[-1]*1.1])\n",
|
|
"resmean = np.sqrt(res[:-1]*res[1:])\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"#ves.showModel(synthModel, ax=ax, label=\"synth\", plot=\"semilogy\", zmax=20)\n",
|
|
"ves.showModel(ves.model, ax=ax, label=\"model\", zmax=50)\n",
|
|
"ax.errorbar(res, mid, marker=\"*\", ls=\"None\", xerr=res*var[nlay-1:])\n",
|
|
"ax.errorbar(resmean, z, marker=\"*\", ls=\"None\", yerr=thk*var[:nlay-1])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"print(res)\n",
|
|
"print(res*var[nlay-1:])\n",
|
|
"print(thk)\n",
|
|
"print(thk*var[:nlay-1])"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "pyGIMLI",
|
|
"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.9"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|