Implement tomography operator

[pavane]

I would like to implement a tomography operator like as follows
tt=F*V (forward Operator) (tt travel time table, v velocity model, F eikonal operator)
V=F^{-1}*tt (inverse operator)
I implemented the forward operator as a 3D eikonal operator (matvec)
I would like help in implementing the adjoint operator ( rmatvec)

[cako]

Hi @pavane, it would be excellent to add more operators to PyLops. In order to implement adjoint operators, I would recommend reading some of the code the library has. In general it is not entirely trivial to implement optimized adjoint operators, but there are definitely some guidelines on it. For example, if you have a for loop like this within your code:

for i in range(n):
    y[i] += A * x[i]

then the adjoint would look something like this:

for i in range(n)[::-1]:
    x[i] += A.T * y[i]

Note that the order of the loop has been interchanged to preserve the order of the summation. So one technique to calculate the adjoint is to go through the code from bottom to top, following rules such as the above. Peng Shen's thesis (PDF) contains this and many other rules in its fourth chapter; Jon Claerbout's Basic Earth Imaging (PDF) and Earth Soundings Analysis contain several examples.

With that said, before you begin, you must ensure that your algorithm is actually linear. In your case I believe you are trying to do linearized traveltime tomography, in which case you must ensure that your algorithm is compatible with that. If you need references I am happy to provide some!

[pavane]

Thank you so much for the references. I am familiar with the theory of adjoint operators, I am struggling with how to implement the operator in pylops. I will try to understand the rules for writing adjoint operators from the references. Yes, I am trying to implement linearized travel time tomography. I would be grateful if you could pass me the references.

[cako]

Great, if you have any specific piece of code you would like us to look at, feel free to share and I can take a look at it.

My understanding of linearized traveltime tomography is the following. Start with a nonlinear relationship between slowness s and traveltime t:

where L(s) is the raypath and ℓ is the arclength. This equation is nonlinear because the integration region (raypath) depends on the slowness itself. The idea is then to linearize this relation near a certain slowness s_0 yielding

or equivalently

Knowing the raypath for s_0, it becomes easy to write the forward operator:

for x0 in raypath_positions_0:
    delta_t +=delta_s[x0]

So it's just summing along certain indices. In this case, the adjoint is quite simple:

for x0 in raypath_positions_0:
    delta_s[x0] = delta_t

So you are just "spreading" the traveltime error along the raypath. For many traveltime measurements, just add another for loop and substitute the equality in the adjoint code for +=

I think the standard tome for this is Gus Nolet's 1987 Seismic Tomography. Another great reference is Rawlinson and Sambridge's Seismic traveltime tomography of the crust and lithosphere which is freely available here (PDF).

However, when you are doing eikonal solving you don't know the raypaths themselves. In this case, I admit I'm not too knowledgeable about how it's done, but maybe the following references will help as they all contain linearization of the eikonal equation:

Fomel, S., 2000, Traveltime computation with the linearized eikonal equation, SEP Report

Alkhalifah, T., 2002, Traveltime computation with the linearized eikonal equation for anisotropic media, Geophysical Prospecting

Aldridge, D., 1994, Linearization of the eikonal equation, Geophysics

[pavane]

Thank you for your detailed response. I will try to come up with a code example based on Fomel 2000 for your review. Thanks again.

…

On Sun, Jun 6, 2021 at 3:25 PM Carlos da Costa ***@***.***> wrote: Great, if you have any specific piece of code you would like us to look at, feel free to share and I can take a look at it. My understanding of linearized traveltime tomography is the following. Start with a nonlinear relationship between slowness s and traveltime t: <https://www.codecogs.com/eqnedit.php?latex=t&space;=&space;%5Cint_%7BL(s)%7D&space;s(x)d%5Cell> where L(s) is the raypath and ℓ is the arclength. This equation is nonlinear because the integration region (raypath) depends on the slowness itself. The idea is then to linearize this relation near a certain slowness s_0 yielding <https://www.codecogs.com/eqnedit.php?latex=t&space;-&space;t_0&space;=&space;%5Cint_%7BL(s_0)%7D(s(x)-s_0(x))%5C,%5Cmathrm%7Bd%7D%5Cell> or equivalently <https://www.codecogs.com/eqnedit.php?latex=%5Cdelta&space;t&space;=&space;%5Cint_%7BL(s_0)%7D%5Cdelta&space;s(x)%5C,%5Cmathrm%7Bd%7D%5Cell> Knowing the raypath for s_0, it becomes easy to write the forward operator: for x0 in raypath_positions_0: delta_t +=delta_s[x0] So it's just summing along certain indices. In this case, the adjoint is quite simple: for x0 in raypath_positions_0: delta_s[x0] = delta_t So you are just "spreading" the traveltime error along the raypath. For many traveltime measurements, just add another for loop and substitute the equality in the adjoint code for += I think the standard tome for this is Gus Nolet's 1987 *Seismic Tomography*. Another great reference is Rawlinson and Sambridge's *Seismic traveltime tomography of the crust and lithosphere* which is freely available here (PDF) <http://web.gps.caltech.edu/classes/ge193.old/papers_books/rawlinson_sambridge2003.pdf> . However, when you are doing eikonal solving you don't know the raypaths themselves. In this case, I admit I'm not too knowledgeable about how it's done, but maybe the following references will help as they all contain linearization of the eikonal equation: Fomel, S., 2000, *Traveltime computation with the linearized eikonal equation* <http://sepwww.stanford.edu/oldsep/sergey/sepsergey/lineiko.pdf>, SEP Report Alkhalifah, T., 2002, *Traveltime computation with the linearized eikonal equation for anisotropic media* <https://onlinelibrary.wiley.com/doi/full/10.1046/j.1365-2478.2002.00322.x>, Geophysical Prospecting Aldridge, D., 1994, *Linearization of the eikonal equation* <https://library.seg.org/doi/10.1190/1.1443552>, Geophysics — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#236 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAYL3HYZHCMJF243HYOYFATTRPRV7ANCNFSM457RRXCQ> .

[mrava87]

Great @pavane, this will make it a great addition. Carlos references are very good ones to understand how to implement the adjoint of the eikonal solver. We are here at any time when you have a first version to review and help :)

[mrava87]

Hi again,
I was actually reading again your initial message and the comment of Carlos and I think there is something important to keep in mind here.

If you are after a traveltime tomography with rays the approach that Carlos describes is perfect and can be easily implemented within the context/framework of PyLops. I have a very simple (non optimized) example in my Seismology teaching material https://github.com/DIG-Kaust/Seismology/blob/main/SeismicTomography/SeismicTomography.ipynb

However, if you are after eikonal-based traveltime tomography I would say that this falls more within the nonlinear adjoint-state type of inversion which doesn't really belong to PyLops (or at least not yet). There has been some discussion in the past to extend the PyLops framework to nonlinear inverse problems where instead of forward and adjoint you would need forward and gradient but this has not yet been done and needs some proper thinking and prototyping of the overall extension before deciding whether it is worth or not having it.

The references Carlos provided discuss the so-called linearized eikonal equation which to my understanding is just a way to make the eikonal solution cheaper by taking advantage of an initial guess. But as you can see the forward and adjoint are linked to the traveltime, not the slowness you wish to solve for. A good reference which I guess goes more in the direction of what you have in mind is https://academic.oup.com/gji/article/209/3/1629/3072667 and allows you to solve the traveltime=Eikonal(slowness) problem via adjoint-state formulation. But similar to Carlos I also don't have first hand experience on implementing this so I may also be missing something here :)

[cako]

Revisiting this to say that the PyKonal project also allows for a posteriori ray-tracing, meaning, you can use the eikonal solver to get the traveltime map, and then use that map to solve the ray equations. If you use those rays for traveltime tomography, effectively you are doing linearized, ray-based, shortest-time, traveltime tomography.

[mrava87]

Nice :) I think the fact you can get rays from the eikonal map is very good and basically the last step would be to grid them into the tomographic grid.

I am not in favour of using PyKonal directly in PyLops as it doesn’t seem to be on PyPi and this would add overhead on our installation pipeline. Said that, if someone is interested to create an example and add it to pylops_notebooks, that would be nice 😊

[cako]

Oh, for sure. Just something I remembered while looking at this issue :)

[mrava87]

If there is interest I think we could revisit this as I see PyKonal is now on PyPi. I haven’t check if the installation works but definitely worth a look as it would allow us to have rays for tt tomo that scikit-fmm cannot provide 😀

[fpicetti]

Not sure if applicable here, but Ettore Biondi has a nice implementation of TT tomography in OccamyPy. It's not based on pykonal, as far as I rembember, but @pavane maybe it's worth having a look at it!

[mrava87]

Thanks @fpicetti! I remember seeing it, indeed this may be a good thing to look at @pavane .if you are still interested :)