{
"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# Reciprocal data analysis of field ERT data\n\nAccording to the principle of reciprocity, the measured resistance of four-point\narrays does not change if the current bipole (A-B) is interchanged with the\npotential bipole (M-N). Already LaBrecque et al. (1996) suggested using the\ndifference between normal and reciprocal data for error estimation fitting\nthe variances through a model linear in the squared data.\nSlater et al. (2000) used the reciprocal error e=a+b*|R| to remove obvious\noutliers (>10%) and computed a and b forming an envelope.\n\nKoestel et al. (2008) suggested a statistical approach by binning the data\ninto intervals of similar magnitude and fitting a linear model through the\nstandard deviations of each bin over its mean value using an error model\n\n\\begin{align}\\delta R=a+b*|R|\\end{align}\n\nconsisting of a percentage error a and a resistance error b (in Ohm).\nThe relative error can be computed by\n\n\\begin{align}\\frac{\\delta R}{|R|} = b + a/|R|\\end{align}\n\nand reflects that the relative error (of either the voltage, the resistance or\napparent resistivity) has a minimum value of $a$ and increases with decreasing\nvoltage gain.\n\nThis procedure, originally designed for lab ERT with a small range of R values,\nwas used by Udphuay et al. (2011) for determining the error model for field data\nwith a much larger range for R values (section 4.4 and Fig. 11).\nFor the two data sets (called OP and RCP) they found rather low percentage errors\nof about 1% or below and ohmic errors of about 1.5m$\\Omega$ (unfortunately the\nmilli is missing in the paper) which is about 100\u00b5V considering a medium value of\nabout 60mA. This value is a good estimate that has been used as default value in\nthe BERT software, which was used to invert the data. To account for other error\nsources like positioning, topography, they increased the percentage error to 2%\nto determine the weighting in inversion as described by G\u00fcnther et al. (2006).\nThis procedure could also be made based on the voltage as this is the critical\nmeasure of signal strength.\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\nfrom pygimli.physics import ert"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After importing the necessary libraries, we load the RCP data from the paper\nof Udphuay et al. (2011).\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"data = pg.getExampleData('ert/reciprocal.ohm')\nprint(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It consist of more than 16000 data acquired on 14 ERT profiles with each 42\nelectrodes. The data container contains only resistances and no voltages or\napparent resistivity, which is computed through numerical geometric factors\nat a later processing stage. We show the electrode positions resembling\nFig. 8 of Udphuay et al. (2011), however rotated.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax, = plt.subplots()\nax.plot(pg.x(data), pg.y(data), \".\", label=\"electrodes\")\nax.set_xlabel(\"x (m)\")\nax.set_ylabel(\"y (m)\")\nax.grid()\nax.legend()\nax.set_aspect(1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before using the ready-made functions, we go through it step by step:\nFirst, we determine the indices for forward and backward measurements and\ngenerate a crossplot.\nIn our case, we have a number of 6143 pairs, large enough for any statistical\nanalyses, plus 4190 single data.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"iF, iB = ert.reciprocalIndices(data, True)\nprint(len(iF), data.size()-len(iF)-len(iB))\nfig, ax = plt.subplots()\nax.loglog(data['r'][iF], data['r'][iB], '.', markersize=2)\nax.set_xlabel(\"R normal (Ohmm)\")\nax.set_ylabel(\"R reciprocal (Ohmm)\")\nax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Except some single outliers, all data are aligned along the identity axis,\nbut clearly with a larger spread for lower resistances as a result of\nincreased (relative as indicated by the logarithmic axis) errors.\nWe look at reciprocal errors a bit more in detail:\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"R = data['r']\nRmean = abs(R[iF] + R[iB]) / 2\ndR = abs(R[iF] - R[iB])\nfig, ax = plt.subplots()\nax.loglog(Rmean, dR, '.', markersize=2)\nax.set_xlabel('R (Ohm)')\nax.set_ylabel('dR(Ohm)')\nax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This plot rather suggests a log-log dependence\n\\begin{align}\\log\\delta R = \\log a + b \\log R \\quad\\Rightarrow\\quad \\delta R = a R^b\\end{align}\n\nas it was suggested by Flores-Oroczo et al. (2018) for phase data.\n\nAlternatively, we can plot the relative error on the y axis.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dRbyR = dR / Rmean\nfig, ax = plt.subplots()\nax.loglog(Rmean, dRbyR, '.', markersize=2)\nax.set_xlabel('R (Ohm)')\nax.set_ylabel('dR / R')\nax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It shows that the majority of data pairs have errors of 1% or below.\nThe statistical analysis (binning, std and mean values) is done by a\ndedicated function `fitReciprocalModel` that can also plot it all.\nNote that the function also takes different signs (here all R are\npositive) and possible NaN values into account.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ab1, ax = ert.fitReciprocalErrorModel(data, show=True)\nprint(ab1)\nax.set_yscale(\"symlog\", linthresh=0.0001)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we end up with $b\\approx$ 1% and $a\\approx$ 0.8 mOhm.\nAlternatively, we can fit the relative errors by `rel=True`.\nThis yields values of about 2.5% and 0.2 mOhm, here the first\none is more realistic but the second appears rather low.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ab2, ax = ert.fitReciprocalErrorModel(data, rel=True, show=True)\nprint(ab2)\nax.set_yscale(\"symlog\", linthresh=0.001)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We plot the two models against each other.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"relerr1 = ab1[1] + ab1[0] / data['r']\nrelerr2 = ab2[0] + ab2[1] / data['r']\nfig, ax = plt.subplots()\nax.loglog(relerr1*100, relerr2*100, \".\", markersize=2)\nax.set_aspect(1.0)\nax.set_xlabel('model 1 error (%)')\nax.set_ylabel('model 2 error (%)')\nax.set_xlim(1, 200)\nax.set_ylim(1, 100)\nax.grid(which='both')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, the second error model (2.3-40% relative error) based on\nfitting relative values, is more compact and appears more realistic than the\nfirst error model (1-180%) based on fitting absolute values.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Further procedure\nA N/R pair of data cannot be independently fitted. Therefore, it is\nreasonable to remove all values above a certail threshold as suggested by\nSlater et al. (2000) who used a threshold of 10% (here `sum(dRbyR>0.1)`=223).\nFurthermore, it makes sense to keep only one quadrupole of each pair by taking\nthe mean resistance, possibly weighted by the current strength or stacking\nerror. This \"stacking\" can improve the data quality. Alternatively, one could\nuse some heuristics to decide for one of the data and to delete the other, but\nthis should rather be done before the N/R analysis.\nThe error of all data is then estimated by the used error model, optionally\nlarge individual values higher than the error model could be taken instead of\nthese. At any rate, I suggest removing data with either large reciprocal error\nor large error misfits above values in the range of 20%.\nAll this is implemented in the function `reciprocalProcessing`, which basically\naverages the N/R pairs and sets the property `data['err']` that is used in\ninversion to weight the data. One can pass a maximum reciprocal error\n(`maxrec`) to eliminate pairs with wrong reciprocity completely, and maximum\nerror (`maxerr`) for large error estimate that might destroy the statistics.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"out = ert.reciprocalProcessing(data, maxrec=0.2, maxerr=0.2)\nprint(out)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we end up in a data container of 9456 data without any pairs that\ncan go into inversion.\n\n## References\n* LaBrecque, D.L., Miletto, M., Daily, W., Ramirez, A. & Owen, E. (1996): The\neffects of noise on Occam\u2019s inversion of resistivity tomography data,\nGeophysics, 61, 538\u2013548.\n* Slater, L., A. M. Binley, W. Daily, and R. Johnson (2000): Cross-hole\nelectrical imaging of a controlled saline tracer injection, J. Appl. Geophys.,\n44, 85\u2013102, doi:10.1016/S0926-9851(00)00002-1.\n* G\u00fcnther, T., R\u00fccker, C. & Spitzer, K. (2006): Three-dimensional modeling and\ninversion of dc resistivity data incorporating topography \u2013 II: Inversion.\nGeophys. J. Int. 166, 506-517, doi:10.1111/j.1365-246X.2006.03011.x.\n* Koestel, J., Kemna, A., Javaux, M., Binley, A. & Vereecken, H. (2008):\nQuantitative imaging of solute transport in an unsaturated and undisturbed soil\nmonolith with 3-D ERT and TDR, Water Resour. Res., 44, W12411,\ndoi:10.1029/2007WR006755.\n* Udphuay, S., G\u00fcnther, T., Everett, M.E., Warden, R.R. & Briaud, J.-L. (2011):\nThree-dimensional resistivity tomography in extreme coastal terrain amidst dense\ncultural signals: application to cliff stability assessment at the historic D-Day\nsite. Geophys. J. Int. 185(1), 201-220, doi:10.1111/j.1365-246X.2010.04915.x.\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.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}