How to get the Jacobian matrix with the honor grid? about fteikpy HOT 6 CLOSED

Hi @yhu5945,

Let's consider the ray traced below (blue solid line). The blue dots are the outputs of the ray tracing using the option honor_grid=True.

1. For each segment along the ray, calculate the midpoint coordinates (red dots),
2. Use the midpoints coordinates to find the indices of the cell transected by the segments (e.g., using numpy.searchsorted),
3. Knowing the indices of the cells, set the values (of the Jacobian matrix) of the cells transected by this raypath to the length of the segments.

This is roughly how I would approach this problem. Obviously you should repeat these steps for each pair of source/receiver (i.e. row of your Jacobian matrix).

from fteikpy.

Hi @keurfonluu,

Thank you so much for your reply. I will work on this and see how to implement it efficiently. :)

from fteikpy.

Hi @yhu5945,

I had a bit of spare time this weekend to test the approach. Please find below a simple code that shows how you can calculate the Jacobian matrix using fteikpy (I only demonstrate the approach for a single pair of source/receiver, you need to repeat that process for each pair you have in your dataset).

import numpy

from fteikpy import Eikonal2D


# Velocity model
vel = numpy.ones((100, 100))
dz, dx = 0.1, 0.1

# Solve Eikonal at source
eik = Eikonal2D(vel, gridsize=(dz, dx))
src = numpy.random.rand(2)
tt = eik.solve(src, return_gradient=True)

# Raytracing
rcv = numpy.random.rand(2) * 10.0
ray = tt.raytrace(rcv, honor_grid=True)

# Find cell indices along raypath
midpoints = 0.5 * (ray[1:] + ray[:-1])
iz = numpy.searchsorted(tt.zaxis, midpoints[:, 0], side="right") - 1
ix = numpy.searchsorted(tt.xaxis, midpoints[:, 1], side="right") - 1

# Integrate along raypath
M = numpy.zeros_like(vel)
M[iz, ix] = numpy.linalg.norm(ray[1:] - ray[:-1], axis=1)

# Check results
print(f"Analytical: {numpy.linalg.norm(src - rcv)}")
print(f"Calculated (interpolation): {tt(rcv)}")
print(f"Calculated (raytracing): {M.ravel() @ (1.0 / vel.ravel())}")

In eik.solve, I enabled the option return_gradient=True as it provides a more accurate traveltime gradient grid (so more accurate ray tracing).

Comparing the results with the analytical solution, you can see that both traveltime values calculated using either the interpolation or raytracing are fairly accurate (and almost exact if source and receiver are placed on a grid point).
M.ravel() is one row of your Jacobian matrix that corresponds to this pair of source/receiver in the linear system Ls=t (where L is the Jacobian matrix, s the slowness vector and t the traveltime vector).

Note that I have fixed a couple of issues I encountered when implemented this code, so please update to fteikpy==2.2.0 before trying this code.

from fteikpy.

Hi @keurfonluu,

I really appreciate your help and detailed explanation! Thanks!

from fteikpy.

Hi @keurfonluu,
When I run the code below, a RuntimeError always appears (see below). It seems not working for any receivers. Would you please help me?

import numpy as np
from fteikpy import Eikonal3D
import matplotlib.pyplot as plt

# The size of Rock
l_z,l_x,l_y = 125,50,50
dims=(l_z,l_x,l_y)

# Velocity model
vel = np.ones((125, 50,50))*2.5
vel[50:100,:,:]=3.0;
vel[100:-1,:,:]=3.5;
dz,dx, dy = 1, 1, 1
z=np.arange(0, 125.+dz,dz)-62.5
x=np.arange(0, 50.+dx,dx)-25
y=np.arange(0, 50.+dy,dy)-25

# Probe position
Nrz,Nrc = 9,10
radius_mm=25
theta=np.linspace(0,2*np.pi,Nrc+1)
theta=np.delete(theta,Nrc)
zr=np.linspace(6.5,118.5,Nrz)-l_z/2
xr=radius_mm*np.cos(theta)
yr=radius_mm*np.sin(theta)
probe_posxy=np.tile(np.vstack((xr,yr)).T,(Nrz,1))
probe_posz=np.tile(zr.T,(Nrc,1)).T.reshape((Nrz*Nrc,1))
probe_pos=np.hstack((probe_posz,probe_posxy))    

# Source coordinates
src=probe_pos
Ns=len(src)

# Receiver coodinates
rcv = list(range(0,Ns))
for i in np.arange(0,Nrc):
    for j in np.arange(0,Nrz):
        temp = np.delete(probe_pos,np.arange(0,Ns,Nrc)+i,axis=0) 
        rcv[Nrc*j+i] = temp

# Solve Eikonal at source
eik = Eikonal3D(vel, gridsize=(dz, dx,dy),origin=(-62.5,-25,-25))
ns=0
tt = eik.solve(src[ns], return_gradient=True)

# Raytracing
ray = tt.raytrace(rcv[ns], honor_grid=True)

from fteikpy.

Hi @David-zw,

Sorry for the delay. The fix for your issue is now available in version 2.4.0. Please update to this version:

pip install fteikpy -U

from fteikpy.