Geoelektrik_Kolleg.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bedienungsanleitung:** Alle Zellen mit der Markierung `In [..]:` enthalten Python-Code. Zum Ausführen des Codes nutzen Sie die Tastenkombination `Shift-Return`. Die jeweils aktive ausführbare Zelle erkennen Sie am Rahmen mit dem grünen vertikalen Balken sowie der Markierung `In [..]:` am linken Rand.\n",
    "\n",
    "Bei der ersten Ausführung der Zelle `In [1]:` kann es zu einer Warnmeldung kommen, welche ignoriert werden kann."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Aufgabenstellung:** Arbeiten Sie das Notebook ab. In der Tabelle vor Code-Zelle 2 `In [2]:` finden Sie die Daten für eine geoelektrische Widerstandstiefensondierung $\\rho_s(AB/2)$. Diese sind in Code-Zelle 2 `In [2]:` bereits in die Python-Arrays `ab2` und `rhoa` eingetragen worden.\n",
    "\n",
    "Damit bleibt für Sie die Aufgabe, spezifische Widerstände $\\rho_i$, $i=1,2,3$ und Mächtigkeiten $h_i$, $i=1,2$ eines horizontalgeschichteten Dreischichtfalls zu bestimmen. Mit vorgegebenen Startwerten erzeugen Sie zunächst eine synthetische Sondierungskurve. Durch gezielte Veränderungen der Werte von $\\rho_i$ und $h_i$ versuchen Sie dann, die gemessene und synthetische Sondierungskurve schrittweise in bestmögliche Übereinstimmung zu bringen.\n",
    "\n",
    "Im Abschluss können Sie eine automatische Anpassung (*geophysikalische Dateninversion*) durchführen, um optimale Werte für $\\rho$ und $h$ zu erhalten.\n",
    "\n",
    "Für das Protokoll verwerten Sie die erzeugten Abbildungen sowie die Zahlenwerte für $\\rho_i$ und $h_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Auswertung einer Widerstandstiefensondierung\n",
    "## Grundlagen der Methode\n",
    "\n",
    "Bei einer Widerstandstiefensondierung werden scheinbare spezifische elektrische Widerstände $\\rho_s$ aufgezeichnet.\n",
    "Dabei wird der Abstand der Stromelektroden A ud B unter Beibehaltung des Mittelpunktes der Anordnung schrittweise vergrößert.\n",
    "Ist der Abstand zwischen den Potentialsonden M und N immer kleiner als AB/3, sprechen wir von einer *Schlumbergeranordnung*.\n",
    "\n",
    "Die Geoelektrik-Apparatur misst dabei elektrische Spannungen zwischen den Sonden M und N sowie den zwischen den Elektroden A und B fließenden elektrischen Strom.\n",
    "Daraus wird mit der Beziehung\n",
    "$$\n",
    "R = \\frac{U}{I} \\quad\\text{  in  } \\Omega\n",
    "$$\n",
    "zunächst der Ohmsche Widerstand ermittelt.\n",
    "Dieser wird mit dem Konfigurationsfaktor für die Schlumbergeranordnung\n",
    "$$\n",
    "k = \\frac{\\pi}{\\text{MN}}\\left( \\frac{\\text{AB}^2}{4} - \\frac{\\text{MN}^2}{4} \\right)\n",
    "$$\n",
    "multipliziert.\n",
    "\n",
    "Wir erhalten auf diese Weise den scheinbaren spezifischen elektrischen Widerstand\n",
    "$$\n",
    "\\rho_s = R \\cdot k \\quad\\text{  in  } \\Omega\\cdot m.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Das folgende Bild zeigt eine typische Sondierungskurve, zu deren Interpretation ein Dreischichtfall hinreichend ist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Sondierungskurve](soundingcurve.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Auswertung\n",
    "Das Ziel der Auswertung besteht in der Ermittlung der *spezifischen Widerstände* und *Mächtigkeiten* von Bodenschichten unter der Annahme einer näherungsweise horizontalen Lagerung.\n",
    "Die Messwerte werden zunächst graphisch als *Sondierungskurve* $\\rho_s = f(AB/2)$ dargestellt.\n",
    "Aus dem Kurvenverlauf schätzt man die minimale Anzahl von Schichten ab.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Für die Auswertung benutzen wir die Python-Bibliothek pygimli ([www.pygimli.org](http://www.pygimli.org)).\n",
    "\n",
    "Dazu importieren wir das Modul `functions`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functions import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pg.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Daten wurden für logarithmisch äquidistante Positionen der Stromelektroden aufgenommen.\n",
    "<!---\n",
    "\n",
    "Wir fassen die Messwerte in einer Tabelle zusammen:\n",
    "\n",
    "| AB/2 in m | $\\rho_s$ in $\\Omega\\cdot m$ |\n",
    "|-------------|-----------------------------|\n",
    "| 1.0 | 195.07 |\n",
    " |1.3   | 197.25 |\n",
    " |1.8  | 186.88 |\n",
    " |2.4   | 162.47 |\n",
    " |3.2   | 127.12 |\n",
    " |4.2   | 89.57 |\n",
    " |5.6   | 55.84 |\n",
    " |7.5   | 33.14 |\n",
    " |10.0   | 29.21 |\n",
    " |13.0   | 31.63 |\n",
    " |18.0   | 42.90 |\n",
    " |24.0   | 57.91 |\n",
    " |32.0   | 72.59 |\n",
    " |42.0   | 96.33 |\n",
    " |56.0   | 121.64 |\n",
    " |75.0   | 168.55 |\n",
    " |100.0   | 204.98 |\n",
    " \n",
    "Für alle Werte von AB/2 war der Abstand der Potentialsonden immer 0.6 m, d.h., MN/2 = 0.3 m.\n",
    "-->\n",
    "Wir fassen die Messwerte in den *Python*-Arrays `ab2`, `mn2` und `rhoa` zusammen. Die Richtung der Profilauslage wird in den Variablennamen durch Anhängen von ``_EW`` bzw. ``_NS`` gekennzeichnet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab2_EW = np.array([1.0, 1.3, 1.8, \n",
    "                2.4, 3.2, 4.2, \n",
    "                4.2, 5.6, 7.5, \n",
    "                10, 13, 18, \n",
    "                24, 32, 42])\n",
    "mn2_EW = np.array([0.3, 0.3, 0.3, \n",
    "                0.3, 0.3, 0.3, \n",
    "                1.0, 1.0, 1.0,\n",
    "                1.0, 1.0, 1.0,\n",
    "                1.0, 1.0, 1.0\n",
    "               ])\n",
    "rhoa_EW = np.array([138.18, 121.56, 106.22,\n",
    "                109.55, 114.79, 118.26,\n",
    "                112.39, 116.84, 125.58,\n",
    "                131.25, 129.04, 97.92,\n",
    "                90.23, 90.63, 76.16\n",
    "                ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab2_NS = np.array([1.00, 1.30, 1.80,\n",
    "        2.40, 3.20, 4.20,\n",
    "        5.60, 7.50, 10.00,\n",
    "        13.00, 18.00, 24.00, 32.00\n",
    "])\n",
    "mn2_NS = np.array([0.30, 0.30, 0.30,\n",
    "    0.30, 0.30, 0.30, 0.30,\n",
    "    1.00, 1.00, 1.00, 1.00,\n",
    "    1.00, 1.00\n",
    "])\n",
    "rhoa_NS = np.array([134.36, 117.62, 114.13,\n",
    "    110.14, 99.12, 98.05,\n",
    "    94.76, 106.13, 111.96,\n",
    "    120.07, 113.65, 92.58, 55.44\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisierung der Sondierungskurven\n",
    "\n",
    "Sondierungskurven werden in doppeltlogarithmischer Skalierung dargestellt. Auf der x-Achse werden die halben Elektrodenabstände (AB/2), auf der y-Achse die scheinbaren spezifischen Widerstände $\\rho_s$ abgetragen.\n",
    "\n",
    "Im Idealfall einer horizontalen, ebenen Schichtung müssen die Messwertkurven der Kreuzsondierungen übereinanderliegen.\n",
    "Abweichungen kommen zustande, wenn der Untergrund komplizierter aufgebaut ist (keine horizontale Lagerung, dreidimensionale Einlagerungen usw.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6,6))\n",
    "ax.loglog(ab2_EW, rhoa_EW, 'rx-', label='EW')\n",
    "ax.loglog(ab2_NS, rhoa_NS, 'bx-', label='NS')\n",
    "ax.set_xlim((1.0, 100.0))\n",
    "ax.set_ylim((10.0,1000.0))\n",
    "ax.set_xlabel('AB/2 in m')\n",
    "ax.set_ylabel(r'$\\rho_s$ in $\\Omega\\cdot m$')\n",
    "ax.legend(loc='best')\n",
    "ax.grid(True, which='both')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modellanpassung\n",
    "\n",
    "Wir versuchen jetzt durch Probieren, die gemessene Sondierungskurve mit einer anhand eines Modells berechneten Sondierungskurve in Übereinstimmung zu bringen.\n",
    "\n",
    "Dafür tragen wir in das Array `res` die Zahlenwerte für die spezifischen elektrischen Widerstände (in $\\Omega\\cdot m$) der drei Schichten beginnend an der Erdoberfläche (vom Hangenden zum Liegenden) ein:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = [180.0, 100.0, 300.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Zahlenwerte der Mächtigkeiten dieser Schichten (in m) fassen wir im Array `thk` zusammen. Beachten Sie, dass die Mächtigkeit der letzten Schicht (des Substratums) nach unten hin unbegrenzt ist und daher in `thk` weggelassen wird."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thk = [1.5, 5.5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Funktion `datenberechnen` berechnet die scheinbaren spezifischen Widerstände, die wir bei einer Messung über einem Dreischichtfall mit den spezifischen Widerständen und Mächtigkeiten für die vorgegebenen Stromelektrodenabstände erhalten würden:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab2 = ab2_EW\n",
    "mn2 = mn2_EW\n",
    "rhoa = rhoa_EW\n",
    "rhoanew = datenberechnen(ab2, mn2, res, thk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die so erhaltenen Ergebnisse (*Modellantwort*) `rhoanew` stellen wir gemeinsam mit den Messwerten (*Daten*) `rhoa` grafisch in Abhängigkeit vom Elektrodenabstand AB/2 `ab2` dar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datenvergleichen(rhoa, rhoanew, ab2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Da wir eine möglichst gute Übereinstimmung beider Kurven anstreben, ist es u.U. nötig, zur Definition von `res` und `thk` zurückzukehren, um die Rechnung mit veränderten Werten zu wiederholen.\n",
    "\n",
    "Notieren Sie die Werte von `res` und `thk`, wenn Sie zufrieden sind mit der Anpassung."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resbest = res\n",
    "thkbest = thk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatische Modellfindung\n",
    "Für eine automatische Ermittlung der spezifischen Widerstände `res` und Mächtigkeiten `thk` wird das Verfahren der geophysikalischen Dateninversion eingesetzt.\n",
    "\n",
    "Dabei wird ein aus den Messwerten erzeugtes Startmodell mit `nl` Schichten systematisch verändert, bis dessen Modellantwort mit den gemessenen Daten bis auf eine vorgegebene Abweichung `errPerc` übereinstimmt.\n",
    "Der Parameter `lam` steuert, wie groß die Sprünge zwischen den spezifischen Widerständen der einzelnen Schichten sein dürfen. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nl = 12\n",
    "lam = 50.0\n",
    "errPerc = 10.0\n",
    "resnew, thknew, rhoaresponse, relrms, chi2 = dateninversion(\n",
    "    ab2, mn2, rhoa, nl, lam, errPerc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotresults(resnew, thknew, ab2, rhoa, rhoaresponse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zusammenfassung der Resultate\n",
    "### Spezifische elektrische Widerstände"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Spezifische elektrische Widerstände in Ohm*m:\")\n",
    "for r in resnew:\n",
    "    print(f'{r:8.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mächtigkeiten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Schichtmächtigkeiten in m:\")\n",
    "for t in thknew:\n",
    "    print(f'{t:8.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bewertung der Anpassug\n",
    "Die Güte der Anpassung wird durch die Summe der Quadrate der Differenzen zwischen den gemessenen und synthetischen Werten der Sondierungskurve bestimmt.\n",
    "\n",
    "Wir unterscheiden zwischen *relativem RMS-Fehler* und $\\chi^2$-*Fehler*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Relativer RMS-Fehler = \", f'{relrms:.2f}')\n",
    "print(\"chi^2-Fehler         = \", f'{chi2:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weiterführende Analyse\n",
    "### Schichtäquivalenz\n",
    "Bei der manuellen bzw. automatischen Anpassung erhalten wir für die zweite Schicht folgende Werte für den spezifischen Widerstand und die Schichtmächtigkeit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Manuelle Anpassung:\")\n",
    "print(\"Spezifischer Widerstand der zweiten Schicht:\" + f'{resbest[1]:8.2f}' + \" Ohm*m\")\n",
    "print(\"Mächtigkeit der zweiten Schicht            :\" + f'{thkbest[1]:8.2f}' + \" m\")\n",
    "print(\"Automatische Anpassung:\")\n",
    "print(\"Spezifischer Widerstand der zweiten Schicht:\" + f'{resnew[1]:8.2f}' + \" Ohm*m\")\n",
    "print(\"Mächtigkeit der zweiten Schicht            :\" + f'{thknew[1]:8.2f}' + \" m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Für eine dünne Schicht mit Mächtigkeit $h$ und niedrigem Widerstand $\\rho$, die in einer Formation mit hohem spezifischen Widerstand eingebettet ist, gilt die Schichtäquivalenz. Dies bedeutet, dass die gesamte Sondierungskurve nahezu unverändert bleibt, sofern das Verhältnis\n",
    "$$\n",
    "S_i = \\frac{h_i}{\\rho_i} = const.\n",
    "$$\n",
    "gleich bleibt.\n",
    "Die Größe $S_i$ wird Längsleitfähigkeit der Schicht $i$ genannt.\n",
    "Ihre physikalische Einheit ist *Siemens*.\n",
    "Die Längsleitfähigkeit ist die einzige Information über die zweite Schicht, die bei der vorliegenden geoelektrischen Widerstandstiefensondierungen sicher bestimmt werden kann.\n",
    "Eine unabhängige Bestimmung von spezifischem Widerstand oder Mächtigkeit der zweiten Schicht ist dagegen nicht möglich.\n",
    "Man spricht vom *Äquivalenzprinzip der Geoelektrik*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Manuelle Anpassung:\")\n",
    "print(\"Längsleitfähigkeit der zweiten Schicht:\" + f'{thkbest[1]/resbest[1]:8.2f}' + \" S\")\n",
    "print(\"Automatische Anpassung:\")\n",
    "print(\"Längsleitfähigkeit der zweiten Schicht:\" + f'{thknew[1]/resnew[1]:8.2f}' + \" S\")"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python [conda env:pg]",
   "language": "python",
   "name": "conda-env-pg-py"
  },
  "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.10"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}