2D Post-stack inversion with different forms of ADMM¶
In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
import os
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import filtfilt
from scipy.linalg import solve
from scipy.sparse.linalg import lsqr
from matplotlib.colors import LinearSegmentedColormap
from pylops.basicoperators import *
from pylops.signalprocessing import *
from pylops.utils.wavelets import *
from pylops.avo.poststack import *
from pylops.optimization.sparsity import *
from pylops.basicoperators import VStack as VStacklop
from pyproximal.proximal import *
from pyproximal import ProxOperator
from pyproximal.optimization.primal import *
In [2]:
def PSNR(x, xinv):
return 10 * np.log10(len(xinv) * np.max(xinv)**2 / np.linalg.norm(x-xinv)**2)
2D Post-stack inversion¶
We consider the following problem
$$ \mathbf{x} = arg min_\mathbf{x} f(x) + g(Dx) $$where $f=\frac{\sigma}{2} ||\mathbf{W}\mathbf{D}\mathbf{x} - \mathbf{y} ||_2^2$ and $g(Dx)=||\mathbf{D} \mathbf{x}||_1$ or $g(Dx)=||\mathbf{D} \mathbf{x}||_{2,1}$.
We can either solve this with:
- Linearized ADMM
- ADMM with the following constraint $\mathbf{D}\mathbf{x}-\mathbf{z}=0$, so $\mathbf{A}=\mathbf{D}$, $\mathbf{B}=-\mathbf{I}$, and $c=0$
Model creation¶
In [3]:
nt0, nx = 301, 251
dt0 = 0.004
t0 = np.arange(nt0)*dt0
dx = 4
x = np.arange(nx)*dx
ai = 1800 * np.ones(nt0)
ai[51:71] = 1900
ai[71:121] = 2100
ai[121:151] = 1600
ai[151:251] = 2000
ai[251:] = 2200
ai1 = np.outer(ai, np.ones(nx))
ai2d = np.roll(ai1, 32, axis=0)
ai2d[:32] = 1800
# fault
tf = 0.3
kf = 5e-4
fault = x * kf + tf
mask = np.zeros_like(ai1, dtype=int)
for ix in range(nx):
ai2d[int(fault[ix]/dt0):, ix] = ai1[int(fault[ix]/dt0):, ix]
ai2d *= 2000
m = np.log(ai2d)
In [4]:
# smooth model
nsmooth = 20
mback = filtfilt(np.ones(nsmooth)/float(nsmooth), 1, m, axis=0)
# wavelet
ntwav = 61
wav, twav, wavc = ricker(t0[:ntwav//2+1], 8)
# operator
Lop = PoststackLinearModelling(wav / 2, nt0=nt0, spatdims=nx)
# data
d = Lop * m.ravel()
d = d.reshape(nt0, nx)
# random noise
#sigman = 2e-2
#n = np.random.normal(0, sigman, d.shape)
# colored noise
sigman = 0 # 6e-2
n = filtfilt(np.ones(10)/10, 1,
filtfilt(np.ones(5)/5, 1, np.random.normal(0, sigman, (nt0, nx)).T, method='gust').T,
method='gust')
dn = d + n
fig, axs = plt.subplots(1, 3, sharey=True, gridspec_kw={'width_ratios': [3, 3, 1]}, figsize=(12, 7))
im = axs[0].imshow(np.exp(m), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[0].axvline(x[nx//2], c='k', ls='--', lw=0.5)
plt.colorbar(im, ax=axs[0])
axs[0].axis('tight')
axs[0].set_ylabel(r'$t [s]$')
axs[0].set_xlabel(r'$x [m]$')
axs[0].set_title(r'AI $[kg/(m^2s)]$', fontsize=17)
im = axs[1].imshow(dn, cmap='seismic', vmin=-0.3, vmax=0.3, extent=(x[0], x[-1], t0[-1], t0[0]))
axs[1].set_xlabel(r'$x [m]$')
axs[1].set_title('Seismic', fontsize=17)
axs[1].axis('tight')
plt.colorbar(im, ax=axs[1])
axs[2].plot(np.exp(m[:, nx//2]), t0, 'k', lw=2)
axs[2].axis('tight')
axs[2].set_title('AI trace', fontsize=17);
In [5]:
# L-ADMM (Blockiness-promoting inversion with isotropic TV)
sigma=0.01
l1 = L21(ndim=2, sigma=sigma)
l2 = L2(Op=Lop, b=dn.ravel(), x0=mback.ravel(), niter=20, warm=True)
Dop = Gradient(dims=(nt0, nx), edge=True, dtype=Lop.dtype, kind='forward')
L = 8. #np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1.
mu = 0.99 * tau / L # optimal mu<=tau/maxeig(Dop^H Dop)
mladmm = LinearizedADMM(l2, l1, Dop, tau=tau, mu=mu, x0=mback.ravel(), niter=100, show=True)[0]
dinv = Lop*mladmm
mladmm = mladmm.reshape(nt0, nx)
dinv = dinv.reshape(nt0, nx)
Linearized-ADMM
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L21.L21'>
Linear operator (A): <class 'pylops.basicoperators.vstack.VStack'>
tau = 1.000000e+00 mu = 1.237500e-01 niter = 100
Itn x[0] f g J = f + g
1 1.50964e+01 1.788e+01 1.572e+00 1.945e+01
2 1.50964e+01 1.310e+01 1.577e+00 1.467e+01
3 1.50964e+01 9.713e+00 1.616e+00 1.133e+01
4 1.50964e+01 7.329e+00 1.685e+00 9.014e+00
5 1.50964e+01 5.636e+00 1.780e+00 7.416e+00
6 1.50964e+01 4.421e+00 1.873e+00 6.294e+00
7 1.50965e+01 3.539e+00 1.953e+00 5.492e+00
8 1.50965e+01 2.889e+00 2.025e+00 4.914e+00
9 1.50965e+01 2.401e+00 2.085e+00 4.486e+00
10 1.50966e+01 2.029e+00 2.134e+00 4.162e+00
11 1.50966e+01 1.740e+00 2.173e+00 3.914e+00
21 1.50969e+01 6.439e-01 2.267e+00 2.911e+00
31 1.50968e+01 3.823e-01 2.199e+00 2.581e+00
41 1.50965e+01 2.723e-01 2.148e+00 2.421e+00
51 1.50961e+01 2.154e-01 2.099e+00 2.314e+00
61 1.50957e+01 1.800e-01 2.061e+00 2.241e+00
71 1.50954e+01 1.557e-01 2.037e+00 2.193e+00
81 1.50953e+01 1.385e-01 2.029e+00 2.167e+00
91 1.50953e+01 1.253e-01 2.011e+00 2.136e+00
92 1.50953e+01 1.242e-01 2.009e+00 2.134e+00
93 1.50953e+01 1.231e-01 2.008e+00 2.131e+00
94 1.50953e+01 1.220e-01 2.007e+00 2.129e+00
95 1.50953e+01 1.210e-01 2.006e+00 2.127e+00
96 1.50953e+01 1.200e-01 2.004e+00 2.124e+00
97 1.50954e+01 1.190e-01 2.002e+00 2.121e+00
98 1.50954e+01 1.180e-01 2.000e+00 2.118e+00
99 1.50954e+01 1.171e-01 1.999e+00 2.116e+00
100 1.50954e+01 1.162e-01 1.997e+00 2.113e+00
Total time (s) = 3.11
---------------------------------------------------------
In [6]:
fig, axs = plt.subplots(1, 3, sharey=True, gridspec_kw={'width_ratios': [3, 3, 1]}, figsize=(12, 7))
im = axs[0].imshow(np.exp(m), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[0].axvline(x[nx//2], c='k', ls='--', lw=0.5)
plt.colorbar(im, ax=axs[0])
axs[0].axis('tight')
axs[0].set_ylabel(r'$t [s]$')
axs[0].set_xlabel(r'$x [m]$')
axs[0].set_title(r'AI True', fontsize=17)
im = axs[1].imshow(np.exp(mladmm), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[1].text(0.5, 1.1, '(b)', ha='center', fontsize=16, fontweight='bold',
family='serif', transform=axs[1].transAxes)
axs[1].set_xlabel(r'$x [m]$')
axs[1].set_title('AI L-ADMM', fontsize=17)
axs[1].axis('tight')
plt.colorbar(im, ax=axs[1])
axs[2].plot(np.exp(m[:, nx//2]), t0, 'k', lw=2)
axs[2].plot(np.exp(mladmm[:, nx//2]), t0, 'r', lw=2)
axs[2].text(0.5, 1.1, '(c)', ha='center', fontsize=16, fontweight='bold',
family='serif', transform=axs[2].transAxes)
axs[2].axis('tight')
axs[2].set_title('AI trace', fontsize=17);
In [7]:
# ADMM (Blockiness-promoting inversion with isotropic TV)
sigma=0.01
l1 = L21(ndim=2, sigma=sigma)
Dop = Gradient(dims=(nt0, nx), edge=True, dtype=Lop.dtype, kind='forward')
L = 8. #np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 0.99 / L
madmm = ADMML2(l1, Lop, dn.ravel(), Dop, x0=mback.ravel(), tau=tau,
niter=100, show=True, **dict(iter_lim=20))[0]
dinv = Lop*madmm
madmm = madmm.reshape(nt0, nx)
dinv = dinv.reshape(nt0, nx)
ADMM
---------------------------------------------------------
Proximal operator (g): <class 'pyproximal.proximal.L21.L21'>
tau = 1.237500e-01 niter = 100
Itn x[0] f g J = f + g
1 1.50992e+01 2.675e+00 2.099e+00 4.774e+00
2 1.50932e+01 1.533e+00 2.226e+00 3.759e+00
3 1.50909e+01 1.036e+00 2.337e+00 3.373e+00
4 1.51019e+01 8.635e-01 2.323e+00 3.187e+00
5 1.51023e+01 7.660e-01 2.280e+00 3.046e+00
6 1.51108e+01 7.085e-01 2.254e+00 2.962e+00
7 1.51063e+01 6.662e-01 2.225e+00 2.891e+00
8 1.51099e+01 6.335e-01 2.211e+00 2.844e+00
9 1.51053e+01 6.099e-01 2.193e+00 2.803e+00
10 1.51054e+01 5.911e-01 2.171e+00 2.763e+00
11 1.51016e+01 5.760e-01 2.153e+00 2.729e+00
21 1.50973e+01 4.986e-01 2.037e+00 2.536e+00
31 1.51026e+01 4.532e-01 1.972e+00 2.425e+00
41 1.51016e+01 4.296e-01 1.936e+00 2.365e+00
51 1.51010e+01 4.135e-01 1.917e+00 2.330e+00
61 1.51003e+01 4.030e-01 1.903e+00 2.306e+00
71 1.50997e+01 3.955e-01 1.892e+00 2.288e+00
81 1.50996e+01 3.906e-01 1.882e+00 2.273e+00
91 1.51005e+01 3.846e-01 1.873e+00 2.258e+00
92 1.51005e+01 3.840e-01 1.873e+00 2.257e+00
93 1.51008e+01 3.834e-01 1.872e+00 2.255e+00
94 1.51008e+01 3.827e-01 1.871e+00 2.254e+00
95 1.51011e+01 3.821e-01 1.871e+00 2.253e+00
96 1.51011e+01 3.813e-01 1.870e+00 2.252e+00
97 1.51014e+01 3.805e-01 1.870e+00 2.250e+00
98 1.51013e+01 3.797e-01 1.869e+00 2.249e+00
99 1.51017e+01 3.790e-01 1.869e+00 2.248e+00
100 1.51016e+01 3.783e-01 1.868e+00 2.246e+00
Total time (s) = 13.02
---------------------------------------------------------
In [8]:
fig, axs = plt.subplots(1, 3, sharey=True, gridspec_kw={'width_ratios': [3, 3, 1]}, figsize=(12, 7))
im = axs[0].imshow(np.exp(m), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[0].axvline(x[nx//2], c='k', ls='--', lw=0.5)
plt.colorbar(im, ax=axs[0])
axs[0].axis('tight')
axs[0].set_ylabel(r'$t [s]$')
axs[0].set_xlabel(r'$x [m]$')
axs[0].set_title(r'AI True', fontsize=17)
im = axs[1].imshow(np.exp(madmm), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[1].text(0.5, 1.1, '(b)', ha='center', fontsize=16, fontweight='bold',
family='serif', transform=axs[1].transAxes)
axs[1].set_xlabel(r'$x [m]$')
axs[1].set_title('AI ADMM', fontsize=17)
axs[1].axis('tight')
plt.colorbar(im, ax=axs[1])
axs[2].plot(np.exp(m[:, nx//2]), t0, 'k', lw=2)
axs[2].plot(np.exp(madmm[:, nx//2]), t0, 'r', lw=2)
axs[2].text(0.5, 1.1, '(c)', ha='center', fontsize=16, fontweight='bold',
family='serif', transform=axs[2].transAxes)
axs[2].axis('tight')
axs[2].set_title('AI trace', fontsize=17);
In [9]:
# ADMM (Blockiness-promoting inversion with isotropic TV)
sigma=1e3
tv = TV(dims=(nt0, nx), sigma=sigma, rtol=1e-2)
l2 = L2(Op=Lop, b=dn.ravel(), x0=mback.ravel(), niter=20, warm=True, sigma=1e3)
L = np.real((Lop.H*Lop).eigs(neigs=1, which='LM')[0])
print(L)
tau = 0.99 / L
madmm1 = ADMM(l2, tv, x0=mback.ravel(), tau=tau, niter=100, show=True)[0]
dinv = Lop*madmm1
madmm1 = madmm1.reshape(nt0, nx)
dinv = dinv.reshape(nt0, nx)
2.0644803395254985
ADMM
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.TV.TV'>
tau = 4.795396e-01 niter = 100
Itn x[0] f g J = f + g
1 1.50998e+01 1.315e+02 3.069e+05 3.070e+05
2 1.50903e+01 5.540e+01 3.100e+05 3.100e+05
3 1.50892e+01 3.561e+01 3.062e+05 3.063e+05
4 1.50879e+01 2.732e+01 3.096e+05 3.096e+05
5 1.50881e+01 2.250e+01 3.098e+05 3.099e+05
6 1.50881e+01 1.941e+01 3.113e+05 3.113e+05
7 1.50881e+01 1.728e+01 3.112e+05 3.112e+05
8 1.50883e+01 1.564e+01 3.117e+05 3.117e+05
9 1.50883e+01 1.440e+01 3.113e+05 3.113e+05
10 1.50884e+01 1.334e+01 3.116e+05 3.116e+05
11 1.50884e+01 1.250e+01 3.112e+05 3.112e+05
21 1.50888e+01 8.114e+00 3.115e+05 3.115e+05
31 1.50892e+01 6.593e+00 3.128e+05 3.128e+05
41 1.50897e+01 5.927e+00 3.131e+05 3.131e+05
51 1.50905e+01 5.658e+00 3.130e+05 3.130e+05
61 1.50916e+01 5.742e+00 3.137e+05 3.137e+05
71 1.50930e+01 6.256e+00 3.157e+05 3.157e+05
81 1.50949e+01 7.773e+00 3.195e+05 3.195e+05
91 1.50970e+01 1.156e+01 3.259e+05 3.259e+05
92 1.50972e+01 1.216e+01 3.267e+05 3.267e+05
93 1.50975e+01 1.281e+01 3.275e+05 3.276e+05
94 1.50977e+01 1.352e+01 3.284e+05 3.284e+05
95 1.50980e+01 1.429e+01 3.294e+05 3.294e+05
96 1.50982e+01 1.513e+01 3.303e+05 3.303e+05
97 1.50985e+01 1.605e+01 3.313e+05 3.314e+05
98 1.50988e+01 1.706e+01 3.324e+05 3.324e+05
99 1.50990e+01 1.816e+01 3.335e+05 3.335e+05
100 1.50993e+01 1.936e+01 3.346e+05 3.347e+05
Total time (s) = 20.37
---------------------------------------------------------
In [10]:
fig, axs = plt.subplots(1, 3, sharey=True, gridspec_kw={'width_ratios': [3, 3, 1]}, figsize=(12, 7))
im = axs[0].imshow(np.exp(m), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[0].axvline(x[nx//2], c='k', ls='--', lw=0.5)
plt.colorbar(im, ax=axs[0])
axs[0].axis('tight')
axs[0].set_ylabel(r'$t [s]$')
axs[0].set_xlabel(r'$x [m]$')
axs[0].set_title(r'AI True', fontsize=17)
im = axs[1].imshow(np.exp(madmm1), vmin=3.2e6, vmax=4.4e6,
cmap='gist_earth_r', extent=(x[0], x[-1], t0[-1], t0[0]))
axs[1].text(0.5, 1.1, '(b)', ha='center', fontsize=16, fontweight='bold',
family='serif', transform=axs[1].transAxes)
axs[1].set_xlabel(r'$x [m]$')
axs[1].set_title('AI ADMMTV', fontsize=17)
axs[1].axis('tight')
plt.colorbar(im, ax=axs[1])
axs[2].plot(np.exp(m[:, nx//2]), t0, 'k', lw=2)
axs[2].plot(np.exp(madmm1[:, nx//2]), t0, 'r', lw=2)
axs[2].text(0.5, 1.1, '(c)', ha='center', fontsize=16, fontweight='bold',
family='serif', transform=axs[2].transAxes)
axs[2].axis('tight')
axs[2].set_title('AI trace', fontsize=17);
In [11]:
PSNR(m, mladmm), PSNR(m, madmm), PSNR(m, madmm1)
Out[11]:
(35.76777249744869, 37.23315467488398, 34.83423962919783)