Attenuation seismic tomography

[kerim371]

Hi,

I'm trying to do some experiments with amplitude tomography.
I've seen topics #311 and #736 and it seems that some expertise is required.

Forward problem:
Lets suppose that we know exactly velocity and attenuation (say Q-factor, bigger Q -> less attenuation).
For such velocity model we calculate raypaths.
For these raypaths we can calculate amplitude of waves for each Shot-Receiver pair using Q-model.

Inverse problem for Q-factor only:
We need to recover Q-model using true raypaths and amplitudes calculated from forward problem.

The difference between this description and #311 is that in the latter there is no raypaths involved in calculation of amplitudes but rather shot-receiver distances.

I'm trying to implement this in a single notebook.
The problems:

1. how to get true raypaths (it should be done using Dijkstra algorithm but I can't find it)
2. how to calculate amplitudes using raypaths. Each raypaths consists of a number of XY points. It should be done similar to calculating traveltimes but the problem is that for each segment of a ray we need to know Q-value from Q-model. Don't know how to get it.

For now I follow crosswhole tutorial and I get:

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

# 1. Generate synthetic model
# Acquisition parameters
bh_spacing = 20.0
bh_length = 25.0
sensor_spacing = 2.5

world = mt.createRectangle(start=[0, -(bh_length + 3)], end=[bh_spacing, 0.0],
                           marker=0)

depth = -np.arange(sensor_spacing, bh_length + sensor_spacing, sensor_spacing)

sensors = np.zeros((len(depth) * 2, 2))  # two boreholes
sensors[len(depth):, 0] = bh_spacing  # x
sensors[:, 1] = np.hstack([depth] * 2)  # y

# Create forward model and mesh
c0 = mt.createCircle(pos=(7.0, -10.0), radius=3, nSegments=25, marker=1)
c1 = mt.createCircle(pos=(12.0, -18.0), radius=4, nSegments=25, marker=2)
geom = world + c0 + c1
for sen in sensors:
    geom.createNode(sen)

mesh_fwd = mt.createMesh(geom, quality=34, area=0.25)
vp = np.array([2000., 2300, 1700])[mesh_fwd.cellMarkers()]      # Velocity
q = np.array([50., 100, 30])[mesh_fwd.cellMarkers()]            # Q-factor
# ax, cb = pg.show(mesh_fwd, vp, logScale=False,
#                  label=pg.unit('vel'), cMap=pg.cmap('vel'), nLevs=3)

# Model visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
pg.show(mesh_fwd, vp, ax=ax1, label="Vp [m/s]", cMap="viridis")
pg.show(mesh_fwd, q, ax=ax2, label="Q factor", cMap="jet_r", logScale=True)
plt.suptitle("Synthetic models Vp and Q")
plt.show()

scheme = tt.createCrossholeData(sensors)

# ====================================
# next I need to calculate raypaths only without inversion
# ====================================

[halbmy]

Just some thoughts to get started (I never did amplitude tomography myself):

1. You need the ray paths for doing amplitude tomography. Therefore, you first have to invert the first-arrivals, thus obtaining a velocity model and an accompagnying way (ray path) matrix J = mgr.fop.jacobian().
2. If you logarithmize the amplitudes, then the forward problem for attenuation is a linear one, i.e. d=J*m where m is the relative energy loss, i.e. $2\pi/Q$. One has to setup a forward operator using pg.core.LinearModelling(mesh, J) and then use an inversion.
3. Formally writing down the important equations here will help discussion.

[halbmy]

Ah yes, I am sorry. I introduced this class in the C++ core but those core classes will not work with the modern inversion classes in pyGIMLi >=1.2. We have to use pg.frameworks.LinearModelling instead, implemented in https://github.com/gimli-org/gimli/blob/6c1c689f5ba473ba2ed0cb464dd5b126c54e0689/pygimli/frameworks/modelling.py#L445

[kerim371]

@halbmy hi, thank you.

I will just leave the formulas while I'm struggling with this.

The travel times are computed according to this formula

$$ T(x) = \sum_{i=1}^{N} \frac{\Delta s_i}{v_i} $$

Where:

• $T(x)$ - travel time at the distance $x$
• $\Delta s_i$ - length of the segment in i-th element
• $v_i$ - velocity in i-th element
• $N$ - the number of elements that the ray travels through

The attenuation through the Q-factor is computed in the following way

It is worth noting that attenuation is a frequency dependent. So we need to know with wich frequency we work. In practice that can be some central frequency of a signal.

$$ A(x) = A_0 \cdot \exp \left( -\pi f \sum_{i=1}^{N} \frac{\Delta s_i}{v_i Q_i} \right) = A_0 \cdot \exp \left( -\pi f \sum_{i=1}^{N} \frac{t_i}{Q_i} \right) $$

Where:

• $A(x)$ - amplitude at the distance $x$
• $A_0$ - source amplitude (at distance 0)
• $f$ - frequency (Hz)
• $\Delta s_i$ - length of the segment in i-th element
• $v_i$ - velocity in i-th element
• $Q_i$ - quality factor in i-th element
• $N$ - the number of elements

We can linearize it by taking the logarithms:

$$ \ln\left(\frac{A(x)}{A_0}\right) = -\pi f \sum_{i=1}^{N} \frac{\Delta s_i}{v_i Q_i} = -\pi f \sum_{i=1}^{N} \frac{t_i}{Q_i} $$

finally we can try make the attenuation formula to very similar to travel times formula:

$$ \frac{\ln\left(\frac{A(x)}{A_0}\right)}{-\pi f } = \sum_{i=1}^{N} \frac{\Delta s_i}{v_i Q_i} = \sum_{i=1}^{N} \frac{t_i}{Q_i} $$

[halbmy]

Very well written. So I see three possibilities:

1. Use pg.frameworks.LinearModelling(mesh, fop.jacobian()) and invert for the parameter $1/vQ$ and convert the result into Q.
2. Generate a matrix A=pg.matrix.MultRightMatrix(fop.jacobian(), slowness) and use this with LinearModelling(mesh, A), inverting for 1/Q.
3. Use an inverse transformation and combine it with LinearModelling through PetroModelling, inverting for Q.

I guess option 2 could be best as Q (and thus 1/Q) can be best constrained.

[kerim371]

Thank you for hints.
The only note is that pg.frameworks.LinearModelling(mesh, fop.jacobian()) takes only one positional argument:

Init signature: pg.frameworks.LinearModelling(A)
Docstring:      Modelling class for linearized problems with a given matrix.
Init docstring: Initialize by storing the (reference to the) matrix.
File:           [c:\users\user\miniconda3\envs\pg\lib\site-packages\pygimli\frameworks\modelling.py](file:///C:/users/user/miniconda3/envs/pg/lib/site-packages/pygimli/frameworks/modelling.py)
Type:           class
Subclasses:

[kerim371]

@halbmy it seems I was able to solve forward problem. But can't understand how to solve inverse problem when we know raypaths and we only need to invert attenuation.

def calculate_attenuation(ray_nodes, mesh, vp_model, q_model, freq=50):
    """
    Calculate attenuation along the ray path considering mesh cells
    :param ray_nodes: ray node coordinates (Nx2 array)
    :param mesh: computational mesh
    :param vp_model: velocity model
    :param q_model: Q-factor model
    :param freq: signal frequency (Hz)
    :return: attenuation coefficient, total path length
    """
    total_attenuation = 1.0
    total_length = 0.0
    
    # Iterate through all ray segments
    for i in range(len(ray_nodes)-1):
        start_point = ray_nodes[i]
        end_point = ray_nodes[i+1]
        segment_length = np.linalg.norm(end_point - start_point)
        
        # Find midpoint of the segment
        mid_point = (start_point + end_point) / 2
        
        # Find cell containing the midpoint
        cell = mesh.findCell(mid_point)
        if cell:
            v = vp_model[cell.id()]
            q = q_model[cell.id()]
            alpha = np.pi * freq / (q * v)
            total_attenuation *= np.exp(-alpha * segment_length)
            total_length += segment_length
    
    return total_attenuation, total_length

# Example usage for all rays
freq = 50  # Signal frequency (Hz)
attenuations = []
lengths = []

# Get ray paths
nodes = mgr.fop._core.mesh().positions(withSecNodes=True)
shots = mgr.fop.data.id("s")  # Shot indices
recei = mgr.fop.data.id("g")  # Receiver indices

for i in range(mgr.fop.data.size()):
    # Get ray nodes
    ray_nodes = nodes[mgr.fop.way(shots[i], recei[i])]
    
    # Calculate attenuation
    att, length = calculate_attenuation(ray_nodes, mesh_fwd, vp_model, q_model, freq)
    attenuations.append(att)
    lengths.append(length)

[kerim371]

Don't know why I cant get satisfiable result:

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

############### FORWARD MODELING ###################

# 1. Generate synthetic model
# Acquisition parameters
bh_spacing = 20.0
bh_length = 25.0
sensor_spacing = 2.5

world = mt.createRectangle(start=[0, -(bh_length + 3)], end=[bh_spacing, 0.0],
                           marker=0)

depth = -np.arange(sensor_spacing, bh_length + sensor_spacing, sensor_spacing)

sensors = np.zeros((len(depth) * 2, 2))  # two boreholes
sensors[len(depth):, 0] = bh_spacing  # x
sensors[:, 1] = np.hstack([depth] * 2)  # y

# Create forward model and mesh
c0 = mt.createCircle(pos=(7.0, -10.0), radius=3, nSegments=25, marker=1)
c1 = mt.createCircle(pos=(12.0, -18.0), radius=4, nSegments=25, marker=2)
geom = world + c0 + c1
for sen in sensors:
    geom.createNode(sen)

mesh_fwd = mt.createMesh(geom, quality=34, area=0.25)
vp_model = np.array([2000., 2300, 2500])[mesh_fwd.cellMarkers()]      # Velocity
q_model = np.array([100, 10, 1])[mesh_fwd.cellMarkers()]            # Q-factor
# ax, cb = pg.show(mesh_fwd, vp_model, logScale=False,
#                  label=pg.unit('vel'), cMap=pg.cmap('vel'), nLevs=3)

# Model visualization
fig = plt.figure(figsize=(10, 8))
gs = fig.add_gridspec(1, 2, width_ratios=[1, 1])
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1])

# Visualize models
pg.show(mesh_fwd, vp_model, ax=ax1, label="Vp [m/s]", cMap="viridis")
pg.show(mesh_fwd, q_model, ax=ax2, label="Q factor", cMap="jet_r", logScale=True)
plt.show()

# set mesh for inverse proble same as forward problem mesh
mesh = mesh_fwd

scheme = tt.createCrossholeData(sensors)
mgr = tt.TravelTimeManager()
mgr.applyData(scheme)
# mgr.setMesh(mesh, secNodes=4)
mgr.applyMesh(mesh, secNodes=4, ignoreRegionManager=True)
mgr.fop.createJacobian(1/vp_model)
t_calc = mgr.fop.response(1.0 / vp_model)

def calculate_attenuation(ray_nodes, mesh, vp_model, q_model, freq=50):
    """
    Calculate attenuation along the ray path considering mesh cells
    :param ray_nodes: ray node coordinates (Nx2 array)
    :param mesh: computational mesh
    :param vp_model: velocity model
    :param q_model: Q-factor model
    :param freq: signal frequency (Hz)
    :return: attenuation coefficient, total path length
    """
    total_attenuation = 1.0
    total_length = 0.0
    
    # Iterate through all ray segments
    for i in range(len(ray_nodes)-1):
        start_point = ray_nodes[i]
        end_point = ray_nodes[i+1]
        segment_length = np.linalg.norm(end_point - start_point)
        
        # Find midpoint of the segment
        mid_point = (start_point + end_point) / 2
        
        # Find cell containing the midpoint
        cell = mesh.findCell(mid_point)
        if cell:
            v = vp_model[cell.id()]
            q = q_model[cell.id()]
            alpha = np.pi * freq / (q * v)
            total_attenuation *= np.exp(-alpha * segment_length)
            total_length += segment_length
    
    return total_attenuation, total_length

# Example usage for all rays
freq = 50  # Signal frequency (Hz)
attenuations = []
lengths = []

# Get ray paths
nodes = mgr.fop._core.mesh().positions(withSecNodes=True)
shots = mgr.fop.data.id("s")  # Shot indices
recei = mgr.fop.data.id("g")  # Receiver indices

for i in range(mgr.fop.data.size()):
    # Get ray nodes
    ray_nodes = nodes[mgr.fop.way(shots[i], recei[i])]
    
    # Calculate attenuation
    att, length = calculate_attenuation(ray_nodes, mesh_fwd, vp_model, q_model, freq)
    attenuations.append(att)
    lengths.append(length)

# Visualization of results
plt.figure(figsize=(10, 5))
plt.scatter(lengths, attenuations, c=t_calc, cmap='viridis')
plt.colorbar(label='Travel time (s)')
plt.xlabel('Ray path length (m)')
plt.ylabel('Attenuation coefficient')
plt.title('Attenuation vs ray path length')
plt.grid(True)
plt.show()

from matplotlib import cm
import matplotlib.colors as colors
from matplotlib.colorbar import ColorbarBase

# Normalize attenuation values for color scale
attenuations_normalized = (np.array(attenuations)-np.min(attenuations))/(np.max(attenuations)-np.min(attenuations))

# Create color map
cmap = cm.plasma_r
norm = colors.Normalize(vmin=np.min(attenuations), vmax=np.max(attenuations))

# Create figure with space for colorbar
fig = plt.figure(figsize=(10, 8))
gs = fig.add_gridspec(1, 3, width_ratios=[1, 1, 0.05])
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1])
cax = fig.add_subplot(gs[2])

# Visualize models
pg.show(mesh_fwd, vp_model, ax=ax1, label="Vp [m/s]", cMap="viridis")
pg.show(mesh_fwd, q_model, ax=ax2, label="Q factor", cMap="jet_r", logScale=True)

# Plot rays with attenuation-based coloring
for ray_idx in range(shots.size()):
    ray_nodes = nodes[mgr.fop.way(shots[ray_idx], recei[ray_idx])]
    ax1.plot(ray_nodes[:,0], ray_nodes[:,1], 
             linewidth=1.5, 
             c=cmap(attenuations_normalized[ray_idx]))
    ax2.plot(ray_nodes[:,0], ray_nodes[:,1], 
             linewidth=1.5, 
             c=cmap(attenuations_normalized[ray_idx]))

# Add colorbar
cb = ColorbarBase(cax, cmap=cmap, norm=norm, orientation='vertical')
cb.set_label('Attenuation Coefficient', rotation=270, labelpad=15)

plt.suptitle("Vp and Q Models with Ray Paths (color = attenuation)")
plt.tight_layout()
plt.show()

############### INVERSE PROBLEM ###################

A = pg.matrix.MultRightMatrix(mgr.fop.jacobian(), 1/vp_model)
fop = pg.frameworks.LinearModelling(A)

dataVec = - np.log(attenuations) / (np.pi * freq)
errorVec = np.ones_like(dataVec) * 0.001
inv = pg.Inversion(fop, maxIter=5, verbose=True)
inv.errorVals = errorVec

inv_model = inv.run(dataVals=dataVec, startModel = 1/vp_model, showProgress=True, verbose=True, lam=20, maxIter=20, mesh=mesh)

# 7. Data fit analysis
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(dataVec, dataVec, 'bo', label='Measured')
ax.plot(dataVec, inv.response, 'rx', label='Predicted')
# ax.plot([min(d), max(d)], [min(d), max(d)], 'k--')
ax.set_xlabel("Measured (-ln(attenuation))")
ax.set_ylabel("Predicted (-ln(attenuation))")
ax.set_title("Data Fit")
ax.legend()
ax.grid(True)
plt.show()

ax,cbar = pg.show(mesh, 1/inv_model)

At the end amplitudes coincide pretty good but the image is bad.

[halbmy]

Thank you for the code. It is just the regularisation. You need to provide the mesh to the forward operator. Then it should work.

[halbmy]

I would try the following

...
velocity = mgr.invert(...)
A = pg.matrix.MultRightMatrix(fop.jacobian(), slowness)
f = pg.frameworks.LinearModelling(A)
f.setMesh(mesh)
inv = pg.Inversion(fop)
...

[kerim371]

Thank you very much!

It seems to be working

[halbmy]

I tried to write it down in my own words: the traveltime for one ray is

$$ t = \sum_{i=1}^{M} \frac{l_i}{v_i}=\sum_{i=1}^{M} l_i s_i $$

expressed for all rays by matrix vector product

$$ {\bf t} = {\bf W} \cdot {\bf s} $$

amplitude for one ray expressed by Q (energy loss in one cycle) in a way element

$$ A = A_0 \cdot \exp \frac{-\pi f t}{Q} = A_0 \cdot \exp \frac{-\pi f l s}{Q} $$

and for a sum

$$ A = A_0 \prod_{i=1}^{N} \exp\left(\frac{-\pi f l_i s_i}{Q_i}\right) = A_0 \cdot \exp \left( \sum_{i=1}^{N} \frac{-\pi f l_i s_i}{Q_i} \right) $$

Normalizing with $A_0$ and taking the logarithm we obtain

$$ \ln \frac{A}{A_0} = \sum_{i=1}^{N} \frac{-\pi f l_i s_i}{Q_i} $$

So it can be expressed by a matrix-vector product using the model parameters $m_i=1/Q_i$:

$$ {\bf f(m)} = {\bf A}\cdot{\bf m} = {\bf W}\cdot\text{diag}(-\pi f s_i)\cdot{\bf m} $$

In a pure forward modelling code this reads:

mm = -np.pi*freq/vp_model/q_model
att = mgr.fop.jacobian().dot(mm)

but the result is

i.e. clearly depending on the traveltime and almost not on the Q factor.

The full forward operator needed for inversion would be:

G = pg.matrix.MultRightMatrix(mgr.fop.jacobian(), -np.pi*freq/vp_model)
f = pg.frameworks.LinearModelling(G)

and one can obtain the same solution by f(1/q_model) and further use it in an inversion

f.setMesh(mesh)
inv = pg.Inversion(fop=f)
q = 1 / inv.run(np.log(A), startModel=0.001, verbose=True)

[halbmy]

Your forward response looks, on the contrary, like this

which makes really sense.

[halbmy]

One more comment: a frequency of 50Hz does not make much sense for this crosshole geometry. A velocity of 2000m/s would imply a wavelength of 40m, much bigger than the size of the model. You would need frequencies in the range of kHz.

[kerim371]

@halbmy thank you for explanation.

There is some things I can't understand.

I have some messed up in my mind, I will try to ask some basics.

The problem we are working on:

1. solve forward problem (compite amplitude of attenuated signal)
2. from the given velocity model and attenuated signals restore original Q-factor model

In my code I compute forward problem by hand using calculate_attenuation function wich is not perfect (it is in python and it iterates through the rays). As I guess we can solve this using forward operators from pygimly. But mm = -np.pi*freq/vp_model/q_model; att = mgr.fop.jacobian().dot(mm) gives the "diagonal" picture you have shown and it is far from more accurate one provided by you later.

Is that ok? Because from my understanding if we apply forward operator on our data we should get the correct attenuation and "diagonal" image doesn't look like correct attenuation.

Basically is there a way to run forward operator to get correct image for attenuated signals llike on this image?