Fat ray tomography

[makeabhishek]

Hi @halbmy,
Thanks for providing the details about fat ray tomography implementation.

I was trying fat ray tomography and following #563 #564 and the references mentioned. I'm a little confuse how the weights are assigned to traveltime inversion. I understood that the weights are calculated to obtained the Fresnel volume, but it is unclear how they are used in slowness estimation. Are there detailed equations/a description of how the slowness field s is calculated using w?

In traveltime tomography the traveltime is calculated using

$$t_q = \int_{L_q} \frac{1}{v(l)} dl \approx {\sum_i} {s_i} {l_i}$$

Where, $v$ is the velocity and $L_q$ the total travel path of this ray, $l_i$ are the individual ray paths in each cell, and $s_i$ are the discrete and constant slowness values. To do a forward computation, the traveltimes through each discrete cell can be determined geometrically and formulated in matrix notation as:

$$ t = G s$$

where, $t$ (i.e., first arrival traveltimes), are related to a set of model parameters $s$ (i.e., slowness),

Here, the cumulated traveltimes of all rays form the vector $t$, the slowness values of all cells are located in s and the ray paths of every ray form the matrix $G$. The matrix $G$ contains one row per ray and one column per model cell. Accordingly, measured traveltimes $t$ are associated with data vector $d$, and slowness values s are associated with the model parameters m. In summary, in the ray-approximation, the Jacobian $(G)$ consists of the ray segments per model cell. Multiplying this with the slowness vector $(s)$ produces the total traveltime of the ray along each segment. And the ray segment in each model cell actually reflects the "sensitivity", as longer ray segments have a larger impact on the traveltime.

As given in Jordi et al [2016], "The main difference between the classical ray-based and fat ray tomography approaches is the calculation of the sensitivities. For ray tomography the sensitivities effectively are the lengths of the ray segments within the respective model cells. In the fat ray tomography scheme the sensitivities are calculated using fat rays without resorting to ray tracing." For a 2D fat ray in a discretized medium the sensitivity of the traveltime $t_i$ with respect to a change in slowness $s_k$ is given by:

$$\frac{\partial t_i}{\partial s_k} = \frac{a_{k}^{i}}{A_{i}} \cdot \frac{t_i}{s_k}$$

Is this sensitivity refer to matrix ${l_i}$

Furthermore, in Pygimli, fatray is implemented by weighting the Jacobian or senitivity, as,

pg.show(mgr.paraDomain, mgr.fop.J[7])
pg.show(mgr.paraDomain, mgr.fop.FresnelWeight[7])

The Jacobian (first command) contains the ray length in meters and the weight (FresnelWeight) is a dimensionless function summing to 1 that is being used to weight the jacobian.

How did the "sensitivity", or "Jacobian" or "Frechet derivatives" incorporated in traveltime equation to obtain the slowness value? Could you refer me to the main source code?

Ref:

1. Jordi, Claudio, Cedric Schmelzbach, and Stewart Greenhalgh. "Frequency-dependent traveltime tomography using fat rays: application to near-surface seismic imaging." Journal of Applied Geophysics 131 (2016): 202-213.

[halbmy]

When using the "fat ray" option, the forward calculation is NOT affected, but is still the same fastest ray approximation. Only the Jacobian matrix is changed in a way that a larger sensitive volume (the fresnel volume) is affected and should stabilize the inversion process but probably not lead to much different velocity models. The whole methodology is explained by Jordi et al. (2016) and exactly implemented in the source code under

gimli/pygimli/physics/traveltime/modelling.py

Line 154 in c8de564

def createJacobian(self, slowness):

[makeabhishek]

Thanks for sharing the details @halbmy. As mentioned in the Pygimli for the Jacobian matrix calculation of fat ray tomography,

self.J[i] = wa * tsr / slowness
self.FresnelWeight[i] = wa

I'm also adding the reference for the notebook to calculate the Jacobian . A good notebook to understand the calculation of Jacobian. Here I do not get, what is $\partial p_j$ in $S_{ij}=\frac{\partial r_i}{\partial p_j}$

1. I did not understand why do we have to multiple wa with tsr / slowness?

2. I tried to find the Jacobian callculation for thin-ray tomography (using Dijkstra Modelling), but unable to find the exact path. May I also know the source code to define Jacobian for thin ray tomography?

[halbmy]

First 2. the Dijkstra algorithm in implemented in the C++ core so you can't see it. But it is just pure geometry, calculating the length of the ray through every cell.

As to 1., it is exactly equation (4) of Jordi et al. (2016): the first term is the (dimensionless) weight and the second ratio is the traveltime between source and receiver divided by the slowness. It is nicely described and easily implemented.

[makeabhishek]

oh yeah, Thanks for pointing this out. I didn't recognize the terms $\frac{a_{k}^{i}}{A_i}$ from Eq. (4) in Jordi et al. (2016), which is given by

$$\frac{\partial t_i}{\partial s_k}= \frac{a_{k}^{i}}{A_i} \cdot \frac{t_i}{s_k}$$

where $a_{k}^{i}$ is the area of the i-th fat ray width on the k-th model cell and $A_i$ is the total area of the i-th fat ray.
So, basically $\frac{a_{k}^{i}}{A_i}$ is the weighting term $w$, which is defined by Eq. (5) in Jordi et al. (2016).

In the ray-approximation, the Jacobian consists of the ray segments per model cell. Multiplying this with the slowness vector produces the total traveltime of the ray along each segment.

1. Are C++ codes not open sourced?
2. As, asked in previous thread, in the notebook to calculate the Jacobian , what is $\partial p_j$ in $S_{ij}=\frac{\partial r_i}{\partial p_j}$. Is there any symbol mismatch in the notebook?

[halbmy]

1. Of course the C++ code is open source, just look into https://github.com/gimli-org/gimli/blob/master/core/src/ttdijkstramodelling.cpp
2. The parameter p is the intrinsic parameter being used, for traveltime this is the slowness, the model m is typically its logarithm or log-based bound function.