Note
Go to the end to download the full example code.
19. Image Domain Least-squares migration#
Seismic migration is the process by which seismic data are manipulated to create an image of the subsurface reflectivity.
In one of the previous tutorials, we have seen how the process can be formulated as an inverse problem, which requires access to a demigration-migration engine. As performing repeated migrations and demigrations can be very expensive, an alternative approach to obtain accurate and high-resolution estimate of the subsurface reflectivity has emerged under the name of image-domain least-squares migration.
In image-domain least-squares migration, we identify a direct, linear link between the migrated image \(\mathbf{m}\) and the sought after reflectivity \(\mathbf{r}\), namely:
Here \(\mathbf{H}\) is the Hessian, which can be written as:
where \(\mathbf{L}\) is the demigration operator, whilst its adjoint \(\mathbf{L}^H\) is the migration operator. In other words, we say that the migrated image can be seen as the result of a pair of demigration/migration of the reflectivity.
Whilst there exists different ways to estimate \(\mathbf{H}\), the approach that we will be using here entails applying demigration and migration to a special reflectivity model composed of regularly space scatterers. What we obtain is the spatially-varying impulse response of the migration operator, where each filter is also usually referred to as local point spread function (PSF).
Once these PSFs are computed (an operation that requires one migration and one
demigration, much cheaper than what we do in LSM), the migrated image can be deconvolved
using the pylops.signalprocessing.NonStationaryConvolve2D
operator.
import matplotlib.pyplot as plt
import numpy as np
import pylops
plt.close("all")
np.random.seed(0)
To start we create a simple model with 2 interfaces (the same we used in the LSM tutorial) and our PSF model with regularly spaced scatteres
# Velocity Model
nx, nz = 81, 60
dx, dz = 4, 4
x, z = np.arange(nx) * dx, np.arange(nz) * dz
v0 = 1000 # initial velocity
kv = 0.0 # gradient
vel = np.outer(np.ones(nx), v0 + kv * z)
# Reflectivity Model
refl = np.zeros((nx, nz))
refl[:, 30] = -1
refl[:, 50] = 0.5
# PSF Reflectivity Model
psfrefl = np.zeros((nx, nz))
psfin = (10, 15)
psfend = (-10, -5)
psfj = (30, 30)
psfx = np.arange(psfin[0], nx + psfend[0], psfj[0])
psfz = np.arange(psfin[1], nz + psfend[1], psfj[1])
Psfx, Psfz = np.meshgrid(psfx, psfz, indexing="ij")
psfrefl[psfin[0] : psfend[0] : psfj[0], psfin[1] : psfend[-1] : psfj[-1]] = 1
# Receivers
nr = 51
rx = np.linspace(10 * dx, (nx - 10) * dx, nr)
rz = 20 * np.ones(nr)
recs = np.vstack((rx, rz))
dr = recs[0, 1] - recs[0, 0]
# Sources
ns = 51
sx = np.linspace(dx * 10, (nx - 10) * dx, ns)
sz = 10 * np.ones(ns)
sources = np.vstack((sx, sz))
ds = sources[0, 1] - sources[0, 0]
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(10, 5))
axs[0].imshow(vel.T, cmap="summer", extent=(x[0], x[-1], z[-1], z[0]))
axs[0].scatter(recs[0], recs[1], marker="v", s=150, c="b", edgecolors="k")
axs[0].scatter(sources[0], sources[1], marker="*", s=150, c="r", edgecolors="k")
axs[0].axis("tight")
axs[0].set_xlabel("x [m]"), axs[0].set_ylabel("z [m]")
axs[0].set_title("Velocity")
axs[0].set_xlim(x[0], x[-1])
axs[1].imshow(refl.T, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]))
axs[1].scatter(recs[0], recs[1], marker="v", s=150, c="b", edgecolors="k")
axs[1].scatter(sources[0], sources[1], marker="*", s=150, c="r", edgecolors="k")
axs[1].axis("tight")
axs[1].set_xlabel("x [m]")
axs[1].set_title("Reflectivity")
axs[1].set_xlim(x[0], x[-1])
axs[2].imshow(psfrefl.T, cmap="gray_r", extent=(x[0], x[-1], z[-1], z[0]))
axs[2].scatter(recs[0], recs[1], marker="v", s=150, c="b", edgecolors="k")
axs[2].scatter(sources[0], sources[1], marker="*", s=150, c="r", edgecolors="k")
axs[2].axis("tight")
axs[2].set_xlabel("x [m]")
axs[2].set_title("PSF Reflectivity")
axs[2].set_xlim(x[0], x[-1])
plt.tight_layout()
We can now create our Kirchhoff modelling object which we will use to model and migrate the data, as well as to model and migrate the PSF model.
nt = 151
dt = 0.004
t = np.arange(nt) * dt
wav, wavt, wavc = pylops.utils.wavelets.ricker(t[:41], f0=20)
kop = pylops.waveeqprocessing.Kirchhoff(
z,
x,
t,
sources,
recs,
v0,
wav,
wavc,
mode="analytic",
dynamic=False,
wavfilter=True,
engine="numba",
)
kopdyn = pylops.waveeqprocessing.Kirchhoff(
z,
x,
t,
sources,
recs,
v0,
wav,
wavc,
mode="analytic",
dynamic=True,
wavfilter=True,
aperture=2,
angleaperture=50,
engine="numba",
)
d = kop @ refl
mmig = kopdyn.H @ d
dpsf = kop @ psfrefl
mmigpsf = kopdyn.H @ dpsf
fig, axs = plt.subplots(1, 2, figsize=(10, 6))
axs[0].imshow(
dpsf[ns // 2, :, :].T,
extent=(rx[0], rx[-1], t[-1], t[0]),
cmap="gray",
vmin=-1e1,
vmax=1e1,
)
axs[0].axis("tight")
axs[0].set_xlabel("x [m]"), axs[0].set_ylabel("t [m]")
axs[0].set_title(r"$d_{psf}$")
axs[1].imshow(
mmigpsf.T, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-5e0, vmax=5e0
)
axs[1].scatter(Psfx.ravel() * dx, Psfz.ravel() * dz, c="r")
axs[1].set_xlabel("x [m]"), axs[1].set_ylabel("z [m]")
axs[1].set_title(r"$m_{psf}$")
axs[1].axis("tight")
plt.tight_layout()
We can now extract the local PSFs and create the 2-dimensional non-stationary filtering operator
psfsize = (21, 21)
psfs = np.zeros((len(psfx), len(psfz), *psfsize))
for ipx, px in enumerate(psfx):
for ipz, pz in enumerate(psfz):
psfs[ipx, ipz] = mmigpsf[
int(px - psfsize[0] // 2) : int(px + psfsize[0] // 2 + 1),
int(pz - psfsize[1] // 2) : int(pz + psfsize[1] // 2 + 1),
]
fig, axs = plt.subplots(2, 1, figsize=(10, 5))
axs[0].imshow(
psfs[:, 0].reshape(len(psfx) * psfsize[0], psfsize[1]).T,
cmap="gray",
vmin=-5e0,
vmax=5e0,
)
axs[0].set_title(r"$m_{psf}$ iz=0")
axs[0].axis("tight")
axs[1].imshow(
psfs[:, 1].reshape(len(psfx) * psfsize[0], psfsize[1]).T,
cmap="gray",
vmin=-5e0,
vmax=5e0,
)
axs[1].set_title(r"$m_{psf}$ iz=1")
axs[1].axis("tight")
plt.tight_layout()
Cop = pylops.signalprocessing.NonStationaryConvolve2D(
hs=psfs, ihx=psfx, ihz=psfz, dims=(nx, nz), engine="numba"
)
mmigpsf = Cop @ refl
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].imshow(
mmig.T, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-5e1, vmax=5e1
)
axs[0].set_title(r"$m_{mig}$")
axs[0].axis("tight")
axs[1].imshow(
mmigpsf.T, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-5e1, vmax=5e1
)
axs[1].set_title(r"$m_{mig, psf}$")
axs[1].axis("tight")
plt.tight_layout()
Finally, we are ready to invert our seismic image for its corresponding
reflectivity using the pylops.optimization.sparsity.fista
solver.
minv, _, resnorm = pylops.optimization.sparsity.fista(
Cop, mmig.ravel(), eps=1e2, niter=100, eigsdict=dict(niter=5, tol=1e-2), show=True
)
minv = minv.reshape(nx, nz)
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].imshow(
mmig.T,
cmap="gray",
extent=(x[0], x[-1], z[-1], z[0]),
vmin=-5e1,
vmax=5e1,
)
axs[0].set_title(r"$m_{mig}$")
axs[0].axis("tight")
axs[1].imshow(minv.T, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-1, vmax=1)
axs[1].set_title(r"$m_{inv}$")
axs[1].axis("tight")
plt.tight_layout()
FISTA (soft thresholding)
--------------------------------------------------------------------------------
The Operator Op has 4860 rows and 4860 cols
eps = 1.000000e+02 tol = 1.000000e-10 niter = 100
alpha = 1.022879e-05 thresh = 5.114397e-04
--------------------------------------------------------------------------------
Itn x[0] r2norm r12norm xupdate
1 -0.0000e+00 2.020e+05 2.074e+05 1.754e+00
2 -0.0000e+00 1.313e+05 1.392e+05 9.219e-01
3 -0.0000e+00 8.785e+04 9.799e+04 8.151e-01
4 -0.0000e+00 6.112e+04 7.309e+04 7.085e-01
5 -0.0000e+00 4.387e+04 5.731e+04 6.224e-01
6 -0.0000e+00 3.228e+04 4.698e+04 5.523e-01
7 -0.0000e+00 2.432e+04 4.016e+04 4.923e-01
8 -0.0000e+00 1.876e+04 3.561e+04 4.391e-01
9 -0.0000e+00 1.481e+04 3.256e+04 3.923e-01
10 -0.0000e+00 1.192e+04 3.043e+04 3.523e-01
11 -0.0000e+00 9.745e+03 2.891e+04 3.186e-01
21 0.0000e+00 2.420e+03 2.408e+04 1.516e-01
31 -0.0000e+00 1.373e+03 2.203e+04 1.134e-01
41 -0.0000e+00 1.165e+03 2.029e+04 1.047e-01
51 -0.0000e+00 1.034e+03 1.899e+04 1.013e-01
61 -0.0000e+00 9.833e+02 1.784e+04 1.014e-01
71 -0.0000e+00 9.334e+02 1.702e+04 9.494e-02
81 -0.0000e+00 8.643e+02 1.650e+04 8.660e-02
91 -0.0000e+00 8.063e+02 1.603e+04 8.294e-02
92 -0.0000e+00 8.047e+02 1.599e+04 8.262e-02
93 -0.0000e+00 8.036e+02 1.595e+04 8.224e-02
94 -0.0000e+00 8.028e+02 1.590e+04 8.196e-02
95 -0.0000e+00 8.024e+02 1.586e+04 8.165e-02
96 -0.0000e+00 8.022e+02 1.582e+04 8.134e-02
97 -0.0000e+00 8.020e+02 1.577e+04 8.108e-02
98 -0.0000e+00 8.019e+02 1.573e+04 8.077e-02
99 -0.0000e+00 8.018e+02 1.569e+04 8.047e-02
100 -0.0000e+00 8.018e+02 1.565e+04 8.034e-02
Iterations = 100 Total time (s) = 2.69
--------------------------------------------------------------------------------
For a more advanced set of examples of both reflectivity and impedance image-domain LSM head over to this notebook.
Total running time of the script: (0 minutes 10.103 seconds)