How to get the Jacobian matrix for least square inversion?

[xichaoqiang]

Thanks for sharing this great package,How can I get the Jacobian matrix from the existing code for tomography inversion?

[malcolmw]

Hi, @xichaoqiang,

Can you please describe in more detail what exactly you are trying to do? You can access the gradient of the traveltime field via the EikonalSolver.traveltime.gradient attribute, but I am not sure if this is what you are looking for.

Thanks, Malcolm

[xichaoqiang]

For tomography inversion problem.How can I get the G matrix （the frechet derivative matrix/ Jacobian matrix） ?

[malcolmw]

@xichaoqiang OK, this will depend on how exactly you are formulating the inverse problem, but, assuming that your G matrix is just the sum of all the path lengths through each cell, you can obtain it as in the example below.

Try running this example, and let me know if anything is unclear.

- Malcolm

import matplotlib.pyplot as plt
import matplotlib.gridspec as gs
import numpy as np
import pykonal
import scipy.ndimage

solver = pykonal.solver.PointSourceSolver(coord_sys="cartesian")
solver.velocity.min_coords = 0, 0, 0
solver.velocity.node_intervals = 0.12524461839530332, 0.12598425196850394, 0.1
solver.velocity.npts = 512, 128, 1

for iy in range(solver.velocity.npts[1]):
    depth = solver.velocity.nodes[0, iy, 0, 1]
    solver.velocity.values[:, iy] = 3.5 + 0.1 * depth

solver.velocity.values += 20 * np.random.randn(*solver.velocity.npts)
solver.velocity.values = scipy.ndimage.gaussian_filter(solver.velocity.values, 8)

solver.velocity.values[116:204, 32:96] += 2
solver.velocity.values[308:396, 32:96] -= 2
solver.velocity.values = scipy.ndimage.gaussian_filter(solver.velocity.values, 2)

solver.velocity.values = np.clip(solver.velocity.values, 1, 8)

solver.src_loc = (32, 12, 0)

solver.solve()

n_cells_x, n_cells_y = 32, 8
cell_edges_x = np.linspace(solver.vv.min_coords[0], solver.vv.max_coords[0], n_cells_x+1)
cell_edges_y = np.linspace(solver.vv.min_coords[1], solver.vv.max_coords[1], n_cells_y+1)

gridspec = gs.GridSpec(nrows=2, ncols=1, height_ratios=[1, 0.04])

plt.close("all")
fig = plt.figure()
ax = fig.add_subplot(gridspec[0])
cax = fig.add_subplot(gridspec[1])

qmesh = ax.pcolormesh(
    solver.velocity.nodes[:, :, 0, 0],
    solver.velocity.nodes[:, :, 0, 1],
    solver.velocity.values[:, :, 0],
    cmap=plt.get_cmap("nipy_spectral_r"),
    vmin=1,
    vmax=8
)
cbar = fig.colorbar(qmesh, cax=cax, orientation="horizontal")
cbar.set_label("Velocity (km/s)")
cbar.ax.xaxis.set_minor_locator(plt.MultipleLocator(0.25))

nodes = solver.velocity.nodes

for cell_edge in cell_edges_x:
    ax.axvline(cell_edge, linewidth=0.5, color="tab:gray")

for cell_edge in cell_edges_y:
    ax.axhline(cell_edge, linewidth=0.5, color="tab:gray")

G_matrix = np.zeros((n_cells_x, n_cells_y))

for rx in range(0, 65, 4):
    ax.scatter(
        rx, 
        0,
        marker=((0, 0), (-2, 4), (2, 4), (0, 0)),
        edgecolor="k",
        facecolor="w",
        linewidth=0.5,
        s=128,
        clip_on=False,
        zorder=3
    )
    ray = solver.trace_ray(np.array([rx, 0., 0.]))
    ax.plot(ray[:, 0], ray[:, 1], color="k", linewidth=0.5)

    cell_counts, _, _ = np.histogram2d(ray[:, 0], ray[:, 1], bins=(cell_edges_x, cell_edges_y))
    G_matrix += cell_counts * solver.tt.step_size

ax.scatter(
    *solver.src_loc[:2],
    marker="*",
    edgecolor="k",
    facecolor="w",
    s=256,
    zorder=3
)

ax.invert_yaxis()
ax.set_xlabel("Horizontal offset (km)")
ax.set_ylabel("Depth (km)")
ax.set_aspect(1)
ax.set_xticks(range(0, 65, 8))
ax.set_yticks(range(0, 17, 4))

[xichaoqiang]

Thanks for your kindly rely.I will try the code.

[malcolmw]

@xichaoqiang Here is an updated example with multiple source locations

import matplotlib.pyplot as plt
import matplotlib.gridspec as gs
import numpy as np
import pykonal
import scipy.ndimage

n_srcs = 12
n_cells_x, n_cells_y = 32, 8
cell_edges_x = np.linspace(solver.vv.min_coords[0], solver.vv.max_coords[0], n_cells_x+1)
cell_edges_y = np.linspace(solver.vv.min_coords[1], solver.vv.max_coords[1], n_cells_y+1)

velocity = pykonal.fields.ScalarField3D(coord_sys="cartesian")
velocity.min_coords = 0, 0, 0
velocity.node_intervals = 0.12524461839530332, 0.12598425196850394, 0.1
velocity.npts = 512, 128, 1

for iy in range(velocity.npts[1]):
    depth = velocity.nodes[0, iy, 0, 1]
    velocity.values[:, iy] = 3.5 + 0.1 * depth

velocity.values += 20 * np.random.randn(*velocity.npts)
velocity.values = scipy.ndimage.gaussian_filter(velocity.values, 8)

velocity.values[116:204, 32:96] += 2
velocity.values[308:396, 32:96] -= 2
velocity.values = scipy.ndimage.gaussian_filter(velocity.values, 2)

velocity.values = np.clip(velocity.values, 1, 8)

plt.close("all")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

qmesh = ax.pcolormesh(
    velocity.nodes[:, :, 0, 0],
    velocity.nodes[:, :, 0, 1],
    velocity.values[:, :, 0],
    cmap=plt.get_cmap("nipy_spectral_r"),
    vmin=1,
    vmax=8
)
cbar = fig.colorbar(qmesh, ax=ax, orientation="horizontal")
cbar.set_label("Velocity (km/s)")
cbar.ax.xaxis.set_minor_locator(plt.MultipleLocator(0.25))

nodes = velocity.nodes

for cell_edge in cell_edges_x:
    ax.axvline(cell_edge, linewidth=0.5, color="tab:gray")

for cell_edge in cell_edges_y:
    ax.axhline(cell_edge, linewidth=0.5, color="tab:gray")

G_matrix = np.zeros((n_cells_x, n_cells_y))

for isrc in range(n_srcs):
    solver = pykonal.solver.PointSourceSolver(coord_sys="cartesian")
    solver.velocity.min_coords = velocity.min_coords
    solver.velocity.npts = velocity.npts
    solver.velocity.node_intervals = velocity.node_intervals
    solver.velocity.values = velocity.values
    solver.src_loc = np.random.rand(3) * [64, 16, 0]
    solver.solve()

    for rx in range(0, 65, 4):
        ax.scatter(
            rx, 
            0,
            marker=((0, 0), (-2, 4), (2, 4), (0, 0)),
            edgecolor="k",
            facecolor="w",
            linewidth=0.5,
            s=128,
            clip_on=False,
            zorder=3
        )
        ray = solver.trace_ray(np.array([rx, 0., 0.]))
        ax.plot(ray[:, 0], ray[:, 1], color="k", linewidth=0.5)

        cell_counts, _, _ = np.histogram2d(ray[:, 0], ray[:, 1], bins=(cell_edges_x, cell_edges_y))
        G_matrix += cell_counts * solver.tt.step_size

    ax.scatter(
        *solver.src_loc[:2],
        marker="*",
        edgecolor="k",
        facecolor="w",
        s=256,
        zorder=3
    )

ax.invert_yaxis()
ax.set_xlabel("Horizontal offset (km)")
ax.set_ylabel("Depth (km)")
ax.set_aspect(1)
ax.set_xticks(range(0, 65, 8))
ax.set_yticks(range(0, 17, 4))

[malcolmw]

@xichaoqiang Any luck with the code above?

[xichaoqiang]

Thanks for your kindly reply.Your code helps me alot.

[malcolmw]

Great. Closing this as solved. Let me know if you have any other questions.