Data error handling in CRTomo

CRTomo implements a last-square inversion with regularization based on the Thikonov appraoch. Herey a cost function is minimized, consisting of a data misfit \(\Psi_d\) and a model misfit \(\Psi_m\). Both terms are weighted against each other using a regularization parameter $lambda$.

\[\begin{split}\Psi &= \Psi_d + \lambda \cdot \Psi_m \\ \Psi_d &= \sqrt{\sum_i \left(\frac{d_i - f_i(m)}{\epsilon_i}\right)^2}\\ \Psi_m &= \left| \left| W_m m \right| \right|^2_2\end{split}\]

hereby \(W_m\) is the regularization matrix, which by default is a first-order approximation of the spatial gradient.

The data misfit term weights the difference between data and response of the complex resistivity model, \(f(m)\) against the data estimate \(\epsilon\). A data point is thus “fitted” if the forward response can reproduce it within the bounds of the error estimate. Numbers below 1 indicate that the model describes the data better than can be expected by the data error estimate, possibly indicating an overfitting in which noise components are treated as data, most often resulting in artifacts (artificial structures that are introduced by the inversion in order to explain the noise).

In principle an individual data error can be assigned to each data point. However, since data errors are not easily determined in geoelectrical applications, it is common to employ various data error models that exploit empirically observed relationships between the data errors and the measurements themselves.

Setup:

import crtomo

mesh = crtomo.crt_grid.create_surface_grid(nr_electrodes=10, spacing=1)
tdm = crtomo.tdMan(grid=mesh)

tdm.add_homogeneous_model(100, 0)
tdm.configs.gen_dipole_dipole(skipc=0)
tdm.measurements()

This grid was sorted using CutMcK. The nodes were resorted!
Triangular grid found
attempting modeling
b' ######### CMod ############\nLicence:\nCopyright \xc2\xa9 1990-2020 Andreas Kemna <kemna@geo.uni-bonn.de>\nCopyright \xc2\xa9 2008-2020 CRTomo development team (see AUTHORS file)\nPermission is hereby granted, free of charge, to any person obtaining a copy of\nthis software and associated documentation files (the \xe2\x80\x9cSoftware\xe2\x80\x9d), to deal in\nthe Software without restriction, including without limitation the rights to\nuse, copy, modify, merge, publish, distribute, sublicense, and/or sell copies\nof the Software, and to permit persons to whom the Software is furnished to do\nso, subject to the following conditions:\n\nThe above copyright notice and this permission notice shall be included in all\ncopies or substantial portions of the Software.\n\nTHE SOFTWARE IS PROVIDED \xe2\x80\x9cAS IS\xe2\x80\x9d, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\nIMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\nFITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\nAUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\nLIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\nOUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE\nSOFTWARE.\n \n \n CRMod Process_ID ::        5073\n OpenMP max threads:   4\n Git-Branch  master\n Commit-ID   cb89ccd2c4341b9b08bbe7bb08f5aea0d60916ee\n Created     Mon-Mar-25-14:34:56-2024\n Compiler    \n OS          GNU/Linux\n\n Reading Input-Files\n++ check 1\r++ check 2 done!\n\rGetting voltage      1\rGetting voltage      2\rGetting voltage      3\rGetting voltage      4\rGetting voltage      5\rGetting voltage      6\rGetting voltage      7\rGetting voltage      8\rGetting voltage      9\rGetting voltage     10\rGetting voltage     11\rGetting voltage     12\rGetting voltage     13\rGetting voltage     14\rGetting voltage     15\rGetting voltage     16\rGetting voltage     17\rGetting voltage     18\rGetting voltage     19\rGetting voltage     20\rGetting voltage     21\rGetting voltage     22\rGetting voltage     23\rGetting voltage     24\rGetting voltage     25\rGetting voltage     26\rGetting voltage     27\rGetting voltage     28 No electrode cap file\n \n Rescheduling..\n less nodes than wavenumbers\n OpenMP threads:   3(  4)\n\n\r Calculating Potentials : Wavenumber          1                                                   \r Calculating Potentials : Wavenumber          3                                                   \r Calculating Potentials : Wavenumber          3                                                   \r Calculating Potentials : Wavenumber          4                                                   \r Calculating Potentials : Wavenumber          5                                                   \r Calculating Potentials : Wavenumber          6                                                   \r Calculating Potentials : Wavenumber          7                                                   \r Calculating Potentials : Wavenumber          9                                                   \r Calculating Potentials : Wavenumber          9                                                    done, now processing\n solution time  0d/  0h/  0m/  0s/  30ms\n\n Modelling completedSTOP 0\n'
reading voltages

Magnitude error model

One common data error model for the resistivity/resistance parameters was formulated by LaBrecque et al 1996 and assumes a linear relationship between data error and the measured resistance:

\(\Delta R = a \cdot R + b\)

Hereby a is a relative component (in CRTomo this is a percentage value, \(a_{crtomo} = a / 100\)) and b and absolute value (in \(\Omega\)).

# set absolute error to 0.01 Ohm
tdm.crtomo_cfg['mag_abs'] = 0.01

# set relative error to 5 percent
tdm.crtomo_cfg['mag_rel'] = 5

Phase error model

Phase errors are implemented in CRTomo as:

\(\delta \phi = A1*abs(R)^B1 + A2*abs(pha) + p0\)

In practise, however, only the relative and absolute parameters A2 and p0 are used: 5% relative phase error

tdm.crtomo_cfg['pha_rel'] = 5

# 1 mrad absolute phase error
tdm.crtomo_cfg['pha_abs'] = 1