Note
Go to the end to download the full example code.
Reflectivity Inversion - 3D#
This is a follow up of the Post Stack Inversion - 3D tutorial. The same seismic data will be considered, however instead of inverting it for an absolute impedance model the reflectivity is targetted.
In other words, the modelling operator becomes
\[d(t) = w(t) *r(t)\]
where \(\text{r}(t)\) is the reflectivity profile and \(w(t)\) is the time domain seismic wavelet. Similar to the other example, being this inherently a 1d operator, we can easily set up a problem where one of the dimensions (here the y-dimension) is distributed across ranks and each of them is in charge of performing modelling for a subvolume of the entire domain.
However, since a reflectivity model is sparse, a sparsity-promoting solver is used here.
We will consider two options:
- pylops_mpi.optimization.ISTA: ad-hoc distributed ISTA solver
acting like PyLops’ ISTA;
pylops_mpi.proximal.optimization.ProximalGradient: general purpose distributed Proximal Gradient solver acting like PyProximal’s ProximalGradient;
import numpy as np
from matplotlib import pyplot as plt
from mpi4py import MPI
import pyproximal
from pylops.utils.wavelets import ricker
from pylops.basicoperators import FirstDerivative
from pylops.signalprocessing import Convolve1D
import pylops_mpi
plt.close("all")
rank = MPI.COMM_WORLD.Get_rank()
size = MPI.COMM_WORLD.Get_size()
Let’s start by defining all the parameters required by the
pylops.avo.poststack.PoststackLinearModelling operator.
# Model
model = np.load("../testdata/avo/poststack_model.npz")
x, z, m = model['x'][::3], model['z'], np.log(model['model'])[:, ::3]
# Making m a 3D model
ny_i = 20 # size of model in y direction for rank i
y = np.arange(ny_i)
m3d_i = np.tile(m[:, :, np.newaxis], (1, 1, ny_i)).transpose((2, 1, 0))
ny_i, nx, nz = m3d_i.shape
# Size of y at all ranks
ny = MPI.COMM_WORLD.allreduce(ny_i)
# Wavelet
dt = 0.004
t0 = np.arange(nz) * dt
ntwav = 41
wav, _, wavc = ricker(t0[:ntwav // 2 + 1], 15)
# Collecting m3d at all ranks
m3d = np.concatenate(MPI.COMM_WORLD.allgather(m3d_i))
We now create the linear operators to model the data (including a time derivative as in the post-stack tutorial) as well as that to invert the data for the underlying reflectivity model.
# Create flattened model data
m3d_dist = pylops_mpi.DistributedArray(global_shape=ny * nx * nz)
m3d_dist[:] = m3d_i.flatten()
# LinearOperator Derivative + Convolve
Dop = FirstDerivative((ny_i, nx, nz), axis=-1)
Cop = Convolve1D((ny_i, nx, nz), wav, offset=wavc, axis=-1)
DDiag = pylops_mpi.basicoperators.MPIBlockDiag(ops=[Dop, ])
CDiag = pylops_mpi.basicoperators.MPIBlockDiag(ops=[Cop, ])
# Reflectivity
r_dist = DDiag @ m3d_dist
r_local = r_dist.local_array.reshape((ny_i, nx, nz))
r = r_dist.asarray().reshape((ny, nx, nz))
# Data
d_dist = CDiag @ r_dist
d_local = d_dist.local_array.reshape((ny_i, nx, nz))
d = d_dist.asarray().reshape((ny, nx, nz))
We now perform sparsity-promotion inversion with the ISTA solver
FISTA (soft thresholding)
--------------------------------------------------------------------------------
The Operator Op has 2937000 rows and 2937000 cols
eps = 2.000000e-02 tol = 1.000000e-08 niter = 200
alpha = 2.150970e-02 thresh = 2.150970e-04
--------------------------------------------------------------------------------
Itn x[0] r2norm r12norm xupdate
1 -4.5123e-04 2.863e+03 3.378e+03 2.448e+01
2 -6.4525e-04 1.471e+03 2.052e+03 6.164e+00
3 -4.0424e-04 8.674e+02 1.496e+03 4.476e+00
4 -0.0000e+00 5.565e+02 1.215e+03 3.542e+00
5 7.4901e-05 3.782e+02 1.056e+03 2.938e+00
6 0.0000e+00 2.690e+02 9.595e+02 2.501e+00
7 -0.0000e+00 1.982e+02 8.964e+02 2.172e+00
8 0.0000e+00 1.503e+02 8.516e+02 1.921e+00
9 0.0000e+00 1.166e+02 8.182e+02 1.725e+00
10 0.0000e+00 9.216e+01 7.919e+02 1.570e+00
11 -0.0000e+00 7.412e+01 7.704e+02 1.443e+00
21 -0.0000e+00 1.886e+01 6.515e+02 8.884e-01
31 -0.0000e+00 1.213e+01 6.005e+02 6.544e-01
41 -0.0000e+00 1.016e+01 5.771e+02 5.203e-01
51 -0.0000e+00 8.928e+00 5.624e+02 4.726e-01
61 -0.0000e+00 8.137e+00 5.505e+02 4.646e-01
71 -0.0000e+00 7.296e+00 5.428e+02 4.203e-01
81 -0.0000e+00 6.926e+00 5.371e+02 3.805e-01
91 -0.0000e+00 6.644e+00 5.339e+02 3.171e-01
101 -0.0000e+00 6.626e+00 5.320e+02 2.625e-01
111 -0.0000e+00 6.528e+00 5.300e+02 2.430e-01
121 -0.0000e+00 6.487e+00 5.287e+02 2.263e-01
131 -0.0000e+00 6.415e+00 5.277e+02 2.083e-01
141 -0.0000e+00 6.440e+00 5.266e+02 2.036e-01
151 -0.0000e+00 6.468e+00 5.256e+02 2.024e-01
161 -0.0000e+00 6.475e+00 5.248e+02 1.985e-01
171 -0.0000e+00 6.505e+00 5.241e+02 1.967e-01
181 -0.0000e+00 6.518e+00 5.234e+02 1.955e-01
191 -0.0000e+00 6.529e+00 5.228e+02 1.943e-01
192 -0.0000e+00 6.531e+00 5.227e+02 1.941e-01
193 -0.0000e+00 6.532e+00 5.227e+02 1.939e-01
194 -0.0000e+00 6.535e+00 5.226e+02 1.937e-01
195 -0.0000e+00 6.538e+00 5.225e+02 1.937e-01
196 -0.0000e+00 6.542e+00 5.225e+02 1.937e-01
197 -0.0000e+00 6.546e+00 5.224e+02 1.937e-01
198 -0.0000e+00 6.550e+00 5.224e+02 1.937e-01
199 -0.0000e+00 6.554e+00 5.223e+02 1.938e-01
200 -0.0000e+00 6.559e+00 5.222e+02 1.939e-01
Iterations = 200 Total time (s) = 62.39
--------------------------------------------------------------------------------
And now with the Proximal Gradient solver
l2d = pylops_mpi.proximal.MPIL2(Op=CDiag, b=d_dist, x0=r0_dist)
l1 = pyproximal.L1(sigma=1e-2)
l1d = pylops_mpi.proximal.MPIProxOperator(l1)
CDiag1 = CDiag.H @ CDiag
maxeig = np.abs(
pylops_mpi.optimization.eigs.power_iteration(
CDiag1,
b_k=r0_dist.empty_like(),
dtype=CDiag1.dtype
)[0]
)
rpg3d_dist = pylops_mpi.proximal.optimization.primal.ProximalGradient(
l2d, l1d, x0=r0_dist, tau=.99 / maxeig, niter=200, acceleration="fista",
show=True
)
rpg3d = rpg3d_dist.asarray().reshape((ny, nx, nz))
Accelerated Proximal Gradient
---------------------------------------------------------
Proximal operator (f): <MPIL2>
Proximal operator (g): <MPIProxOperator (L1)>
tau = 0.021296768599487922 epsg = 1.0
niter = 200 tol = None
niterback = 100 acceleration = fista
Itn x[0] f g J=f+eps*g tau
1 -4.46760e-04 2.909e+03 2.551e+02 3.164e+03 2.130e-02
2 -6.41387e-04 1.491e+03 2.896e+02 1.780e+03 2.130e-02
3 -4.12640e-04 8.781e+02 3.139e+02 1.192e+03 2.130e-02
4 -0.00000e+00 5.633e+02 3.289e+02 8.923e+02 2.130e-02
5 8.40024e-05 3.829e+02 3.385e+02 7.215e+02 2.130e-02
6 0.00000e+00 2.724e+02 3.451e+02 6.175e+02 2.130e-02
7 -0.00000e+00 2.008e+02 3.490e+02 5.498e+02 2.130e-02
8 0.00000e+00 1.522e+02 3.507e+02 5.029e+02 2.130e-02
9 0.00000e+00 1.181e+02 3.509e+02 4.690e+02 2.130e-02
10 0.00000e+00 9.339e+01 3.500e+02 4.434e+02 2.130e-02
21 -0.00000e+00 1.901e+01 3.167e+02 3.357e+02 2.130e-02
41 -0.00000e+00 1.019e+01 2.837e+02 2.939e+02 2.130e-02
61 -0.00000e+00 8.156e+00 2.713e+02 2.795e+02 2.130e-02
81 -0.00000e+00 6.940e+00 2.652e+02 2.721e+02 2.130e-02
101 -0.00000e+00 6.629e+00 2.627e+02 2.694e+02 2.130e-02
121 -0.00000e+00 6.493e+00 2.611e+02 2.676e+02 2.130e-02
141 -0.00000e+00 6.436e+00 2.601e+02 2.666e+02 2.130e-02
161 -0.00000e+00 6.468e+00 2.592e+02 2.657e+02 2.130e-02
181 -0.00000e+00 6.518e+00 2.585e+02 2.650e+02 2.130e-02
192 -0.00000e+00 6.531e+00 2.581e+02 2.647e+02 2.130e-02
193 -0.00000e+00 6.533e+00 2.581e+02 2.646e+02 2.130e-02
194 -0.00000e+00 6.534e+00 2.581e+02 2.646e+02 2.130e-02
195 -0.00000e+00 6.536e+00 2.580e+02 2.646e+02 2.130e-02
196 -0.00000e+00 6.540e+00 2.580e+02 2.645e+02 2.130e-02
197 -0.00000e+00 6.544e+00 2.580e+02 2.645e+02 2.130e-02
198 -0.00000e+00 6.547e+00 2.579e+02 2.645e+02 2.130e-02
199 -0.00000e+00 6.551e+00 2.579e+02 2.645e+02 2.130e-02
200 -0.00000e+00 6.555e+00 2.579e+02 2.644e+02 2.130e-02
Total time (s) = 28.07
---------------------------------------------------------
Finally, we display the modeling and inversion results
if rank == 0:
# Check the distributed implementation gives the same result
# as the one running only on rank0
Dop0 = FirstDerivative((ny, nx, nz), axis=-1)
Cop0 = Convolve1D((ny, nx, nz), wav, offset=wavc, axis=-1)
r0 = Dop0 @ m3d
d0 = Cop0 @ r0
# Check the two distributed implementations give the same modelling results
print('Reflectivity Distr == Local', np.allclose(r, r0))
print('Data Distr == Local', np.allclose(d, d0))
# Visualize
fig, axs = plt.subplots(nrows=5, ncols=3, figsize=(9, 14), constrained_layout=True)
axs[0][0].imshow(m3d[5, :, :].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
axs[0][0].set_title("Model x-z")
axs[0][0].axis("tight")
axs[0][1].imshow(m3d[:, 200, :].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
axs[0][1].set_title("Model y-z")
axs[0][1].axis("tight")
axs[0][2].imshow(m3d[:, :, 220].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
axs[0][2].set_title("Model y-z")
axs[0][2].axis("tight")
axs[1][0].imshow(r[5, :, :].T, cmap="gray", vmin=-.1, vmax=.1)
axs[1][0].set_title("Reflectivity Model x-z")
axs[1][0].axis("tight")
axs[1][1].imshow(r[:, 200, :].T, cmap="gray", vmin=-.1, vmax=.1)
axs[1][1].set_title("Reflectivity Model y-z")
axs[1][1].axis("tight")
axs[1][2].imshow(r[:, :, 220].T, cmap="gray", vmin=-.1, vmax=.1)
axs[1][2].set_title("Reflectivity Model y-z")
axs[1][2].axis("tight")
axs[2][0].imshow(d[5, :, :].T, cmap="gray", vmin=-1, vmax=1)
axs[2][0].set_title("Data x-z")
axs[2][0].axis("tight")
axs[2][1].imshow(d[:, 200, :].T, cmap='gray', vmin=-1, vmax=1)
axs[2][1].set_title('Data y-z')
axs[2][1].axis('tight')
axs[2][2].imshow(d[:, :, 220].T, cmap='gray', vmin=-1, vmax=1)
axs[2][2].set_title('Data x-y')
axs[2][2].axis('tight')
axs[3][0].imshow(rfista3d[5, :, :].T, cmap='gray', vmin=-.1, vmax=.1)
axs[3][0].set_title("FISTA Reflectivity iter x-z")
axs[3][0].axis("tight")
axs[3][1].imshow(rfista3d[:, 200, :].T, cmap='gray', vmin=-.1, vmax=.1)
axs[3][1].set_title('FISTA Reflectivity iter y-z')
axs[3][1].axis('tight')
axs[3][2].imshow(rfista3d[:, :, 220].T, cmap='gray', vmin=-.1, vmax=.1)
axs[3][2].set_title('FISTA Reflectivity iter x-y')
axs[3][2].axis('tight')
axs[4][0].imshow(rpg3d[5, :, :].T, cmap='gray', vmin=-.1, vmax=.1)
axs[4][0].set_title("PG Reflectivity iter x-z")
axs[4][0].axis("tight")
axs[4][1].imshow(rpg3d[:, 200, :].T, cmap='gray', vmin=-.1, vmax=.1)
axs[4][1].set_title('PG Reflectivity iter y-z')
axs[4][1].axis('tight')
axs[4][2].imshow(rpg3d[:, :, 220].T, cmap='gray', vmin=-.1, vmax=.1)
axs[4][2].set_title('PG Reflectivity iter x-y')
axs[4][2].axis('tight')
plt.savefig('Reflectivity')
Reflectivity Distr == Local True
Data Distr == Local True
To run this tutorial with our NCCL backend, refer to Reflectivity Inversion with NCCL tutorial in the repository.
Total running time of the script: (1 minutes 35.396 seconds)