Attenuation tomography

[sminize]

Hello everyone,
is there a way to do attenuation tomography in pygimli? I want to use the calculated ray paths from travel time tomography (see example below) and invert the attenuation with these ray paths. Is this somehow possible?

https://www.pygimli.org/_examples_auto/seismics/plot_3_crosshole_tomography.html

Best regards

Simon

[halbmy]

Of course, attenuation tomography should be possible. If you take the logarithmized amplitudes, the forward problem is a linear one. So we just have to extract and set up a linear modelling operator by fAtt=pg.core.LinearModelling(fTT.jacobian()). Would be great to test this if you have an example data set. However, often the source amplitude is not well known and would have to be added to inversion.

[sminize]

Hi Mr. Günther,
Thank you for your reply. The travel time tomography problem can be reconstructed with the following code (it is not the latest pygimili version, sorry for that) :

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics.traveltime import [TravelTimeManager](url)

bh_spacing = 500
bh_length = 500
sensor_spacing = 10
world = mt.createRectangle(start=[0, (bh_length + 0.01)], end=[bh_spacing, 0.0],
                           marker=0)
scheme.load("travel_time.dat")

data=scheme
tt = TravelTimeManager()
tt.applyData(data)
refinement = 1
x = np.arange(0, bh_spacing + refinement, sensor_spacing * refinement)
y = np.arange(0.0, bh_length + 3, sensor_spacing * refinement)
mesh = pg.meshtools.createGrid(x, y)

ax, _ = pg.show(mesh, hold=True)
ax.plot(A[:, 0], A[:, 1], "ro")

tt.inv.inv.setRecalcJacobian(False) #With RAY BENDING  

invmodel = tt.invert(data, mesh=mesh, secNodes=3, lam=800, zWeight=1,
                     useGradient=False, verbose=False)

fig, (ax2) = plt.subplots(1, 1, figsize=(8, 7), sharex=True, sharey=True)

tt.showResult(ax=ax2, logScale=False, nLevs=3)
tt.drawRayPaths(ax=ax2, color="0.8", alpha=0.3)
fig.tight_layout()

The travel time data, with source, receiver positions and the error is loaded with scheme.load("travel_time.dat"). The data is given in the appendix.

Now i can set up a linear model operator using:

fAtt=pg.core.LinearModelling(tt.fop.jacobian())
scheme.load("atten.dat")
data=scheme

Furthermore, i can load the attenuation data together with the receiver, source location and error with scheme.load("atten.dat").
(the data format is the same as before, only the travel time data is replaced by the attenuation values).

I am not sure how to proceed now. I want to use the same Jacobian matrix (and ray paths) and reconstruct the attenuation on the same mesh.
Can you help me with this? Thank you in advance.

Best regards

Simon

[sminize]

The data can be found in these dropbox links. The data comes from a simulation

https://www.dropbox.com/s/oap44vlfi33pwpn/travel_time.dat?dl=0

https://www.dropbox.com/s/yxzzvvzgu1pfos5/atten.dat?dl=0

[halbmy]

I don't understand why you need two files for it. Can't you just add the amplitude data as a data column to the travel times? The shot/geophone combinations need to be the same anyway. In the two file I found:

(traveltime)
# g s err t valid 
11	1	1.06916526174729e-02	7.50443089958069e-02	1
12	1	1.03458263087364e-02	7.43526563783340e-02	1
...

and

(attenuation)
# g s err t valid 
11	1	7.01652617472872e-04	7.50443089958069e-02	1
12	1	3.55826308736440e-04	7.43526563783340e-02	1
...

The traveltime is identical but what does the error mean?

[sminize]

Thank you very much for your reply. I changed the data file in such a way that now the attenuation values (indicated by 'a') are also present in the datacontainer. I included before only the travel time because my procedure is based on the crosshole tomography example (https://www.pygimli.org/_examples_auto/seismics/plot_3_crosshole_tomography.html) and my datacontainer had the same structure.

I recognized that the attenuation values before were not correct and i changed them in the new file. Since i calculated them analytically there is no noise in the data of the attenuation values.
The error in the data refers to the travel times. I picked the p wave on the measured x,y data and calculated the standard deviation. I don't know if this is a valid way to estimate the error. :)

https://www.dropbox.com/s/joewciqirqggmdf/data.dat?dl=0

Best regards and thank you for your help

Simon Schmid

[sminize]

[sminize]

Here ist what the model and the inversion of the travel times looks like

[sminize]

The inversion is very prone to how the error is set.

[halbmy]

I am a bit confused. In the file there are amplitude data:

...
# g s a err t valid 
11	1	5.00000000000000e+01	8.74565771841089e-04	7.46120261098863e-02	1
12	1	5.00000000000000e+01	6.15196040288768e-04	7.43526563783340e-02	1
13	1	5.00000000000000e+01	5.71967751696715e-04	7.53901353045433e-02	1
14	1	5.00000000000000e+01	5.28739463104663e-04	7.73786365797778e-02	1
15	1	5.00000000000000e+01	3.12598020144380e-04	8.04046167812216e-02	1
16	1	5.00000000000000e+01	9.64565771841047e-05	8.36035101370337e-02	1
17	1	5.00000000000000e+01	5.32282885920524e-05	8.74075995331345e-02	1

but they seem not to be depending on the travel path length. How did you generate them? If assuming a damping coefficient in dB/m (i.e. a relative amplitude decrease per meter ray path) one can use the Jacobian matrix (J=fop.jacobian()) of the traveltime simulation (or inversion) directly to compute a relative amplitude a=np.exp(J*np.alpha). Or am I wrong?

[sminize]

I only took into account the length of the path (with straight rays) intersecting the anomaly. I corrected this in the new data set. :). I also increased the error a little bit (the inversion result was better).
https://www.dropbox.com/s/pk4fwtlmybq3q1d/data_new.dat?dl=0

The attenuation values were calculated as follows:
d=legth of the ray
f=fraction of the ray intersection the anomaly
a1=attenuation coefficient of material 1
a2=attenuation coefficient of material 2
att=d_i*(f_i*a1+(1-f_i)*a2)

I hope that clarifies it. This is just a "toy example", which should illustrate the general workflow of attenuation tomography based on the rays determined with travel time tomography.

Alternatively of course the jacobian matrix can be used. Which would also take into account bent rays.
Of course, it does not make sense to generate data with straight rays and invert it with bent rays. However, i just want to understand the workflow for doing the inversion.

Thank you very much for your help.

[halbmy]

Just to get things right, the attenuation coefficient a defines the relative loss of amplitude A per ray path length

a = - 1/A dA/dl

So the amplitude after passing a cell is A=A0 exp(- a l) (A0 being the initial amplitude).

The logarithmized relative amplitude is

log(A/A0) = - integral a dl

and can be expressed as matrix-vector product of the attenuation coefficient vector by the Jacobian (way path) matrix -W a.

But how do you know the initial amplitude A0? And do you have to take into account spherical attenuation (i.e. energy distributes on spheres of increasing size)?

[sminize]

Thank you for clarifying. I think I understood it now. I calculated the relative amplitude using the jacobian. As an initial amplitude, i assumed 1. That makes it easier I think. In the data, the logarithmized relative amplitude is given. I am not sure if i have to take into account spherical attenuation through geometrical spreading. Is that not automatically included in the attenuation coefficient, since it is one of the attenuation mechanisms.

https://www.dropbox.com/s/mgr4z0tnto3x5ah/data_with_Jacobian.dat?dl=0

However, I have still problems doing the inversion with this data using the jacobian from the travel time tomography.
Can you help me with this?

Best regards

Simon Schmid

[halbmy]

The point about the initial amplitude A0 is purely practical, because you usually do cannot measure it at the transmitter. So it would be an additional unknown for the inversion but it could be well constrained in case of many receivers per shot.

The geometric spreading (energy balance) is of course not included in the normal attenuation coefficient, because even in case of no attenuation you have a geometric spreading.

So for the synthetics it all does not matter and you can make your modellings, but it will not be applicable to real data.

[sminize]

Thank you for that remark. When i want to use it on real data i have to put more thoughts into it. :)

[halbmy]

Can you attach what you have right now and I'm gonna have a look?

[sminize]

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics.traveltime import TravelTimeManager
import numpy as np
import matplotlib.pyplot as plt
import pickle


bh_spacing = 500
bh_length = 500
sensor_spacing = 10
world = mt.createRectangle(start=[0, (bh_length + 0.01)], end=[bh_spacing, 0.0],
                           marker=0)

file_name = "datalist.pkl"

open_file = open(file_name, "rb")
loaded_list = pickle.load(open_file)
open_file.close()

scheme = pg.DataContainer()

A=np.concatenate((loaded_list[6],loaded_list[7]), axis=0)

for sen in A:
    scheme.createSensor(sen)

scheme.resize(len(loaded_list[0]))
scheme["s"] = loaded_list[0]
scheme["g"] = loaded_list[1]
scheme["valid"] = loaded_list[2]
scheme.registerSensorIndex("s")
scheme.registerSensorIndex("g")
scheme['t']=loaded_list[3]
scheme['a']=loaded_list[4]
scheme['err']=loaded_list[5]



data=scheme
tt = TravelTimeManager()
tt.applyData(data)
refinement = 1
x = np.arange(0, bh_spacing + refinement, sensor_spacing * refinement)
y = np.arange(0.0, bh_length + 3, sensor_spacing * refinement)
mesh = pg.meshtools.createGrid(x, y)

ax, _ = pg.show(mesh, hold=True)
ax.plot(A[:, 0], A[:, 1], "ro")

tt.inv.inv.setRecalcJacobian(False) #With RAY BENDING  

invmodel = tt.invert(data, mesh=mesh, secNodes=3, lam=800, zWeight=1,
                     useGradient=False, verbose=False)


fig, (ax2) = plt.subplots(1, 1, figsize=(8, 7), sharex=True, sharey=True)

tt.showResult(ax=ax2, logScale=False, nLevs=3)
tt.drawRayPaths(ax=ax2, color="0.8", alpha=0.3)
fig.tight_layout()

#the travel time tomography could probably be improved a bit adapting the inversion parameters
#Here the attenuation tomography should be done. The data
fAtt=pg.core.LinearModelling(tt.fop.jacobian())  #Get jacobian from Travel Time

[sminize]

I uploaded the data in a pickle file. It should be possible to load it with the given code.

https://www.dropbox.com/s/9eqdegqs8wsao50/datalist.pkl?dl=0

[halbmy]

I am still not getting the point in your amplitude data. They are between 1.4 and 18, but they should be below 1. Looking at the data using ax, cb = pg.viewer.mpl.showDataContainerAsMatrix(data, "g", "s", "a") yields

A forward calculation would look like:

J = mgr.fop.jacobian()
att = np.ones(mesh.cellCount()) * 0.005
amp = np.exp(-J.mult(att))
ax, cb = pg.viewer.mpl.showDataContainerAsMatrix(data, "g", "s",
                                                  amp, label="amplitude") 

yielding

i.e. the amplitude dropping to values between 4% and 70% of the original amplitude, which is looking similar to your values but the other way round.

[sminize]

I think i found the mistake. I forgot the - in the formula. Thank you for clarifying that. Here is the corrected data:

https://www.dropbox.com/s/lot1nbkiipjuwoo/datalist.pkl?dl=0

Best regards

Simon

[halbmy]

For sake of open science and user friendliness, I would prefer to generate an Ascii-readable data container and attach this file here instead of generating some binary data and upload it to some proprietary service where it is probably not available anymore soon.

[sminize]

That makes sense. It will be easier to reproduce for others in that way. Do you know how to export the data container as an Ascii-readable data container?

[halbmy]

data.save("filename.dat")

[sminize]

data.txt

[sminize]

I uploaded the data in a .txt file. :)

[sminize]

data.zip

[sminize]

The data is also uploaded as a zip file because the .dat seems not to work with github.

In order to unzip following can be used:

import zipfile 
with zipfile.ZipFile("data.zip","r") as zip_ref: 
                    zip_ref.extractall("./")

However, i am not sure how to load it. The following does not work correctly:

data=pg.load("./data.dat")

[halbmy]

OK, here's the code for the attenuation tomography

fAtt=pg.core.LinearModelling(mgr.fop.jacobian(), -1)  # take traveltime jacobian multiplied by -1 (attenuation)
fAtt.setMesh(mesh)
att = pg.Vector(mesh.cellCount(), 0.01)
tLog = pg.trans.TransLog()
fAtt.modelTrans = tLog
inv = pg.Inversion(fop=fAtt)
inv.transModel = tLog
dataVec = - np.log(data["a"])  # a linear problem if log-scaled amplitudes are used
error = np.ones_like(dataVec) * 0.01
model = inv.run(dataVec, error, startModel=att)
pg.show(mesh, model)

yielding

Here's the complete code: pg311.zip

[sminize]

Thank you very much for your help. ;) This was very useful for me.