startModel for TravelTimeManager.invert()

[Olga5991]

Hello,
Could you please provide me with information what should be used as startModel in TravelTimeManager.invert()? In description it is said to be float (one float number?). I thought it as a numpy array with start velocity model, and as i see it should be constructed based on the mesh. I have preliminary information about velocity model and I want to use it as a start model for inversion. How do you think, is it possible to obtain better result through it?

[halbmy]

You can provide either a floating number for a homogeneous slowness or a self-defined vector. However, you can also work with the vTop and vBottom keywords to generate a gradient model as typical in refraction tomography, or (if both are constant) a homogeneous model as typical in crosshole tomography.

[anonymouscapibara]

I think I have got the same question of @Olga5991
I created a simple 4-layer model
left = 0; right = 66; depth = 25;
geometry= mt.createWorld(start = [left, 8], end = [right, -depth], layers=[5, 0, -19]);

the I created a mesh
meshVel = mt.createMesh(geometry, quality=33.5, area=2, smooth=True);

and finally populated the layers with seismic velocities (though I am afraid these velocities are not assigned to the cells of the mesh)
velMap = [[1, 200], [2, 600], [3, 1500], [4, 2000]];

I would like to set the above mesh, with the related seismic velocities, as starting model for the traveltime inversion
pygimli.physics.traveltime.TravelTimeManager.invert

I know how to set either a homogeneous or a gradient starting model. On the contrary, it is not clear to me how to "import" a user-defined starting mesh (see code above) and how to use a self-defined vector.
I also read that

The inversion will always start from the starting model taken from the forward operator. If you want to run the inversion from a specified prior model, e.g., from a other run, set this model as starting model to the FOP (fop.setStartModel). Any self.inv.setModel() settings will be overwritten.

but could not figure it out anyway.
Could you please give some hints?
Thank you very much.
Diego

[halbmy]

import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt

From the mesh

left = 0; right = 66; depth = 25;
geometry= mt.createWorld(start = [left, 8], end = [right, -depth], layers=[5, 0, -19])
meshVel = mt.createMesh(geometry, quality=33.5, area=2, smooth=True)

one can generate a velocity vector by mapping the map

velMap = [[1, 200], [2, 600], [3, 1500], [4, 2000]];
vel = pg.solver.parseMapToCellArray(velMap, meshVel)
pg.show(meshVel, vel)

that can be used as a starting model (slowness!) in inversion

mgr = tt.TravelTimeManager(data)
mgr.invert(mesh=meshVel, startModel=1/vel)

Or did I get it wrong?

[anonymouscapibara]

Thank you for your prompt reply.
Unfortunately your code does not run fine, at least on my PC.
I copy-and-pasted your very code above, and, before initializing the TravelTimeManager with the data, I loaded my picking data (also attached to this message)
data = tt.load('picking_NS_20240607.sgt', verbose=True)
print(data)
fig, ax = plt.subplots();
lines = tt.drawFirstPicks(ax, data);
data['err'] = data['t'] * 0.03 # Adding 3% picking error to data
print(data)

However, when I run the inversion, I get this error

So it seems there is a mismatch in the number of cells. And even if a constant model is created, there is a runtime error.

I thank you for any suggestions.
Diego
picking_NS_20240607.zip

[halbmy]

I think the problem are the different regions (layers) of which one is taken as background. You can either remove this prior to inversion (but after generating the starting velocityvector) by resetting the region markers:

meshVel.setCellMarkers(pg.Vector(meshVel.cellCount()))

(from v1.4.7 also by meshVel['marker']=0) OR setting all regions explicitly NOT to background)

mgr = tt.TravelTimeManager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)

[anonymouscapibara]

I tried both strategies.

OPTION 1 (Working - but found 1 region)
meshVel.setCellMarkers(pg.Vector(meshVel.cellCount()))
mgr = tt.TravelTimeManager(data)
mgr.invert(mesh=meshVel, startModel=1/vel, useGradient=False, verbose=True)

OPTION 2 (Not working - but found 4 regions)
mgr = tt.TravelTimeManager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)
mgr.invert(mesh=meshVel, startModel=1/vel, useGradient=False, verbose=True)

BTW, my pygimli is 1.4.6

[halbmy]

I could not read the whole error message (i.e. where exactly the error occurs).

[anonymouscapibara]

My apologies. Below I paste the complete error message.
Thank you

[halbmy]

The problem is in the regions of the mesh. Actually, it does not make sense to use the mesh from the synthetic modelling also in inversion (this is an inverse crime), because then you will of course reconstruct the layers as these are already decoupled. So better don't specify a mesh or use a newly created that does not already include your result.

Independent on this, it should work. The problem is most probably a bad default region setting (using one as background). Can you try

mgr = tt.TravelTimeManager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)
mgr.invert(startModel=1/vel)

[anonymouscapibara]

First, regarding the inverse crime, actually I am not inverting synthetic modelling data but just field data. So any mesh I give as input to the inversion should be fine. Or I did not get it?
Second, your last suggestion seems to work. Better, now there are no errors. But the result is unacceptable: something is very bad with the first layer, which is the region with the smaller marker (i.e., 1), that the software sets to background when inversion is run. Please see image below.

In the end my goal is very simple: run the inversion with a regularization (0th-order) given by the reference layered model that I created. I also tried to manually set 0-th order regularization with cType argument, but results are even worse.
Below my complete code ant attached the file with sources-geophones geometry and picked traveltimes.
Thank you for any hints.

import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt
import numpy as np
import matplotlib.pyplot as plt

left = 0;
right = 66;
depth = 25;
world = mt.createWorld(start = [left, 7.7], end = [right, -depth], layers=[5, 0, -19]);
print(world)
pg.show(world, boundaryMarkers=True);

meshVel = mt.createMesh(world, quality=33.5, area=2, smooth=True);
print(meshVel)
pg.show(meshVel, markers=True, showMesh=True);

velMap = [[1, 500], [2, 700], [3, 1600], [4, 2500]];
vel = pg.solver.parseMapToCellArray(velMap, meshVel);
pg.show(meshVel, vel, label='Velocity [m/s]');
print(vel)

data = tt.load('picking_NS_20240607.sgt', verbose=True)
print(data)
fig, ax = plt.subplots();
lines = tt.drawFirstPicks(ax, data);
# ADDING % PICKING ERRORS TO DATA
data['err'] = data['t'] * 0.03    # 3%

# DATA INVERSION
mgr = tt.TravelTimeManager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)
mgr.invert(startModel=1/vel, lam=50, verbose=True)
rrms = mgr.inv.relrms();
chi2stat = mgr.inv.chi2();

# PLOT INVERSION RESULTS
plt.rcParams['figure.figsize'] = [12, 6]
figInvResults = plt.figure(constrained_layout=True)
gs = figInvResults.add_gridspec(5, 5)

# Plot inverted model ------------------------------------------------------------
figInvResults_ax1 = figInvResults.add_subplot(gs[:-2, :])
figInvResults_ax1.set_title('INVERTED MODEL')
_ , cbar = mgr.showResult(ax=figInvResults_ax1, logScale=False, cMin=100, cMax=2000,
                            contour=True, cMap="jet", nLevs=6, orientation='vertical') #, coverage=mgr.standardizedCoverage())
                            # nLevs  --> number of colormap labels
# Plot rays
rays = mgr.drawRayPaths(ax=figInvResults_ax1, color="k", lw=0.3, alpha=0.5)
# mgr.coverage() yields the ray coverage in m and standardizedCoverage as 0/1
# Plot sensor positions
pg.viewer.mpl.drawSensors(figInvResults_ax1, data.sensors(), diam=.5, facecolor='red', edgecolor='black')
# Set axes limits
figInvResults_ax1.set_xlim(0, 66)
figInvResults_ax1.set_ylim(-18, 8)
figInvResults_ax1.set_xlabel('X [m]')
figInvResults_ax1.set_ylabel('Elevation a.s.l. [m]')
# Add RRMS text
figInvResults_ax1.text(55, -15, 'RRMS = ' + str(np.round(rrms, 2)) + '%', color="black", va="center", ha="center");
figInvResults_ax1.text(55, -16.5, r'$\chi^2$ = ' + str(np.round(chi2stat, 2)), color="black", va="center", ha="center");

# Plot goodness of fit (lines) ------------------------------------------------------------
figInvResults_ax2 = figInvResults.add_subplot(gs[3:, :-3])
figInvResults_ax2.set_title('GOODNESS OF FIT [1]')
mgr.showFit(ax=figInvResults_ax2, firstPicks=True)
figInvResults_ax2.set_xlabel('X [m]')
figInvResults_ax2.set_ylabel('Traveltimes [s]')

# Plot goodness of fit (points) ------------------------------------------------------------
figInvResults_ax3 = figInvResults.add_subplot(gs[3:, 3:])
figInvResults_ax3.set_title('GOODNESS OF FIT [2]')
reference= np.linspace(0,1.05*max(max(data["t"]),max(data["t"])))
plt.plot(data['t'], mgr.inv.response, 'ko-', linewidth = 0, markersize = 3);  #dati_calcolati
plt.plot(reference, reference, 'r:', linewidth = 0.5)
plt.axis([0,1.05*max(max(data["t"]),max(data["t"])),0,1.05*max(max(data["t"]),max(data["t"]))]);
plt.grid()
plt.ylabel('Computed traveltimes [s]');
plt.xlabel('Measured traveltimes [s]');
# Add R^2 text
ss_res = np.sum((data['t'] - mgr.inv.response)**2);
ss_tot = np.sum((data['t'] - np.mean(data['t']))**2);
Rsquared = 1 - (ss_res/ss_tot);
figInvResults_ax3.text(0.06, 0.02, '$R^2$ = ' + str(np.round(Rsquared, 3)), color="black", va="center", ha="center", size= '11');

picking_NS_20240607.zip

[halbmy]

Sorry, I didn't mean to affront you. Of course you can use any mesh.

As to the error, I think an absolute error (e.g. 0.001) makes more sense, otherwise the small offsets with small traveltimes will be overly weighted and produce artifacts like in your case.

Furthermore, the mesh might be a bit too coarse, so maybe decrease area or increase quality.

The apparent velocity (tt.showVA(data, pseudosection=True)

Already shows minimum values below air speed for small but not near-zero offsets, maybe a hint to a hidden layer.

data["err"] = 0.001
mgr = tt.Manager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)
mgr.invert(verbose=True)
mgr.showResult()

fits the data with chi-sqaure 1.1, i.e. RMS=1ms

Additionally specifying startModel=1/vel yields

a similar and equivalent (chi2=1.2) model.

[halbmy]

I put the whole notebook to https://github.com/gimli-org/notebooks/blob/main/issues/642/issue642.ipynb

[anonymouscapibara]

Dear Thomas,
absolutely no problem! In fact, I really thank you for your support.
Good point about the absolute error in place of the relative ones. I should have figured it out myself. I do apologize.

Regarding your first inversion approach (GRADIENT STARTING MODEL)

mgr = tt.Manager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)
mgr.invert(verbose=True)

I see the starting model is a gradient one (from 500 to 5000m/s) and, accordingly, one region is found. However, when you show the results, the model shows the layers that have been created in the meshVel mesh. So I am a bit puzzled about what the command mgr.setMesh(meshVel) really does. The command docstring is Set a mesh and distribute it to the forward operator. In the end, are the synthetic traveltimes computed considering the gradient or the layered starting model in the iterative (Gauss-Newton) inversion process?

In addition to what above, if I plot the results masked by coverage and with rays on top (see image below), I think I clearly see that there is something wrong with the rays (and with coverage accordingly), as they seem to start from z = 0. Or, better, there seems to be no rays in the topmost layer, where sources and geophones are placed. Or am I missing something?

Regarding your second inversion approach (LAYERED STARTING MODEL)

mgr.invert(startModel=1/vel, verbose=True)

It seems that the correct layered starting model is given as input, though there is no info printed about the different model regions. I think the starting model is correct because of the printed slowness range (see image below - velocity from 200 to 2000m/s)

But, again, there seems to be that strange ray behavior

[halbmy]

Obviously there is a mismatch between data and mesh.

You should add the data sensors to the world so that they appear as nodes. Formally, this is not required but it seems like wrong ray paths are constructed.

[anonymouscapibara]

Dear Thomas,
I tried to add sources' and geophones' positions as nodes, but seismic rays still have that odd behavior and they start at z = 0 (i.e., within the model).
Here is my code

left = 0; right = 66; depth = 25;
geometry= mt.createWorld(start = [left, 7.7], end = [right, -depth], layers=[5, 0, -19]);

nsensors = 24;
sources = np.ones((9,2))*7.7;
sources[:,0] = [0, 8, 16, 24, 32, 40, 48, 58, 66];
geophones = np.ones((nsensors,2))*7.7;
geophones[:,0] = np.linspace(10, 56, nsensors)
#
for pos in sources:
    geometry.createNode(pos)
for posGeo in geophones:
    geometry.createNode(posGeo)
#
ax, _ = pg.show(geometry)
ax.plot(sources[:,0],sources[:,1],'ro')
ax.plot(geophones[:,0],geophones[:,1],'bv')

meshVel = mt.createMesh(geometry, quality=34.3, area=1.5, smooth=True);
velMap = [[1, 200], [2, 600], [3, 1500], [4, 2000]];
vel = pg.solver.parseMapToCellArray(velMap, meshVel)
# vel = meshVel.populate("vel", velMap)                          ####### DOES NOT WORK!
meshVel.setCellMarkers(pg.Vector(meshVel.cellCount()))
pg.show(meshVel, vel)

data = tt.load("picking_NS_20240607.sgt")
data["err"] = 0.001
mgr = tt.Manager(data)
mgr.setMesh(meshVel)
mgr.inv.setRegularization('*', background=False)
mgr.invert(verbose=True)
#
ax, cbar = mgr.showResult(coverage=mgr.standardizedCoverage())
rays = mgr.drawRayPaths(ax=ax, color="r", lw=0.3, alpha=0.5)

It seems that the problem is related to the non-zero topography (though flat) of the model. If I set the maximum elevation to zero and set to zero also the elevation of sources and geophones, rays start from the topographic surface as they should.

[halbmy]

The problem is the z column in the data file

33 # shot/geophone points
# x z [m]
66.00	7.70

If you change the second line to

#x y

all should be well. Actually, we should check this like we do in ERT.

Then all is well

I updated the notebook https://github.com/gimli-org/notebooks/blob/main/issues/642/issue642.ipynb

[anonymouscapibara]

In the end it was a trivial mistake.
Thank you very much for your help!