Poststack

Post-stack inversion of Volve data

Author: M. Ravasi, KAUST

Welcome to the Matrix-free inverse problems with PyLops tutorial!

The aim of this tutorial is to:

  • Utilize PyLops to perform seismic post-stack inversion of the Volve dataset;
  • Familiarize with novel optimization techniques that can enhance the quality of your inversion using L1 regularization;
  • Get introduced to PyProximal, PyLops little sister library to perform even more advanced optimization by means of so-called proximal algorithms.

Post-stack seismic modelling is the process of constructing seismic post-stack data from an acoustic impedance profile as function of time (or depth). From a mathematical point of view, this can be easily achieved using the following forward model:

d(t,θ)=12w(t)∗dln(m(t))dtd(t, \theta) = \frac{1}{2} w(t) * \frac{dln(m(t))}{dt}
(1)

where m(t)m(t) is the acoustic impedance profile and w(t)w(t) is the time domain seismic wavelet.

In compact form:

d=WDm\mathbf{d}= \mathbf{W} \mathbf{D} \mathbf{m}
(2)

The associated inverse problems we wish to solve can written as:

  • Least-squares

    J=argminm∣∣d−WDm∣∣2+ϵ1∣∣∇x,ym∣∣2 J = arg min_{\mathbf{m}} ||\mathbf{d} - \mathbf{W} \mathbf{D} \mathbf{m}||_2 + \epsilon_1||\nabla_{x,y}\mathbf{m}||_2
    (3)

  • L1-regularized

    J=argminm∣∣d−WDm∣∣2+ϵ1∣∣∇x,ym∣∣2+ϵ2∣∣Dm∣∣1 J = arg min_{\mathbf{m}} ||\mathbf{d} - \mathbf{W} \mathbf{D} \mathbf{m}||_2 + \epsilon_1||\nabla_{x,y}\mathbf{m}||_2 + \epsilon_2||\mathbf{D}\mathbf{m}||_1
    (4)

  • TV-regularized

    J=argminm∣∣d−WDm∣∣2+ϵTV(m) J = arg min_{\mathbf{m}} ||\mathbf{d} - \mathbf{W} \mathbf{D} \mathbf{m}||_2 + \epsilon TV(\mathbf{m})
    (5)

where ∇x,y\nabla_{x,y} is the Laplacian along the x,y plane.

Data retrieval

In order to be able to run this notebook we need to first make sure we have access to the Volve dataset.

We cannot provide the data directly given its size, but the Volve data is hosted on a Azure Blob storage so it is very easy to download it using the Azure CLI (https://github.com/Azure/azure-cli or pip install azure-cli would do).

First of all let's investigate what is present in the Seismic directory

#!az storage blob list --account-name datavillagesa --container-name volve --prefix Seismic/ --sas-token "$YOURTOKEN" > ../data/seismicinversion/list_seismic.txt
!head -n 60 data/list_seismic.txt
[
  {
    "container": "volve",
    "content": "",
    "deleted": null,
    "encryptedMetadata": null,
    "encryptionKeySha256": null,
    "encryptionScope": null,
    "isAppendBlobSealed": null,
    "isCurrentVersion": null,
    "lastAccessedOn": null,
    "metadata": {},
    "name": "Seismic/2D/ST0299/Navigation/ST0299-CMP.p190",
    "objectReplicationDestinationPolicy": null,
    "objectReplicationSourceProperties": [],
    "properties": {
      "appendBlobCommittedBlockCount": null,
      "blobTier": "Hot",
      "blobTierChangeTime": null,
      "blobTierInferred": true,
      "blobType": "BlockBlob",
      "contentLength": 237136,
      "contentRange": null,
      "contentSettings": {
        "cacheControl": null,
        "contentDisposition": null,
        "contentEncoding": null,
        "contentLanguage": null,
        "contentMd5": "YycxlvGGZqp4s7nW0TH3vQ==",
        "contentType": "text/plain; charset=utf-8"
      },
      "copy": {
        "completionTime": null,
        "destinationSnapshot": null,
        "id": null,
        "incrementalCopy": null,
        "progress": null,
        "source": null,
        "status": null,
        "statusDescription": null
      },
      "creationTime": "2020-10-09T11:29:25+00:00",
      "deletedTime": null,
      "etag": "0x8D86C4693A8213F",
      "lastModified": "2020-10-09T11:29:25+00:00",
      "lease": {
        "duration": null,
        "state": "available",
        "status": "unlocked"
      },
      "pageBlobSequenceNumber": null,
      "pageRanges": null,
      "rehydrationStatus": null,
      "remainingRetentionDays": null,
      "serverEncrypted": true
    },
    "rehydratePriority": null,
    "requestServerEncrypted": null,
    "snapshot": null,
    "tagCount": null,

where you will need to substitute $YOURTOKEN with your personal token. To get a token, simply register at https://data.equinor.com/authenticate and find the token in the red text in the Usage section. Ensure to copy everything from ?sv= to =rl in place of $YOURTOKEN.

We can now download the file of interest

# TO BE RUN ONLY ONCE TO RETRIEVE THE DATA
#!az storage blob download --account-name datavillagesa --container-name volve --name Seismic/ST10010/Stacks/ST10010ZC11_PZ_PSDM_KIRCH_FULL_D.MIG_FIN.POST_STACK.3D.JS-017536.segy --file ST10010ZC11_PZ_PSDM_KIRCH_FULL_D.MIG_FIN.POST_STACK.3D.JS-017536.segy --sas-token "$YOURTOKEN"
#!mv ST10010ZC11_PZ_PSDM_KIRCH_FULL_D.MIG_FIN.POST_STACK.3D.JS-017536.segy ../data/seismicinversion
# TO BE RUN ONLY ONCE TO RETRIEVE THE DATA
#!az storage blob download --account-name datavillagesa --container-name volve --name Seismic/ST10010/Velocities/ST10010ZC11-MIG-VEL.MIG_VEL.VELOCITY.3D.JS-017527.segy --file ST10010ZC11-MIG-VEL.MIG_VEL.VELOCITY.3D.JS-017527.segy --sas-token "$YOURTOKEN"
#!mv ST10010ZC11-MIG-VEL.MIG_VEL.VELOCITY.3D.JS-017527.segy ../data/seismicinversion

Once the data is downloaded on your machine you are ready for the fun part of this notebook. More specifically, we will perform the following step:

  • Data and velocity model are read from SEG-Y file using segyio (note that for the Volve data we will have to deal with irregular geometry);
  • Velocity model is resampled to the data grid (both time and spatial sampling) and scaled to become the background AI for inversion;
  • Absolute inversion is applied by means of pylops.avo.poststack.PoststackInversion;
  • TV regularized inversion is performed using PyProximal.

Data loading

Let's first import all the libraries we need

# Run this when using Colab (will install the missing libraries)
# !pip install pylops pyproximal segyio scooby
%load_ext autoreload
%autoreload 2
%matplotlib inline

import os
import sys
import numpy as np
import segyio
import pylops
import matplotlib.pyplot as plt

from scipy.interpolate import RegularGridInterpolator
from pylops.basicoperators import *
from pylops.avo.poststack import *
from pyproximal.proximal import *
from pyproximal import ProxOperator
from pyproximal.optimization.primal import *
from pyproximal.optimization.primaldual import *

And define some of the input parameters we will use later (refer to the following for a more detailed description of the types of inversion)

# Data
itmin = 600 # index of first time/depth sample in data used in colored inversion
itmax = 800 # index of last time/depth sample in data used in colored inversion

# Subsampling (can save file at the end only without subsampling)
jt = 1
jil = 2
jxl = 2

# Wavelet estimation
nt_wav = 21 # number of samples of statistical wavelet
nfft = 512 # number of samples of FFT

# Trace-by-Trace Inversion
epsI_tt = 1e-3 # damping

# Spatially simultaneous
niter_sr = 3 # number of iterations of lsqr
epsI_sr = 1e-4 # damping
epsR_sr = 1e2 # spatial regularization

# Blocky simultaneous
niter_out_b = 3 # number of outer loop iterations
niter_in_b = 1 # number of inner loop iterations
niter_b = 10 # number of iterations of lsqr
mu_b = 1e-1 # damping for data term
epsI_b = 1e-4 # damping
epsR_b = 0.1 # spatial regularization
epsRL1_b = 1. # blocky regularization

Let's now read the Volve data.

Note that we add the ignore_geometry=True parameter when we open the file. As we will see the geometry in this file is not regular, so we cannot rely on the inner working of segyio to get our data into a 3d numpy.

We thus need to read all ILs, XLs and traces and reorganize them ourselves into a regular grid. No worries, numpy will do most of the hard work.

segyfile = 'data/ST10010ZC11_PZ_PSDM_KIRCH_FULL_D.MIG_FIN.POST_STACK.3D.JS-017536.segy'
f = segyio.open(segyfile, ignore_geometry=True)

traces = segyio.collect(f.trace)[:]
traces = traces[:, itmin:itmax]
ntraces, nt = traces.shape

t = f.samples[itmin:itmax]
il = f.attributes(segyio.TraceField.INLINE_3D)[:]
xl = f.attributes(segyio.TraceField.CROSSLINE_3D)[:]

# Define regular IL and XL axes
il_unique = np.unique(il)
xl_unique = np.unique(xl)

il_min, il_max = min(il_unique), max(il_unique)
xl_min, xl_max = min(xl_unique), max(xl_unique)

dt = t[1] - t[0]
dil = min(np.unique(np.diff(il_unique)))
dxl = min(np.unique(np.diff(xl_unique)))

ilines = np.arange(il_min, il_max + dil, dil)
xlines = np.arange(xl_min, xl_max + dxl, dxl)
nil, nxl = ilines.size, xlines.size

ilgrid, xlgrid = np.meshgrid(np.arange(nil),
                             np.arange(nxl),
                             indexing='ij')

# Look-up table
traces_indeces = np.full((nil, nxl), np.nan)
iils = (il - il_min) // dil
ixls = (xl - xl_min) // dxl
traces_indeces[iils, ixls] = np.arange(ntraces)
traces_available = np.logical_not(np.isnan(traces_indeces))

# Reorganize traces in regular grid
d = np.zeros((nil, nxl, nt))
d[ilgrid.ravel()[traces_available.ravel()],
  xlgrid.ravel()[traces_available.ravel()]] = traces

# Subsample
d = d[::jil, ::jxl, ::jt]
ilines = ilines[::jil]
xlines = xlines[::jxl]
t = t[::jt]
nil, nxl, nt = len(ilines), len(xlines), len(t)

# Display data
plt.figure(figsize=(12, 6))
plt.imshow(d[nil//2].T, cmap='RdYlBu', vmin=-5, vmax=5,
           extent=(xlines[0], xlines[-1], t[-1], t[0]))
plt.title('Seismic data - section')
plt.colorbar()
plt.axis('tight')
(1961.0, 2679.0, 3200.0, 2404.0)
<Figure size 864x432 with 2 Axes>

We read also the migration velocity model. In this case, the SEG-Y file is in a regular grid, but the grid is different from that of the data.

Let's resample the velocity model to the grid of the data.

segyfilev = 'data/ST10010ZC11-MIG-VEL.MIG_VEL.VELOCITY.3D.JS-017527.segy'
fv = segyio.open(segyfilev)
v = segyio.cube(fv)

IL, XL, T = np.meshgrid(ilines, xlines, t, indexing='ij')

vinterp = RegularGridInterpolator((fv.ilines, fv.xlines, fv.samples), v, 
                                  bounds_error=False, fill_value=0)
vinterp = vinterp(np.vstack((IL.ravel(), XL.ravel(), T.ravel())).T)
vinterp = vinterp.reshape(nil, nxl, nt)
# Display data
plt.figure(figsize=(12, 6))
plt.imshow(v[len(fv.ilines)//2].T, cmap='terrain',
           extent=(fv.xlines[0], fv.xlines[-1], fv.samples[-1], fv.samples[0]))
plt.title('Velocity - section')
plt.colorbar()
plt.axis('tight')

# Display data
plt.figure(figsize=(12, 6))
plt.imshow(d[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
           extent=(xlines[0], xlines[-1], t[-1], t[0]))
plt.title('Data - section')
plt.colorbar()
plt.axis('tight');

# Display data
plt.figure(figsize=(12, 6))
plt.imshow(vinterp[nil//2].T, cmap='terrain',
           extent=(xlines[0], xlines[-1], t[-1], t[0]))
plt.title('Regridded velocity - section')
plt.colorbar()
plt.axis('tight');
<Figure size 864x432 with 2 Axes><Figure size 864x432 with 2 Axes><Figure size 864x432 with 2 Axes>

Velocity model preparation

We need now to scale this model to its acoustic impedance equivalent.

This calibration step was performed outside of this notebook using a welllog and this velocity model along a well trajectory. In this example we will simply use a scaling (gradient) and a shift (intercept) from that study.

intercept = -3218.0003362662665
gradient = 3.2468122679241023

aiinterp = intercept + gradient*vinterp

# Display data
plt.figure(figsize=(12, 6))
plt.imshow(aiinterp[nil//2].T, cmap='terrain',
           extent=(xlines[0], xlines[-1], t[-1], t[0]))
plt.title('Regridded AI - section')
plt.colorbar()
plt.axis('tight');
<Figure size 864x432 with 2 Axes>

Statistical wavelet estimation

Let's now try to get a quick estimate of the wavelet in our data using a simple statistical wavelet estimation in frequency domain.

Note that this notebook is not focused on the pre-processing but we will need access to this to apply a colored inversion.

# Wavelet time axis
t_wav = np.arange(nt_wav) * (dt/1000)
t_wav = np.concatenate((np.flipud(-t_wav[1:]), t_wav), axis=0)

# Estimate wavelet spectrum
wav_est_fft = np.mean(np.abs(np.fft.fft(d[::2, ::2], nfft, axis=-1)), axis=(0, 1))
fwest = np.fft.fftfreq(nfft, d=dt/1000)

# Create wavelet in time
wav_est = np.real(np.fft.ifft(wav_est_fft)[:nt_wav])
wav_est = np.concatenate((np.flipud(wav_est[1:]), wav_est), axis=0)
wav_est = wav_est / wav_est.max()

# Display wavelet
fig, axs = plt.subplots(1, 2, figsize=(20, 5))
fig.suptitle('Statistical wavelet estimate')
axs[0].plot(fwest[:nfft//2], wav_est_fft[:nfft//2], 'k')
axs[0].set_title('Frequency')
axs[1].plot(t_wav, wav_est, 'k')
axs[1].set_title('Time');
<Figure size 1440x360 with 2 Axes>

2D Inversion

We can apply seismic inversion using a starting background AI model. The inverted AI model will therefore have the correct physical quantities of acoustic impedance.

# Swap time axis back to first dimension
d2d = d[nil//2].T
aiinterp2d = aiinterp[nil//2].T

m0 = np.log(aiinterp2d)
m0[np.isnan(m0)] = 0

Inversion with least-squares solver

wav_amp = 1e1 # guessed as we have estimated the wavelet statistically
mls, rls = \
    pylops.avo.poststack.PoststackInversion(d2d, wav_amp*wav_est, m0=m0, explicit=False, 
                                            epsR=epsR_b,
                                            **dict(show=True, iter_lim=niter_b, damp=epsI_b))
mls = np.exp(mls).T
 
LSQR            Least-squares solution of  Ax = b
The matrix A has 144000 rows and 72000 columns
damp = 1.00000000000000e-04   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       10
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   4.781e+02  4.781e+02    1.0e+00  3.2e-02
     1 -9.98364e-04   2.124e+02  2.124e+02    4.4e-01  5.4e-01   1.7e+01  1.0e+00
     2 -1.75151e-03   1.450e+02  1.450e+02    3.0e-01  3.5e-01   2.2e+01  2.2e+00
     3 -2.21677e-02   1.103e+02  1.103e+02    2.3e-01  2.6e-01   2.6e+01  3.6e+00
     4 -4.91789e-02   8.817e+01  8.817e+01    1.8e-01  2.1e-01   2.9e+01  5.3e+00
     5 -7.23960e-02   7.337e+01  7.337e+01    1.5e-01  1.7e-01   3.2e+01  7.0e+00
     6 -6.59324e-02   6.330e+01  6.330e+01    1.3e-01  1.5e-01   3.5e+01  8.9e+00
     7 -5.25435e-02   5.494e+01  5.494e+01    1.1e-01  1.3e-01   3.8e+01  1.1e+01
     8 -3.74643e-02   4.824e+01  4.824e+01    1.0e-01  1.2e-01   4.0e+01  1.3e+01
     9 -2.70929e-02   4.252e+01  4.252e+01    8.9e-02  1.1e-01   4.2e+01  1.6e+01
    10 -1.60231e-02   3.796e+01  3.796e+01    7.9e-02  1.0e-01   4.4e+01  1.8e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 3.8e+01   anorm = 4.4e+01   arnorm = 1.7e+02
itn   =      10   r2norm = 3.8e+01   acond = 1.8e+01   xnorm  = 4.8e+01
 

# Visualize
fig, axs = plt.subplots(5, 1, figsize=(12, 25))
fig.suptitle('Least-squares inversion - iline section',
             y=0.91, fontweight='bold', fontsize=18)
axs[0].imshow(d2d, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[0].set_title('Seismic data')
axs[0].axis('tight')
axs[1].imshow(rls, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[1].set_title('Residual')
axs[1].axis('tight')
axs[2].imshow(aiinterp2d, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[2].set_title('Background AI')
axs[2].axis('tight')
axs[3].imshow(mls.T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[3].set_title('Inverted AI')
axs[3].axis('tight');
axs[4].imshow(mls.T - aiinterp2d, cmap='seismic',
              vmin=-0.7*(mls.T-aiinterp2d).max(), vmax=0.7*(mls.T-aiinterp2d).max(),
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[4].set_title('Inversion dAI')
axs[4].axis('tight');
<Figure size 864x1800 with 5 Axes>

Inversion with Split-Bregman

msb, rsb = \
    pylops.avo.poststack.PoststackInversion(d2d, wav_amp*wav_est, m0=m0, explicit=False, 
                                            epsR=epsR_b, epsRL1=epsRL1_b,
                                            **dict(mu=mu_b, niter_outer=10, 
                                                   niter_inner=1, show=True,
                                                   iter_lim=40, damp=epsI_b))
msb = np.exp(msb).T
Split-Bregman optimization
---------------------------------------------------------
The Operator Op has 72000 rows and 72000 cols
niter_outer =  10     niter_inner =   1   tol = 1.00e-10
mu = 1.00e-01         epsL1 = [1.0]	  epsL2 = [0.1]     
---------------------------------------------------------

   Itn          x[0]           r2norm          r12norm
     1   9.09400e+00        1.668e+01        2.434e+03
     2   9.33999e+00        2.761e+01        2.328e+03
     3   9.28567e+00        4.375e+01        2.264e+03
     4   9.30212e+00        6.070e+01        2.224e+03
     5   9.30282e+00        7.773e+01        2.193e+03
     6   9.30234e+00        9.419e+01        2.168e+03
     7   9.30080e+00        1.095e+02        2.147e+03
     8   9.29695e+00        1.229e+02        2.134e+03
     9   9.29200e+00        1.345e+02        2.119e+03
    10   9.28627e+00        1.446e+02        2.100e+03

Iterations = 10        Total time (s) = 1.67
---------------------------------------------------------


# Visualize
fig, axs = plt.subplots(5, 1, figsize=(12, 25))
fig.suptitle('Split-Bregman inversion - iline section',
             y=0.91, fontweight='bold', fontsize=18)
axs[0].imshow(d2d, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[0].set_title('Seismic data')
axs[0].axis('tight')
axs[1].imshow(rsb, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[1].set_title('Residual')
axs[1].axis('tight')
axs[2].imshow(aiinterp2d, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[2].set_title('Background AI')
axs[2].axis('tight')
axs[3].imshow(msb.T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[3].set_title('Inverted AI')
axs[3].axis('tight');
axs[4].imshow(msb.T - aiinterp2d, cmap='seismic',
              vmin=-0.7*(mls.T-aiinterp2d).max(), vmax=0.7*(mls.T-aiinterp2d).max(),
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[4].set_title('Inversion dAI')
axs[4].axis('tight');
<Figure size 864x1800 with 5 Axes>

Inversion with Primal-Dual

# Modelling operator
Lop = PoststackLinearModelling(wav_amp*wav_est, nt0=nt, spatdims=nxl)
l2 = L2(Op=Lop, b=d2d.ravel(), niter=70, warm=True)

# Regularization
sigma = 0.5
l1 = L21(ndim=2, sigma=sigma)
Dop = Gradient(dims=(nt, nxl), edge=True, dtype=Lop.dtype, kind='forward')

# Steps 
L = 8. # np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1. / np.sqrt(L)
mu = 0.99 / (tau * L)

mpd = PrimalDual(l2, l1, Dop, m0.ravel(), tau=tau, mu=mu, theta=1., niter=50, show=True)
rpd = d2d.ravel() - Lop * mpd

mpd = np.exp(mpd).T
mpd = mpd.reshape(aiinterp2d.shape)
rpd = rpd.reshape(d2d.shape)
Primal-dual: min_x f(Ax) + x^T z + g(x)
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L21.L21'>
Linear operator (A): <class 'pylops.basicoperators.VStack.VStack'>
Additional vector (z): None
tau = 3.535534e-01	mu = 3.500179e-01
theta = 1.00		niter = 50

   Itn       x[0]          f           g          z^x       J = f + g + z^x
     1   8.49050e+00   1.102e+03   2.015e+03   0.000e+00       3.117e+03
     2   8.56054e+00   2.794e+02   1.880e+03   0.000e+00       2.160e+03
     3   8.61955e+00   1.782e+02   1.834e+03   0.000e+00       2.013e+03
     4   8.67312e+00   1.458e+02   1.806e+03   0.000e+00       1.952e+03
     5   8.72681e+00   1.340e+02   1.785e+03   0.000e+00       1.919e+03
     6   8.78197e+00   1.302e+02   1.765e+03   0.000e+00       1.895e+03
     7   8.83552e+00   1.298e+02   1.747e+03   0.000e+00       1.876e+03
     8   8.88222e+00   1.309e+02   1.729e+03   0.000e+00       1.860e+03
     9   8.91678e+00   1.326e+02   1.715e+03   0.000e+00       1.848e+03
    10   8.93686e+00   1.346e+02   1.703e+03   0.000e+00       1.838e+03
    11   8.94368e+00   1.367e+02   1.694e+03   0.000e+00       1.830e+03
    16   8.91421e+00   1.440e+02   1.657e+03   0.000e+00       1.801e+03
    21   8.91309e+00   1.470e+02   1.628e+03   0.000e+00       1.775e+03
    26   8.90161e+00   1.482e+02   1.608e+03   0.000e+00       1.756e+03
    31   8.89860e+00   1.488e+02   1.593e+03   0.000e+00       1.742e+03
    36   8.90426e+00   1.490e+02   1.582e+03   0.000e+00       1.731e+03
    41   8.91364e+00   1.491e+02   1.573e+03   0.000e+00       1.722e+03
    42   8.91532e+00   1.491e+02   1.571e+03   0.000e+00       1.721e+03
    43   8.91678e+00   1.491e+02   1.570e+03   0.000e+00       1.719e+03
    44   8.91795e+00   1.492e+02   1.569e+03   0.000e+00       1.718e+03
    45   8.91879e+00   1.492e+02   1.567e+03   0.000e+00       1.717e+03
    46   8.91928e+00   1.492e+02   1.566e+03   0.000e+00       1.715e+03
    47   8.91948e+00   1.492e+02   1.565e+03   0.000e+00       1.714e+03
    48   8.91946e+00   1.492e+02   1.564e+03   0.000e+00       1.713e+03
    49   8.91921e+00   1.493e+02   1.563e+03   0.000e+00       1.712e+03
    50   8.91875e+00   1.493e+02   1.562e+03   0.000e+00       1.711e+03

Total time (s) = 16.62
---------------------------------------------------------


# Visualize
fig, axs = plt.subplots(5, 1, figsize=(12, 25))
fig.suptitle('Primal-Dual inversion - iline section',
             y=0.91, fontweight='bold', fontsize=18)
axs[0].imshow(d2d, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[0].set_title('Seismic data')
axs[0].axis('tight')
axs[1].imshow(rpd, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[1].set_title('Residual')
axs[1].axis('tight')
axs[2].imshow(aiinterp2d, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[2].set_title('Background AI')
axs[2].axis('tight')
axs[3].imshow(mpd, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[3].set_title('Inverted AI')
axs[3].axis('tight');
axs[4].imshow(mpd - aiinterp2d, cmap='seismic',
              vmin=-0.7*(mls-aiinterp2d.T).max(), vmax=0.7*(mls-aiinterp2d.T).max(),
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[4].set_title('Inversion dAI')
axs[4].axis('tight');
<Figure size 864x1800 with 5 Axes>

Let's now visualize a single trace for the three inversions

plt.figure(figsize=(3, 12))
plt.plot(aiinterp[nil//2,nxl//2], t, 'k', lw=2, label='back')
plt.plot(mls[nxl//2], t, 'r', lw=2, label='LS')
plt.plot(msb[nxl//2], t, 'g', lw=2, label='Split-Bregman')
plt.plot(mpd[:, nxl//2], t, 'b', lw=2, label='PD')
plt.gca().invert_yaxis()
plt.legend();
<Figure size 216x864 with 1 Axes>

3D Inversion

Inversion with least-squares solver

# Swap time axis back to first dimension
d = np.swapaxes(d, -1, 0)
aiinterp = np.swapaxes(aiinterp, -1, 0)

m0 = np.log(aiinterp)
m0[np.isnan(m0)] = 0

# Inversion
wav_amp = 1e1 # guessed as we have estimated the wavelet statistically
minv, rinv = \
    pylops.avo.poststack.PoststackInversion(d, wav_amp*wav_est, m0=m0, explicit=False, 
                                            epsR=epsR_b,
                                            **dict(show=True, iter_lim=niter_b, damp=epsI_b))
#minv, rinv = \
#    pylops.avo.poststack.PoststackInversion(d, wav_amp*wav_est, m0=m0, explicit=False, 
#                                            epsR=epsR_b, epsRL1=epsRL1_b,
#                                            **dict(mu=mu_b, niter_outer=niter_out_b, 
#                                                   niter_inner=niter_in_b, show=True,
#                                                   iter_lim=niter_b, damp=epsI_b))

# Swap time axis back to last dimension
aiinterp = np.swapaxes(aiinterp, 0, -1)
d = np.swapaxes(d, 0, -1)
minv = np.swapaxes(np.exp(minv), 0, -1)
rinv = np.swapaxes(rinv, 0, -1)
<ipython-input-21-780e0c8ab0cb>:5: RuntimeWarning: invalid value encountered in log
  m0 = np.log(aiinterp)

 
LSQR            Least-squares solution of  Ax = b
The matrix A has 28944000 rows and 14472000 columns
damp = 1.00000000000000e-04   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       10
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   6.941e+03  6.941e+03    1.0e+00  2.3e-03
     1  2.91663e-04   2.948e+03  2.948e+03    4.2e-01  5.5e-01   1.7e+01  1.0e+00
     2  8.48306e-04   2.004e+03  2.004e+03    2.9e-01  3.5e-01   2.2e+01  2.2e+00
     3  1.68877e-03   1.518e+03  1.518e+03    2.2e-01  2.5e-01   2.6e+01  3.6e+00
     4  2.73650e-03   1.240e+03  1.240e+03    1.8e-01  2.0e-01   2.9e+01  5.2e+00
     5  4.07560e-03   1.049e+03  1.049e+03    1.5e-01  1.7e-01   3.2e+01  7.0e+00
     6  5.66421e-03   9.081e+02  9.081e+02    1.3e-01  1.4e-01   3.5e+01  9.0e+00
     7  7.50241e-03   8.009e+02  8.009e+02    1.2e-01  1.2e-01   3.8e+01  1.1e+01
     8  9.50411e-03   7.207e+02  7.207e+02    1.0e-01  1.1e-01   4.0e+01  1.3e+01
     9  1.17558e-02   6.568e+02  6.568e+02    9.5e-02  9.4e-02   4.2e+01  1.6e+01
    10  1.42242e-02   6.058e+02  6.058e+02    8.7e-02  8.5e-02   4.4e+01  1.8e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 6.1e+02   anorm = 4.4e+01   arnorm = 2.3e+03
itn   =      10   r2norm = 6.1e+02   acond = 1.8e+01   xnorm  = 6.6e+02
 

# Visualize
fig, axs = plt.subplots(5, 1, figsize=(12, 25))
fig.suptitle('Least-squares inversion - iline section',
             y=0.91, fontweight='bold', fontsize=18)
axs[0].imshow(d[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[0].set_title('Seismic data')
axs[0].axis('tight')
axs[1].imshow(rinv[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[1].set_title('Residual')
axs[1].axis('tight')
axs[2].imshow(aiinterp[nil//2].T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[2].set_title('Background AI')
axs[2].axis('tight')
axs[3].imshow(minv[nil//2].T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[3].set_title('Inverted AI')
axs[3].axis('tight');
axs[4].imshow(minv[nil//2].T -  aiinterp[nil//2].T, cmap='seismic',
              vmin=-0.7*(minv-aiinterp).max(), vmax=0.7*(minv-aiinterp).max(),
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[4].set_title('Inversion dAI')
axs[4].axis('tight')

fig, axs = plt.subplots(5, 1, figsize=(12, 25))
fig.suptitle('Least-squares inversion - time section',
             y=0.91, fontweight='bold', fontsize=18)
axs[0].imshow(d[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
              extent=(ilines[0], xlines[-1], xlines[0], xlines[-1]))
axs[0].set_title('Seismic data')
axs[0].axis('tight')
axs[1].imshow(rinv[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
              extent=(ilines[0], xlines[-1], xlines[0], xlines[-1]))
axs[1].set_title('Residual')
axs[1].axis('tight')
axs[2].imshow(aiinterp[nil//2].T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(ilines[0], xlines[-1], xlines[0], xlines[-1]))
axs[2].set_title('Background AI')
axs[2].axis('tight')
axs[3].imshow(minv[..., nt//2], cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(ilines[0], xlines[-1], xlines[0], xlines[-1]))
axs[3].set_title('Inverted AI')
axs[3].axis('tight');
axs[4].imshow(minv[nil//2].T -  aiinterp[nil//2].T, cmap='seismic',
              vmin=-0.7*(minv-aiinterp).max(), vmax=0.7*(minv-aiinterp).max(),
              extent=(ilines[0], xlines[-1], xlines[0], xlines[-1]))
axs[4].set_title('Inversion dAI')
axs[4].axis('tight');
<Figure size 864x1800 with 5 Axes><Figure size 864x1800 with 5 Axes>

Inversion with Primal-Dual

# Swap time axis back to first dimension
d = np.swapaxes(d, -1, 0)
aiinterp = np.swapaxes(aiinterp, -1, 0)

m0 = np.log(aiinterp)
m0[np.isnan(m0)] = 0
print(np.max(m0))

# Modeling operator
Lop = PoststackLinearModelling(wav_amp*wav_est, nt0=nt, spatdims=(nxl, nil))
l2 = L2(Op=Lop, b=d.ravel(), niter=70, warm=True)

# Regularization
sigma = 0.2
l1 = L21(ndim=3, sigma=sigma)
Dop = Gradient(dims=(nt, nxl, nil), edge=True, dtype=Lop.dtype, kind='forward')

# Steps
L = 12. # np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1. / np.sqrt(L)
mu = 0.99 / (tau * L)

mpd = PrimalDual(l2, l1, Dop, m0.ravel(), tau=tau, mu=mu, theta=1., niter=10, show=True)
rpd = d.ravel() - Lop * mpd

mpd = mpd.reshape(aiinterp.shape)
rpd = rpd.reshape(d.shape)

# Swap time axis back to last dimension
aiinterp = np.swapaxes(aiinterp, 0, -1)
d = np.swapaxes(d, 0, -1)
rpd = np.swapaxes(rpd, 0, -1)
mpd = np.swapaxes(np.exp(mpd), 0, -1)
<ipython-input-23-5364b53a5b0a>:5: RuntimeWarning: invalid value encountered in log
  m0 = np.log(aiinterp)

9.55030812870091
Primal-dual: min_x f(Ax) + x^T z + g(x)
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L21.L21'>
Linear operator (A): <class 'pylops.basicoperators.VStack.VStack'>
Additional vector (z): None
tau = 2.886751e-01	mu = 2.857884e-01
theta = 1.00		niter = 10

   Itn       x[0]          f           g          z^x       J = f + g + z^x
     1   5.77063e-02   1.836e+05   3.027e+05   0.000e+00       4.863e+05
     2   8.85981e-02   6.492e+04   3.044e+05   0.000e+00       3.693e+05
     3   1.27763e-01   4.572e+04   3.022e+05   0.000e+00       3.479e+05
     4   1.71610e-01   3.751e+04   2.992e+05   0.000e+00       3.367e+05
     5   2.17243e-01   3.316e+04   2.963e+05   0.000e+00       3.294e+05
     6   2.64567e-01   3.062e+04   2.937e+05   0.000e+00       3.244e+05
     7   3.13475e-01   2.898e+04   2.916e+05   0.000e+00       3.206e+05
     8   3.63843e-01   2.782e+04   2.899e+05   0.000e+00       3.177e+05
     9   4.15464e-01   2.692e+04   2.883e+05   0.000e+00       3.152e+05
    10   4.68104e-01   2.618e+04   2.868e+05   0.000e+00       3.130e+05

Total time (s) = 1759.56
---------------------------------------------------------


# Visualize
fig, axs = plt.subplots(5, 1, figsize=(12, 25))
fig.suptitle('Primal-Dual inversion - iline section',
             y=0.91, fontweight='bold', fontsize=18)
axs[0].imshow(d[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[0].set_title('Seismic data')
axs[0].axis('tight')
axs[1].imshow(rpd[nil//2].T, cmap='seismic', vmin=-10, vmax=10,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[1].set_title('Residual')
axs[1].axis('tight')
axs[2].imshow(aiinterp[nil//2].T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[2].set_title('Background AI')
axs[2].axis('tight')
axs[3].imshow(mpd[nil//2].T, cmap='terrain',
              vmin=6000, vmax=18000,
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[3].set_title('Inverted AI')
axs[3].axis('tight');
axs[4].imshow(mpd[nil//2].T-aiinterp[nil//2].T, cmap='seismic',
              vmin=-0.7*(minv-aiinterp).max(), vmax=0.7*(minv-aiinterp).max(),
              extent=(xlines[0], xlines[-1], t[-1], t[0]))
axs[4].set_title('Inversion dAI')
axs[4].axis('tight');
<Figure size 864x1800 with 5 Axes>
Notebooks
Horizons