Adding data noise to measurements

Especially for synthetic modeling/inversion studies it is important to add data noise to the data.

While CRTomo contains built-in functionality to add data noise, we now recommend to add noise manually for better control.

Note that error parameters in the inversion context usually refer to standard deviations in the statistical sense. This implies that the actual realization of a given data noise component can be smaller or larger than the specified noise level. As such it may be possible/useful to sometimes reduce a given noise estimate below the actual noise level added to synthetic data.

For further reading, see:

https://en.wikipedia.org/wiki/Pseudorandom_number_generator

Imports

import numpy as np
import crtomo

Setup: Generate some synthetic data

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)
rmag = tdm.measurements()[:, 0]
rpha = tdm.measurements()[:, 1]

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 ::        5094\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          2                                                   \r Calculating Potentials : Wavenumber          3                                                   \r Calculating Potentials : Wavenumber          3                                                   \r Calculating Potentials : Wavenumber          4                                                   \r Calculating Potentials : Wavenumber          6                                                   \r Calculating Potentials : Wavenumber          6                                                   \r Calculating Potentials : Wavenumber          7                                                   \r Calculating Potentials : Wavenumber          8                                                   \r Calculating Potentials : Wavenumber          9                                                    done, now processing\n solution time  0d/  0h/  0m/  0s/  30ms\n\n Modelling completedSTOP 0\n'
reading voltages

Generate data noise

For synthetic studies a good starting point is that the structure of the actual noise should be the same as used in the inversion error model. As such we use a linear model for the magnitude noise components, and an absolute standard deviation for the phase values.

# Important: ALWAYS initialize the random number generator using a seed!
np.random.seed(2048)

# absolute component in [Ohm ]
noise_level_rmag_absolute = 0.01
# relative component [0, 1]
noise_level_rmag_relative = 0.05

noise_rmag = np.random.normal(
    loc=0,
    scale=rmag * noise_level_rmag_relative + noise_level_rmag_absolute
)

rmag_with_noise = rmag + noise_rmag

# 0.5 mrad absolute noise level
noise_level_phases = 0.5

noise_rpha = np.random.normal(
    loc=0,
    scale=noise_level_phases
)
rpha_with_noise = rpha + noise_rpha

# register the noise-added data as new measurements and mark them for use in a
# subsequent inversion
tdm.register_measurements(rmag_with_noise, rpha_with_noise)

Remove physically implausible negative magnitude values

indices = np.where(rmag_with_noise <= 0)[0]
tdm.configs.delete_data_points(indices)

Deleting configurations:
[]