Issue to invert real traveltime data

[OlgaPuntous]

Problem description

Hello,

I am trying to make cross hole travel time tomography Pygimli like in example. In example programm invert synthetic data, produces by function simulate data

tt.simulate(mesh=mesh_fwd, scheme=scheme, slowness=1. / model, noiseLevel=0.001, noiseAbs=1e-5), where is 0.1% relative and 10 microseconds of absolute noise.

In my case I would like to invert real traveltime data. I noticed that invertion result depend a lot from noise level input.... Which noise level I should introduce in my real data to get moreless stable result?

Thanks a lot for your advice.

Your environment

Windows

Python version: 3.6 pyGIMLi version: 2019 Way of installation: e.g. Conda package

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

import pygimli as pg
import pygimli.meshtools as mt

from pygimli.physics.traveltime import Refraction

from datetime import datetime
start_time = datetime.now()

mpl.rcParams['image.cmap'] = 'inferno_r'

################################################################################
# Next, we build the crosshole acquisition geometry with two shallow boreholes.

# Acquisition parameters
bh_spacing = 38.7
bh_length = 130

ssensor_spacing = 20
rsensor_spacing = 6.8

world = mt.createRectangle(start=[0, -80], end=[bh_spacing, -220],marker=0)
sources = -np.arange(110, 170 + ssensor_spacing, ssensor_spacing)
receivers=-np.arange(88, 196.9, rsensor_spacing)

sensors = np.zeros((len(sources) +len(receivers) , 2))  # two boreholes
sensors[len(sources):, 0] = bh_spacing  # x
sensors[:len(sources), 1] = np.hstack(sources)
sensors[len(sources):, 1] = np.hstack (receivers)

################################################################################
# Traveltime calculations work on unstructured meshes and structured grids. We
# demonstrate this here by simulating the synthetic data on an unstructured mesh
# and inverting it on a simple structured grid.

# Create forward model and mesh

c0 = mt.createRectangle(start=[0, -80], end=[bh_spacing, -102], marker=1)
c1 = mt.createRectangle(start=[0, -102], end=[bh_spacing, -136],marker=2)
c3= mt.createRectangle(start=[30, -80], end=[bh_spacing, -102],marker=3)
c4= mt.createRectangle(start=[30, -102], end=[bh_spacing, -136],marker=3)
c5= mt.createRectangle(start=[30, -136], end=[bh_spacing, -210],marker=3)

geom = mt.mergePLC([world,c0,c1,c3,c4,c5])

for sen in sensors:
    geom.createNode(sen)

mesh_fwd = mt.createMesh(geom, quality=34, area=.25)
model = np.array([15000., 40000, 26000,34000,34000,34000])[mesh_fwd.cellMarkers()]
pg.show(mesh_fwd, model, label="Velocity (m/s)", nLevs=2, logScale=False)

################################################################################
# Create inversion mesh
refinement = 0.25
x = np.arange(0, bh_spacing + refinement, rsensor_spacing * refinement)
y = -np.arange(80, bh_length + 100, rsensor_spacing* refinement)
mesh = pg.createMesh2D(x, y)
model1= np.array([15000.,40000,24000,34000,34000,34000])[mesh.cellMarkers()]
ax, _ = pg.show(mesh, hold=True)
ax.plot(sensors[:, 0], sensors[:, 1], "ro")

################################################################################
# Next, we create an empty DataContainer and fill it with sensor positions and
# all possible shot-recevier pairs for the two-borehole scenario using the
# product function in the itertools module (Python standard library).

from itertools import product
snumbers = np.arange(len(sources))
rnumbers= np.arange(len(receivers))
rays = list(product(snumbers,rnumbers+len(snumbers)))

# Empty container
scheme = pg.DataContainer()

# Add sensors
for sen in sensors:
    scheme.createSensor(sen)

# Add measurements
rays = np.array(rays)
scheme.resize(len(rays))
scheme.add("s", rays[:, 0])
scheme.add("g", rays[:, 1])
scheme.add("valid", np.ones(len(rays)))
scheme.registerSensorIndex("s")
scheme.registerSensorIndex("g")

################################################################################
# The forward simulation is performed with a few lines of code. We initialize an
# instance of the Refraction manager and call its `simulate` function with the
# mesh, the scheme and the slowness model (1 / velocity). We also add 0.1%
# relative and 10 microseconds of absolute noise.
#
# Secondary nodes allow for more accurate forward simulations. Check out the
# paper by `Giroux & Larouche (2013)
# <https://doi.org/10.1016/j.cageo.2012.12.005>`_ to learn more about it.

tt = Refraction()
#mesh_fwd.createSecondaryNodes(5)
#data = tt.simulate(mesh=mesh_fwd, scheme=scheme, slowness=1. / model, noiseLevel=0.001, noiseAbs=1e-5)
#data.save('mydataSM.dat')

data = tt.loadData('mydataVSP_M10.dat') 
mesh_fwd.createSecondaryNodes(5)
data = pg.DataContainer('mydataVSP_M10.dat', "s g")

ttinv = Refraction()
ttinv.setData(data)  # Set previously simulated data
ttinv.setMesh(mesh, secNodes=5)
invmodel = ttinv.invert(lam=1200, zWeight=5.0, startModel=1. / model1, startModelIsReference= True, )
print("chi^2 = %.2f" % ttinv.inv.getChi2())  # Look at the data fit

################################################################################
# Finally, we visualize the true model and the inversion result next to each
# other.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 7), sharex=True, sharey=True)
ax1.set_title("True model")
ax2.set_title("Inversion result")

pg.show(mesh_fwd, model, ax=ax1, showMesh=True, label="Velocity (m/s)",
        logScale=False, nLevs=1)

for ax in (ax1, ax2):
    ax.plot(sensors[:, 0], sensors[:, 1], "wo")

ttinv.showResult(ax=ax2,nLevs=2)
ttinv.showRayPaths(ax=ax2, color="0.8", alpha=0.5, nLevs=3)
fig.tight_layout()

end_time = datetime.now()
print('Duration: {}'.format(end_time - start_time))

################################################################################
# Note how the rays are attracted by the high velocity anomaly while
# circumventing the low velocity region. This is also reflected in the coverage,
# which can be visualized as follows:

#fig, ax = plt.subplots()
#ttinv.showCoverage(ax=ax, cMap="Greens")
#ttinv.showRayPaths(ax=ax, color="k", alpha=0.3)
#ax.plot(sensors[:, 0], sensors[:, 1], "ko")

################################################################################
# White regions indicate the model null space, i.e. cells that are not traversed
# by any ray.

Expected behavior

Actual behavior

Paste output here.
``
> If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

Estimating an error model (noise level) for real data is an essential step as it defines the regularization strength to be chosen and therefore the smoothness of the solution and the resolution properties. Looking at the data is therefore essential, but often you have rules of thumb from experience, by picking first arrivals or from considerations regarding frequency content and fresnel volumes. In Refraction seismics you often have an impression of the noise from how oscillating the traveltime curves are, in crosshole tomography it is a bit more complicated. What you can do is to look at reciprocal data, i.e. differences between traveltimes with interchanged shot/receiver positions, if available. At any rate it is nothing we can answer here as it is specific to your data, but any hints on estimating noise from others are welcome. In general, I can only recommend starting with a conservative guess and become more and more optimistic until artifacts show up in the model, then you're a bit too far.

Note that in the above code there is some inconsistency:

noiseLevel=0.001, noiseAbs=1e-5), where is 0.1% relative and 10 microseconds of absolute noise.

The relative noise level is only for stability and not important here.

[OlgaPuntous]

Thanks a lot Thomas for such a quick and clear answer.

[layufuji]

how to display data and response values ​​in 2D ERT modelling and inversion?

[halbmy]

how to display data and response values ​​in 2D ERT modelling and inversion?

This question does not seem to fit here. There are examples of 2D ERT including data and response at pygimli.org