PyLops

This Python library is inspired by the MATLAB Spot – A Linear-Operator Toolbox project.

Linear operators and inverse problems are at the core of many of the most used algorithms in signal processing, image processing, and remote sensing. When dealing with small-scale problems, the Python numerical scientific libraries numpy and scipy allow to perform most of the underlying matrix operations (e.g., computation of matrix-vector products and manipulation of matrices) in a simple and expressive way.

Many useful operators, however, do not lend themselves to an explicit matrix representation when used to solve large-scale problems. PyLops operators, on the other hand, still represent a matrix and can be treated in a similar way, but do not rely on the explicit creation of a dense (or sparse) matrix itself. Conversely, the forward and adjoint operators are represented by small pieces of codes that mimic the effect of the matrix on a vector or another matrix.

Luckily, many iterative methods (e.g. cg, lsqr) do not need to know the individual entries of a matrix to solve a linear system. Such solvers only require the computation of forward and adjoint matrix-vector products as done for any of the PyLops operators.

Here is a simple example showing how a dense first-order first derivative operator can be created, applied and inverted using numpy/scipy commands:

import numpy as np
from scipy.linalg import lstsq

nx = 7
x = np.arange(nx) - (nx-1)/2

D = np.diag(0.5*np.ones(nx-1), k=1) - \
    np.diag(0.5*np.ones(nx-1), k=-1)
D[0] = D[-1] = 0 # take away edge effects

# y = Dx
y = np.dot(D, x)
# x = D'y
xadj = np.dot(D.T, y)
# xinv = D^-1 y
xinv = lstsq(D, y)[0]

and similarly using PyLops commands:

from pylops import FirstDerivative

nx = 7
x = range(-(nx // 2), nx // 2 + (1 if nx % 2 else 0))

Dlop = FirstDerivative(nx, dtype='float64')

# y = Dx
y = Dlop*x
# x = D'y
xadj = Dlop.H*y
# xinv = D^-1 y
xinv = Dlop / y

Note how this second approach does not require creating a dense matrix, reducing both the memory load and the computational cost of applying a derivative to an input vector \(\mathbf{x}\). Moreover, the code becomes even more compact and espressive than in the previous case letting the user focus on the formulation of equations of the forward problem to be solved by inversion.

Terminology

A common terminology is used within the entire documentation of PyLops. Every linear operator and its application to a model will be referred to as forward model (or operation)

\[\mathbf{y} = \mathbf{A} \mathbf{x}\]

while its application to a data is referred to as adjoint modelling (or operation)

\[\mathbf{x} = \mathbf{A}^H \mathbf{y}\]

where \(\mathbf{x}\) is called model and \(\mathbf{y}\) is called data. The operator \(\mathbf{A}:\mathbb{F}^m \to \mathbb{F}^n\) effectively maps a vector of size \(m\) in the model space to a vector of size \(n\) in the data space, conversely the adjoint operator \(\mathbf{A}^H:\mathbb{F}^n \to \mathbb{F}^m\) maps a vector of size \(n\) in the data space to a vector of size \(m\) in the model space. As linear operators mimics the effect a matrix on a vector we can also loosely refer to \(m\) as the number of columns and \(n\) as the number of rows of the operator.

Ultimately, solving an inverse problems accounts to removing the effect of \(\mathbf{A}\) from the data \(\mathbf{y}\) to retrieve the model \(\mathbf{x}\).

For a more detailed description of the concepts of linear operators, adjoints and inverse problems in general, you can head over to one of Jon Claerbout’s books such as Basic Earth Imaging.

Implementation

PyLops is build on top of the scipy class scipy.sparse.linalg.LinearOperator.

This class allows in fact for the creation of objects (or interfaces) for matrix-vector and matrix-matrix products that can ultimately be used to solve any inverse problem of the form \(\mathbf{y}=\mathbf{A}\mathbf{x}\).

As explained in the scipy LinearOperator official documentation, to construct a scipy.sparse.linalg.LinearOperator, a user is required to pass appropriate callables to the constructor of this class, or subclass it. More specifically one of the methods _matvec and _matmat must be implemented for the forward operator and one of the methods _rmatvec or _adjoint may be implemented to apply the Hermitian adjoint. The attributes/properties shape (pair of integers) and dtype (may be None) must also be provided during __init__ of this class.

Any linear operator developed within the PyLops library follows this philosophy. As explained more in details in Implementing new operators section, a linear operator is created by subclassing the scipy.sparse.linalg.LinearOperator class and _matvec and _rmatvec are implemented.

History

PyLops was initially written by Equinor It is a flexible and scalable python library for large-scale optimization with linear operators that can be tailored to our needs, and as contribution to the free software community. Since June 2021, PyLops is a NUMFOCUS Affiliated Project.

Installation

The PyLops project strives to create a library that is easy to install in any environment and has a very limited number of dependencies. However, since Python2 will retire soon, we have decided to only focus on a Python3 implementation. If you are still using Python2, hurry up!

For this reason you will need Python 3.6 or greater to get started.

Dependencies

Our mandatory dependencies are limited to:

• numpy
• scipy

We advise using the Anaconda Python distribution to ensure that these dependencies are installed via the Conda package manager. This is not just a pure stylistic choice but comes with some hidden advantages, such as the linking to Intel MKL library (i.e., a highly optimized BLAS library created by Intel).

If you simply want to use PyLops for teaching purposes or for small-scale examples, this should not really affect you. However, if you are interested in getting better code performance, read carefully the Advanced installation page.

Optional dependencies

PyLops’s optional dependencies refer to those dependencies that we do not include in our requirements.txt and environment.yml files and thus are not strictly needed nor installed directly as part of a standar installation (see below for details)

However, we sometimes implement additional back-ends (referred to as engine in the code) for some of our operators in order to improve their performance. To do so, we rely on third-party libraries. Those libraries are generally added to the list of our optional dependencies. If you are not after code performance, you may simply stick to the mandatory dependencies and pylops will ensure to always fallback to one of those for any linear operator.

If you are instead after code performance, take a look at the Optional dependencies section in the Advanced installation page.

Step-by-step installation for users

Python environment

Activate your Python environment, and simply type the following command in your terminal to install the PyPi distribution:

>> pip install pylops

If using Conda, you can also install our conda-forge distribution via:

>> conda install -c conda-forge pylops

Note that using the conda-forge distribution is recommended as all the dependencies (both mandatory and optional) will be correctly installed for you, while only mandatory dependencies are installed using the pip distribution.

Alternatively, to access the latest source from github:

>> pip install git+https://git@github.com/PyLops/pylops.git@master

or just clone the repository

>> git clone https://github.com/PyLops/pylops.git

or download the zip file from the repository (green button in the top right corner of the main github repo page) and install PyLops from terminal using the command:

>> make install

Docker

If you want to try PyLops but do not have Python in your local machine, you can use our Docker image instead.

After installing Docker in your computer, type the following command in your terminal (note that this will take some time the first time you type it as you will download and install the docker image):

>> docker run -it -v /path/to/local/folder:/home/jupyter/notebook -p 8888:8888 mrava87/pylops:notebook

This will give you an address that you can put in your browser and will open a jupyter-notebook enviroment with PyLops and other basic Python libraries installed. Here /path/to/local/folder is the absolute path of a local folder on your computer where you will create a notebook (or containing notebooks that you want to continue working on). Note that anything you do to the notebook(s) will be saved in your local folder.

A larger image with Conda distribution is also available. Simply use conda_notebook instead of notebook in the previous command.

Step-by-step installation for developers

Fork and clone the repository by executing the following in your terminal:

>> git clone https://github.com/your_name_here/pylops.git

The first time you clone the repository run the following command:

>> make dev-install

If you prefer to build a new Conda enviroment just for PyLops, run the following command:

>> make dev-install_conda

To ensure that everything has been setup correctly, run tests:

>> make tests

Make sure no tests fail, this guarantees that the installation has been successfull.

If using Conda environment, always remember to activate the conda environment every time you open a new bash shell by typing:

>> source activate pylops

Advanced installation

In this section we discuss some important details regarding code performance when using PyLops.

To get the most out of PyLops operators in terms of speed you will need to follow these guidelines as much as possible or ensure that the Python libraries used by PyLops are efficiently installed (e.g., allow multithreading) in your systemt.

Dependencies

PyLops relies on the numpy and scipy libraries and being able to link these to the most performant BLAS will ensure optimal performance of PyLops when using only required dependencies.

As already mentioned in the Installation page, we strongly encourage using the Anaconda Python distribution as numpy and scipy will be automatically linked to the Intel MKL library, which is per today the most performant library for basic linear algebra operations (if you don’t believe it, take a read at this blog post).

The best way to understand which BLAS library is currently linked to your numpy and scipy libraries is to run the following commands in ipython:

import numpy as np
import scipy as sp
print(np.__config__.show())
print(sp.__config__.show())

You should be able to understand if your numpy and scipy are linked to Intel MKL or something else.

Note

Unfortunately, PyLops is so far only shipped with PyPI, meaning that if you have not already installed numpy and scipy in your environment they will be installed as part of the installation process of the pylops library, all of those using pip. This comes with the disadvantage that numpy and scipy are linked to OpenBlas instead of Intel MKL, leading to a loss of performance. To prevent this, we suggest the following strategy:

• create conda environment, e.g. conda create -n envname python=3.6.4 numpy scipy
• install pylops using pip install pylops

Finally, it is always important to make sure that your environment variable OMP_NUM_THREADS is correctly set to the maximum number of threads you would like to use in your code. If that is not the case numpy and scipy will underutilize your hardware even if linked to a performant BLAS library.

For example, first set OMP_NUM_THREADS=1 (single-threaded) in your terminal:

>> export OMP_NUM_THREADS=1

and run the following code in python:

import os
import numpy as np
from timeit import timeit

size = 4096
A = np.random.random((size, size)),
B = np.random.random((size, size))
print('Time with %s threads: %f s' \
      %(os.environ.get('OMP_NUM_THREADS'),
        timeit(lambda: np.dot(A, B), number=4)))

Subsequently set OMP_NUM_THREADS=2, or any higher number of threads available in your hardware (multi-threaded):

>> export OMP_NUM_THREADS=2

and run the same python code. By both looking at your processes (e.g. using top) and at the python print statement you should see a speed-up in the second case.

Alternatively, you could set the OMP_NUM_THREADS variable directly inside your script using os.environ['OMP_NUM_THREADS']=str(2). Moreover, note that when using Intel MKL you can alternatively set the MKL_NUM_THREADS instead of OMP_NUM_THREADS: this could be useful if your code runs other parallel processes which you can control indipendently from the Intel MKL ones using OMP_NUM_THREADS.

Note

Always remember to set OMP_NUM_THREADS (or MKL_NUM_THREADS) in your enviroment when using PyLops

Optional dependencies

To avoid increasing the number of required dependencies, which may lead to conflicts with other libraries that you have in your system, we have decided to build some of the additional features of PyLops in such a way that if an optional dependency is not present in your python environment, a safe fallback to one of the required dependencies will be enforced.

When available in your system, we reccomend using the Conda package manager and install all the mandatory and optional dependencies of PyLops at once using the command:

>> conda install -c conda-forge pylops

in this case all dependencies will be installed from their conda distributions.

Alternatively, from version 1.4.0 optional dependencies can also be installed as part of the pip installation via:

>> pip install pylops[advanced]

Dependencies are however installed from their PyPI wheels.

An exception is however represented by cupy. This library is NOT installed automatically. Users interested to accelerate their compuations with the aid of GPUs should install it prior to installing pylops (see below for more details).

numba

Although we always stive to write code for forward and adjoint operators that takes advantage of the perks of numpy and scipy (e.g., broadcasting, ufunc), in some case we may end up using for loops that may lead to poor performance. In those cases we may decide to implement alternative (optional) back-ends in numba.

In this case a user can simply switch from the native, always available implementation to the numba implementation by simply providing the following additional input parameter to the operator engine='numba'. This is for example the case in the pylops.signalprocessing.Radon2D.

If interested to use numba backend from conda, you will need to manually install it:

>> conda install numba

Finally, it is also advised to install the additional package icc_rt.

>> conda install -c numba icc_rt

or pip equivalent. Similarly to Intel MKL, you need to set the environment variable NUMBA_NUM_THREADS to tell numba how many threads to use. If this variable is not present in your environment, numba code will be compiled with parallel=False.

fft routines

Two different engines are provided by the pylops.signalprocessing.FFT operator for fft and ifft routines in the forward and adjoint modes: engine='numpy' (default) and engine='fftw'.

The first engine comes as default as numpy is part of the dependencies of PyLops and automatically installed when PyLops is installed if not already available in your Python distribution.

The second engine implements the well-known FFTW via the python wrapper pyfftw.FFTW. This optimized fft tends to outperform the one from numpy in many cases, however it has not been inserted in the mandatory requirements of PyLops, meaning that when installing PyLops with pip, pyfftw.FFTW will not be installed automatically.

Again, if interested to use FFTW backend from conda, you will need to manually install it:

>> conda install -c conda-forge pyfftw

or pip equivalent.

skfmm

This library is used to compute traveltime tables with the fast-marching method in the initialization of the pylops.waveeqprocessing.Demigration operator when choosing mode == 'eikonal'.

As this may not be of interest for many users, this library has not been inserted in the mandatory requirements of PyLops. If interested to use skfmm, you will need to manually install it:

>> conda install -c conda-forge scikit-fmm

or pip equivalent.

spgl1

This library is used to solve sparsity-promoting BP, BPDN, and LASSO problems in pylops.optimization.sparsity.SPGL1 solver.

If interested to use spgl1, you can manually install it:

>> pip install spgl1

pywt

This library is used to implement the Wavelet operators.

If interested to use pywt, you can manually install it:

>> conda install pywavelets

or pip equivalent.

Note

If you are a developer, all the above optional dependencies can also be installed automatically by cloning the repository and installing pylops via make dev-install or make dev-install_conda.

Optional Dependencies (GPU)

cupy

This library is used as a drop-in replacement to numpy for GPU-accelerated computations. Since many different versions of cupy exist (based on the CUDA drivers of the GPU), users must install cupy prior to installing pylops. PyLops will automatically check if cupy is installed and in that case use it any time the input vector passed to an operator is of cupy type. Users can however disabilitate this option even if cupy is installed. For more details of GPU-accelerated PyLops read the GPU support section.

cusignal

This library is used as a drop-in replacement to scipy.signal for GPU-accelerated computations. Similar to cupy, users must install cusignal prior to installing pylops. PyLops will automatically check if cusignal is installed and in that case use it any time the input vector passed to an operator is of cusignal type. Users can however disabilitate this option even if cupy is installed. For more details of GPU-accelerated PyLops read the GPU support section.

GPU support

From v1.12.0, PyLops supports computations on GPUs powered by cupy and cusignal.

Note

For this to work, ensure to have installed cupy-cudaXX>=8.1.0 and cusignal>=0.16.0. If neither cupy nor cusignal is installed, pylops will work just fine using its CPU backend. Note, however, that setting the environment variable CUPY_PYLOPS (or CUSIGNAL_PYLOPS) to 0 will force pylops not to use the cupy backend. This can be also used if a previous version of cupy or cusignal is installed in your system, otherwise you will get an error when importing pylops.

Apart from a few exceptions, all operators and solvers in PyLops can seamlessly work with numpy arrays on CPU as well as with cupy arrays on GPU. Users do simply need to consistently create operators and provide data vectors to the solvers - e.g., when using pylops.MatrixMult the input matrix must be a cupy array if the data provided to a solver is a cupy array.

pylops.LinearOperator methods that are currently not available for GPU computations are:

• eigs, cond, and tosparse, and estimate_spectral_norm

Operators that are currently not available for GPU computations are:

• pylops.Spread
• pylops.signalprocessing.Radon2D
• pylops.signalprocessing.Radon3D
• pylops.signalprocessing.DWT
• pylops.signalprocessing.DWT2D
• pylops.signalprocessing.Seislet
• pylops.waveeqprocessing.Demigration
• pylops.waveeqprocessing.LSM

Solvers that are currently not available for GPU computations are:

• pylops.optimization.sparsity.SPGL1

Example

Finally, let’s briefly look at an example. First we write a code snippet using numpy arrays which PyLops will run on your CPU:

ny, nx = 400, 400
G = np.random.normal(0,1,(ny, nx)).astype(np.float32)
x = np.ones(nx, dtype=np.float32)

Gop = MatrixMult(G, dtype='float32')
y = Gop * x
xest = Gop / y

Now we write a code snippet using cupy arrays which PyLops will run on your GPU:

ny, nx = 400, 400
G = cp.random.normal(0,1,(ny, nx)).astype(np.float32)
x = cp.ones(nx, dtype=np.float32)

Gop = MatrixMult(G, dtype='float32')
y = Gop * x
xest = Gop / y

The code is almost unchanged apart from the fact that we now use cupy arrays, PyLops will figure this out! For more examples head over to these notebooks.

Extensions

PyLops brings to users the power of linear operators in a simple and easy to use programming interface.

While very powerful on its own, this library is further extended to take advantage of more advanced computational resources, either in terms of multiple-node clusters or GPUs. Moreover, some indipendent small libraries are created to use third party softwares that cannot be included as dependencies to our main library for licensing issues but may be useful for academic purposes.

Spin-off projects that aim at extending the capabilities of PyLops are:

• PyLops-GPU
• PyLops-Distributed
• curvelops

Tutorials

01. The LinearOpeator

01. The LinearOpeator

This first tutorials is aimed at easing the use of the PyLops library for both new users and developers.

Since PyLops heavily relies on the use of the scipy.sparse.linalg.LinearOperator class of scipy, we will start by looking at how to initialize a linear operator as well as different ways to apply the forward and adjoint operations. Finally we will investigate various special methods, also called magic methods (i.e., methods with the double underscores at the beginning and the end) that have been implemented for such a class and will allow summing, subtractring, chaining, etc. multiple operators in very easy and expressive way.

Let’s start by defining a simple operator that applies element-wise multiplication of the model with a vector d in forward mode and element-wise multiplication of the data with the same vector d in adjoint mode. This operator is present in PyLops under the name of pylops.Diagonal and its implementation is discussed in more details in the Implementing new operators page.

import timeit
import matplotlib.pyplot as plt
import numpy as np
import pylops

n = 10
d = np.arange(n) + 1.
x = np.ones(n)
Dop = pylops.Diagonal(d)

First of all we apply the operator in the forward mode. This can be done in four different ways:

• _matvec: directly applies the method implemented for forward mode
• matvec: performs some checks before and after applying _matvec
• *: operator used to map the special method __matmul__ which checks whether the input x is a vector or matrix and applies _matvec or _matmul accordingly.
• @: operator used to map the special method __mul__ which performs like the * opetator

We will time these 4 different executions and see how using _matvec (or matvec) will result in the faster computation. It is thus advised to use * (or @) in examples when expressivity has priority but prefer _matvec (or matvec) for efficient implementations.

# setup command
cmd_setup ="""\
import numpy as np
import pylops
n = 10
d = np.arange(n) + 1.
x = np.ones(n)
Dop = pylops.Diagonal(d)
DopH = Dop.H
"""

# _matvec
cmd1 = 'Dop._matvec(x)'

# matvec
cmd2 = 'Dop.matvec(x)'

# @
cmd3 = 'Dop@x'

# *
cmd4 = 'Dop*x'

# timing
t1 = 1.e3 * np.array(timeit.repeat(cmd1, setup=cmd_setup,
                                   number=500, repeat=5))
t2 = 1.e3 * np.array(timeit.repeat(cmd2, setup=cmd_setup,
                                   number=500, repeat=5))
t3 = 1.e3 * np.array(timeit.repeat(cmd3, setup=cmd_setup,
                                   number=500, repeat=5))
t4 = 1.e3 * np.array(timeit.repeat(cmd4, setup=cmd_setup,
                                   number=500, repeat=5))

plt.figure(figsize=(7, 2))
plt.plot(t1, 'k', label=' _matvec')
plt.plot(t2, 'r', label='matvec')
plt.plot(t3, 'g', label='@')
plt.plot(t4, 'b', label='*')
plt.axis('tight')
plt.legend();

Out:

<matplotlib.legend.Legend object at 0x7fc20ad23080>

Similarly we now consider the adjoint mode. This can be done in three different ways:

• _rmatvec: directly applies the method implemented for adjoint mode
• rmatvec: performs some checks before and after applying _rmatvec
• .H*: first applies the adjoint .H which creates a new scipy.sparse.linalg._CustomLinearOperator` where _matvec and _rmatvec are swapped and then applies the new _matvec.

Once again, after timing these 3 different executions we can see see how using _rmatvec (or rmatvec) will result in the faster computation while .H* is very unefficient and slow. Note that if the adjoint has to be applied multiple times it is at least advised to create the adjoint operator by applying .H only once upfront. Not surprisingly, the linear solvers in scipy as well as in PyLops actually use matvec and rmatvec when dealing with linear operators.

# _rmatvec
cmd1 = 'Dop._rmatvec(x)'

# rmatvec
cmd2 = 'Dop.rmatvec(x)'

# .H* (pre-computed H)
cmd3 = 'DopH*x'

# .H*
cmd4 = 'Dop.H*x'

# timing
t1 = 1.e3 * np.array(timeit.repeat(cmd1, setup=cmd_setup,
                                   number=500, repeat=5))
t2 = 1.e3 * np.array(timeit.repeat(cmd2, setup=cmd_setup,
                                   number=500, repeat=5))
t3 = 1.e3 * np.array(timeit.repeat(cmd3, setup=cmd_setup,
                                   number=500, repeat=5))
t4 = 1.e3 * np.array(timeit.repeat(cmd4, setup=cmd_setup,
                                   number=500, repeat=5))

plt.figure(figsize=(7, 2))
plt.plot(t1, 'k', label=' _rmatvec')
plt.plot(t2, 'r', label='rmatvec')
plt.plot(t3, 'g', label='.H* (pre-computed H)')
plt.plot(t4, 'b', label='.H*')
plt.axis('tight')
plt.legend();

Out:

<matplotlib.legend.Legend object at 0x7fc21377bf60>

Just to reiterate once again, it is advised to call matvec and rmatvec unless PyLops linear operators are used for teaching purposes.

We now go through some other methods and special methods that are implemented in scipy.sparse.linalg.LinearOperator (and pylops.LinearOperator):

• Op1+Op2: maps the special method __add__ and performs summation between two operators and returns a pylops.LinearOperator
• -Op: maps the special method __neg__ and performs negation of an operators and returns a pylops.LinearOperator
• Op1-Op2: maps the special method __sub__ and performs summation between two operators and returns a pylops.LinearOperator
• Op1**N: maps the special method __pow__ and performs exponentiation of an operator and returns a pylops.LinearOperator
• Op/y (and Op.div(y)): maps the special method __truediv__ and performs inversion of an operator
• Op.eigs(): estimates the eigenvalues of the operator
• Op.cond(): estimates the condition number of the operator
• Op.conj(): create complex conjugate operator

# +
print(Dop + Dop)

# -
print(-Dop)
print(Dop - 0.5 * Dop)

# **
print(Dop ** 3)

#* and /
y = Dop * x
print(Dop/y)

# eigs
print(Dop.eigs(neigs=3))

# cond
print(Dop.cond())

# conj
print(Dop.conj())

Out:

<10x10 LinearOperator with dtype=float64>
<10x10 LinearOperator with dtype=float64>
<10x10 LinearOperator with dtype=float64>
<10x10 LinearOperator with dtype=float64>
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[10.+0.j  9.+0.j  8.+0.j]
(9.999999999999986+0j)
<10x10 _ConjLinearOperator with dtype=float64>

To understand the effect of conj we need to look into a problem with an operator in the complex domain. Let’s create again our pylops.Diagonal operator but this time we populate it with complex numbers. We will see that the action of the operator and its complex conjugate is different even if the model is real.

n = 5
d = 1j*(np.arange(n) + 1.)
x = np.ones(n)
Dop = pylops.Diagonal(d)

print('y = Dx = ', Dop*x)
print('y = conj(D)x = ', Dop.conj()*x)

Out:

y = Dx =  [0.+1.j 0.+2.j 0.+3.j 0.+4.j 0.+5.j]
y = conj(D)x =  [0.-1.j 0.-2.j 0.-3.j 0.-4.j 0.-5.j]

At this point, the concept of linear operator may sound abstract. The convinience method pylops.LinearOperator.todense can be used to create the equivalent dense matrix of any operator. In this case for example we expect to see a diagonal matrix with d values along the main diagonal

D = Dop.todense()

plt.figure(figsize=(5, 5))
plt.imshow(np.abs(D))
plt.title('Dense representation of Diagonal operator')
plt.axis('tight')
plt.colorbar()

Out:

<matplotlib.colorbar.Colorbar object at 0x7fc20b677780>

Finally it is worth reiterating that if two linear operators are combined by means of the algebraical operations shown above, the resulting operator is still a pylops.LinearOperator operator. This means that we can still apply any of the methods implemented in the original scipy class definition like *, as well as those in our class definition like /

Dop1 = Dop - Dop.conj()

y = Dop1 * x
print('x = (Dop - conj(Dop))/y = ', Dop1 / y)

D1 = Dop1.todense()

plt.figure(figsize=(5, 5))
plt.imshow(np.abs(D1))
plt.title(r'Dense representation of $|D + D^*|$')
plt.axis('tight')
plt.colorbar()

Out:

x = (Dop - conj(Dop))/y =  [1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j]

<matplotlib.colorbar.Colorbar object at 0x7fc20ad54a58>

This first tutorial is completed. You have seen the basic operations that can be performed using scipy.sparse.linalg.LinearOperator and our overload of such a class pylops.LinearOperator and you should be able to get started combining various PyLops operators and solving your own inverse problems.

02. The Dot-Test

02. The Dot-Test

One of the most important aspect of writing a Linear operator is to be able to verify that the code implemented in forward mode and the code implemented in adjoint mode are effectively adjoint to each other. If this is the case, your Linear operator will successfully pass the so-called dot-test. Refer to the Notes section of pylops.utils.dottest) for a more detailed description.

In this example, I will show you how to use the dot-test for a variety of operator when model and data are either real or complex numbers.

# pylint: disable=C0103
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as pltgs

import pylops
from pylops.utils import dottest

plt.close('all')

Let’s start with something very simple. We will make a pylops.MatrixMult operator and verify that its implementation passes the dot-test. For this time, we will do this step-by-step, replicating what happens in the pylops.utils.dottest routine.

N, M = 5, 3
Mat = np.arange(N*M).reshape(N, M)
Op = pylops.MatrixMult(Mat)

v = np.random.randn(N)
u = np.random.randn(M)

# Op * u
y = Op.matvec(u)
# Op'* v
x = Op.rmatvec(v)

yy = np.dot(y, v)  # (Op  * u)' * v
xx = np.dot(u, x)  # u' * (Op' * v)

print('Dot-test %e' % np.abs((yy - xx) / ((yy + xx + 1e-15) / 2)))

Out:

Dot-test 1.366060e-16

And here is a visual intepretation of what a dot-test is

gs = pltgs.GridSpec(1, 9)
fig = plt.figure(figsize=(7, 3))
ax = plt.subplot(gs[0, 0:2])
ax.imshow(Op.A, cmap='rainbow')
ax.set_title(r'$(Op*$', size=20, fontweight='bold')
ax.set_xticks(np.arange(M-1)+0.5)
ax.set_yticks(np.arange(N-1)+0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis('tight')
ax = plt.subplot(gs[0, 2])
ax.imshow(u[:, np.newaxis], cmap='rainbow')
ax.set_title(r'$u)^T$', size=20, fontweight='bold')
ax.set_xticks([])
ax.set_yticks(np.arange(M-1)+0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis('tight')
ax = plt.subplot(gs[0, 3])
ax.imshow(v[:, np.newaxis], cmap='rainbow')
ax.set_title(r'$v$', size=20, fontweight='bold')
ax.set_xticks([])
ax.set_yticks(np.arange(N-1)+0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 4])
ax.text(0.35, 0.5, '=', horizontalalignment='center',
        verticalalignment='center', size=40, fontweight='bold')
ax.axis('off')
ax = plt.subplot(gs[0, 5])
ax.imshow(u[:, np.newaxis].T, cmap='rainbow')
ax.set_title(r'$u^T$', size=20, fontweight='bold')
ax.set_xticks(np.arange(M-1)+0.5)
ax.set_yticks([])
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 6:8])
ax.imshow(Op.A.T, cmap='rainbow')
ax.set_title(r'$(Op^T*$', size=20, fontweight='bold')
ax.set_xticks(np.arange(N-1)+0.5)
ax.set_yticks(np.arange(M-1)+0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis('tight')
ax = plt.subplot(gs[0, 8])
ax.imshow(v[:, np.newaxis], cmap='rainbow')
ax.set_title(r'$v)$', size=20, fontweight='bold')
ax.set_xticks([])
ax.set_yticks(np.arange(N-1)+0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])

Out:

[Text(0, 0.5, ''), Text(0, 1.5, ''), Text(0, 2.5, ''), Text(0, 3.5, '')]

From now on, we can simply use the pylops.utils.dottest implementation of the dot-test and pass the operator we would like to validate, its size in the model and data spaces and optionally the tolerance we will be accepting for the dot-test to be considered succesfull. Finally we need to specify if our data or/and model vectors contain complex numbers using the complexflag parameter. While the dot-test will return True when succesfull and False otherwise, we can also ask to print its outcome putting the verb parameters to True.

N = 10
d = np.arange(N)
Dop = pylops.Diagonal(d)

dottest(Dop, N, N, tol=1e-6, complexflag=0, verb=True)

Out:

Dot test passed, v^T(Opu)=-6.919654 - u^T(Op^Tv)=-6.919654

True

We move now to a more complicated operator, the pylops.signalprocessing.FFT operator. We use once again the pylops.utils.dottest to verify its implementation and since we are dealing with a transform that can be applied to both real and complex array, we try different combinations using the complexflag input.

dt = 0.005
nt = 100
nfft = 2**10

FFTop = pylops.signalprocessing.FFT(dims=(nt,), nfft=nfft,
                                    sampling=dt, dtype=np.complex128)
dottest(FFTop, nfft, nt, complexflag=2, verb=True)
dottest(FFTop, nfft, nt, complexflag=3, verb=True)

Out:

Dot test passed, v^T(Opu)=-17.139104+3.828834i - u^T(Op^Tv)=-17.139104+3.828834i
Dot test passed, v^T(Opu)=11.420800-3.688151i - u^T(Op^Tv)=11.420800-3.688151i

True

03. Solvers

03. Solvers

This tutorial will guide you through the pylops.optimization module and show how to use various solvers that are included in the PyLops library.

The main idea here is to provide the user of PyLops with very high-level functionalities to quickly and easily set up and solve complex systems of linear equations as well as include regularization and/or preconditioning terms (all of those constructed by means of PyLops linear operators).

To make this tutorial more interesting, we will present a real life problem and show how the choice of the solver and regularization/preconditioning terms is vital in many circumstances to successfully retrieve an estimate of the model. The problem that we are going to consider is generally referred to as the data reconstruction problem and aims at reconstructing a regularly sampled signal of size \(M\) from \(N\) randomly selected samples:

\[\mathbf{y} = \mathbf{R} \mathbf{x}\]

where the restriction operator \(\mathbf{R}\) that selects the \(M\) elements from \(\mathbf{x}\) at random locations is implemented using pylops.Restriction, and

\[\mathbf{y}= [y_1, y_2,...,y_N]^T, \qquad \mathbf{x}= [x_1, x_2,...,x_M]^T, \qquad\]

with \(M>>N\).

# pylint: disable=C0103
import numpy as np
import matplotlib.pyplot as plt

import pylops

plt.close('all')
np.random.seed(10)

Let’s first create the data in the frequency domain. The data is composed by the superposition of 3 sinusoids with different frequencies.

# Signal creation in frequency domain
ifreqs = [41, 25, 66]
amps = [1., 1., 1.]
N = 200
nfft = 2**11
dt = 0.004
t = np.arange(N)*dt
f = np.fft.rfftfreq(nfft, dt)

FFTop = 10*pylops.signalprocessing.FFT(N, nfft=nfft, real=True)

X = np.zeros(nfft//2+1, dtype='complex128')
X[ifreqs] = amps
x = FFTop.H*X

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
axs[0].plot(f, np.abs(X), 'k', LineWidth=2)
axs[0].set_xlim(0, 30)
axs[0].set_title('Data(frequency domain)')
axs[1].plot(t, x, 'k', LineWidth=2)
axs[1].set_title('Data(time domain)')
axs[1].axis('tight')

Out:

(-0.0398, 0.8358000000000001, -0.7554257721191237, 1.4249324045744862)

We now define the locations at which the signal will be sampled.

# subsampling locations
perc_subsampling = 0.2
Nsub = int(np.round(N*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(N))[:Nsub])

# Create restriction operator
Rop = pylops.Restriction(N, iava, dtype='float64')

y = Rop*x
ymask = Rop.mask(x)

# Visualize data
fig = plt.figure(figsize=(12, 4))
plt.plot(t, x, 'k', lw=3)
plt.plot(t, x, '.k', ms=20, label='all samples')
plt.plot(t, ymask, '.g', ms=15, label='available samples')
plt.legend()
plt.title('Data restriction')

Out:

Text(0.5, 1.0, 'Data restriction')

To start let’s consider the simplest ‘solver’, i.e., least-square inversion without regularization. We aim here to minimize the following cost function:

\[J = ||\mathbf{y} - \mathbf{R} \mathbf{x}||_2^2\]

Depending on the choice of the operator \(\mathbf{R}\), such problem can be solved using explicit matrix solvers as well as iterative solvers. In this case we will be using the latter approach (more specifically the scipy implementation of the LSQR solver - i.e., scipy.sparse.linalg.lsqr) as we do not want to explicitly create and invert a matrix. In most cases this will be the only viable approach as most of the large-scale optimization problems that we are interested to solve using PyLops do not lend naturally to the creation and inversion of explicit matrices.

This first solver can be very easily implemented using the / for PyLops operators, which will automatically call the scipy.sparse.linalg.lsqr with some default parameters.

xinv = Rop / y

We can also use pylops.optimization.leastsquares.RegularizedInversion (without regularization term for now) and customize our solvers using kwargs.

xinv = \
    pylops.optimization.leastsquares.RegularizedInversion(Rop, [], y,
                                                          **dict(damp=0,
                                                                 iter_lim=10,
                                                                 show=1))

Out:

LSQR            Least-squares solution of  Ax = b
The matrix A has 40 rows and 200 columns
damp = 0.00000000000000e+00   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   3.759e+00  3.759e+00    1.0e+00  2.7e-01
     1  0.00000e+00   0.000e+00  0.000e+00    0.0e+00  0.0e+00   0.0e+00  0.0e+00

LSQR finished
Ax - b is small enough, given atol, btol

istop =       1   r1norm = 0.0e+00   anorm = 0.0e+00   arnorm = 0.0e+00
itn   =       1   r2norm = 0.0e+00   acond = 0.0e+00   xnorm  = 3.8e+00

Finally we can select a different starting guess from the null vector

xinv_fromx0 = \
    pylops.optimization.leastsquares.RegularizedInversion(Rop, [], y,
                                                          x0=np.ones(N),
                                                          **dict(damp=0,
                                                                 iter_lim=10,
                                                                 show=0))

The cost function above can be also expanded in terms of its normal equations

\[\mathbf{x}_{ne}= (\mathbf{R}^T \mathbf{R})^{-1} \mathbf{R}^T \mathbf{y}\]

The method pylops.optimization.leastsquares.NormalEquationsInversion implements such system of equations explicitly and solves them using an iterative scheme suitable for square matrices (i.e., \(M=N\)).

While this approach may seem not very useful, we will soon see how regularization terms could be easily added to the normal equations using this method.

xne = pylops.optimization.leastsquares.NormalEquationsInversion(Rop, [], y)

Let’s now visualize the different inversion results

fig = plt.figure(figsize=(12, 4))
plt.plot(t, x, 'k', lw=2, label='original')
plt.plot(t, xinv, 'b', ms=10, label='inversion')
plt.plot(t, xinv_fromx0, '--r', ms=10, label='inversion from x0')
plt.plot(t, xne, '--g', ms=10, label='normal equations')
plt.legend()
plt.title('Data reconstruction without regularization')

Out:

Text(0.5, 1.0, 'Data reconstruction without regularization')

Regularization

You may have noticed that none of the inversion has been successfull in recovering the original signal. This is a clear indication that the problem we are trying to solve is highly ill-posed and requires some prior knowledge from the user.

We will now see how to add prior information to the inverse process in the form of regularization (or preconditioning). This can be done in two different ways

• regularization via pylops.optimization.leastsquares.NormalEquationsInversion or pylops.optimization.leastsquares.RegularizedInversion)
• preconditioning via pylops.optimization.leastsquares.PreconditionedInversion

Let’s start by regularizing the normal equations using a second derivative operator

\[\mathbf{x} = (\mathbf{R^TR}+\epsilon_\nabla^2\nabla^T\nabla)^{-1} \mathbf{R^Ty}\]

# Create regularization operator
D2op = pylops.SecondDerivative(N, dims=None, dtype='float64')

# Regularized inversion
epsR = np.sqrt(0.1)
epsI = np.sqrt(1e-4)

xne = \
    pylops.optimization.leastsquares.NormalEquationsInversion(Rop, [D2op], y,
                                                              epsI=epsI,
                                                              epsRs=[epsR],
                                                              returninfo=False,
                                                              **dict(maxiter=50))

Note that in case we have access to a fast implementation for the chain of forward and adjoint for the regularization operator (i.e., \(\nabla^T\nabla\)), we can modify our call to pylops.optimization.leastsquares.NormalEquationsInversion as follows:

ND2op = pylops.MatrixMult((D2op.H * D2op).tosparse()) # mimic fast D^T D

xne1 = \
    pylops.optimization.leastsquares.NormalEquationsInversion(Rop, [], y,
                                                              NRegs=[ND2op],
                                                              epsI=epsI,
                                                              epsNRs=[epsR],
                                                              returninfo=False,
                                                              **dict(maxiter=50))

We can do the same while using pylops.optimization.leastsquares.RegularizedInversion which solves the following augmented problem

\[\begin{split}\begin{bmatrix} \mathbf{R} \\ \epsilon_\nabla \nabla \end{bmatrix} \mathbf{x} = \begin{bmatrix} \mathbf{y} \\ 0 \end{bmatrix}\end{split}\]

xreg = \
    pylops.optimization.leastsquares.RegularizedInversion(Rop, [D2op], y,
                                                          epsRs=[np.sqrt(0.1)],
                                                          returninfo=False,
                                                          **dict(damp=np.sqrt(1e-4),
                                                                 iter_lim=50,
                                                                 show=0))

We can also write a preconditioned problem, whose cost function is

\[J= ||\mathbf{y} - \mathbf{R} \mathbf{P} \mathbf{p}||_2^2\]

where \(\mathbf{P}\) is the precondioned operator, \(\mathbf{p}\) is the projected model in the preconditioned space, and \(\mathbf{x}=\mathbf{P}\mathbf{p}\) is the model in the original model space we want to solve for. Note that a preconditioned problem converges much faster to its solution than its corresponding regularized problem. This can be done using the routine pylops.optimization.leastsquares.PreconditionedInversion.

# Create regularization operator
Sop = pylops.Smoothing1D(nsmooth=11, dims=[N], dtype='float64')

# Invert for interpolated signal
xprec = \
    pylops.optimization.leastsquares.PreconditionedInversion(Rop, Sop, y,
                                                             returninfo=False,
                                                             **dict(damp=np.sqrt(1e-9),
                                                                    iter_lim=20,
                                                                    show=0))

Let’s finally visualize these solutions

# sphinx_gallery_thumbnail_number=4
fig = plt.figure(figsize=(12, 4))
plt.plot(t[iava], y, '.k', ms=20, label='available samples')
plt.plot(t, x, 'k', lw=3, label='original')
plt.plot(t, xne, 'b', lw=3, label='normal equations')
plt.plot(t, xne1, '--c', lw=3, label='normal equations (with direct D^T D)')
plt.plot(t, xreg, '-.r', lw=3, label='regularized')
plt.plot(t, xprec, '--g', lw=3, label='preconditioned equations')
plt.legend()
plt.title('Data reconstruction with regularization')

subax = fig.add_axes([0.7, 0.2, 0.15, 0.6])
subax.plot(t[iava], y, '.k', ms=20)
subax.plot(t, x, 'k', lw=3)
subax.plot(t, xne, 'b', lw=3)
subax.plot(t, xne1, '--c', lw=3)
subax.plot(t, xreg, '-.r', lw=3)
subax.plot(t, xprec, '--g', lw=3)
subax.set_xlim(0.05, 0.3)

Out:

(0.05, 0.3)

Much better estimates! We have seen here how regularization and/or preconditioning can be vital to succesfully solve some ill-posed inverse problems.

We have however so far only considered solvers that can include additional norm-2 regularization terms. A very active area of research is that of sparsity-promoting solvers (also sometimes referred to as compressive sensing): the regularization term added to the cost function to minimize has norm-p (\(p \le 1\)) and the problem is generally recasted by considering the model to be sparse in some domain. We can follow this philosophy as our signal to invert was actually created as superposition of 3 sinusoids (i.e., three spikes in the Fourier domain). Our new cost function is:

\[J_1 = ||\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{p}||_2^2 + \epsilon ||\mathbf{p}||_1\]

where \(\mathbf{F}\) is the FFT operator. We will thus use the pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA solvers to estimate our input signal.

pista, niteri, costi = \
    pylops.optimization.sparsity.ISTA(Rop*FFTop.H, y, niter=1000,
                                      eps=0.1, tol=1e-7, returninfo=True)
xista = FFTop.H*pista

pfista, niterf, costf = \
    pylops.optimization.sparsity.FISTA(Rop*FFTop.H, y, niter=1000,
                                       eps=0.1, tol=1e-7, returninfo=True)
xfista = FFTop.H*pfista

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle('Data reconstruction with sparsity', fontsize=14,
             fontweight='bold', y=0.9)
axs[0].plot(f, np.abs(X), 'k', lw=3)
axs[0].plot(f, np.abs(pista), '--r', lw=3)
axs[0].plot(f, np.abs(pfista), '--g', lw=3)
axs[0].set_xlim(0, 30)
axs[0].set_title('Frequency domain')
axs[1].plot(t[iava], y, '.k', ms=20, label='available samples')
axs[1].plot(t, x, 'k', lw=3, label='original')
axs[1].plot(t, xista, '--r', lw=3, label='ISTA')
axs[1].plot(t, xfista, '--g', lw=3, label='FISTA')
axs[1].set_title('Time domain')
axs[1].axis('tight')
axs[1].legend()
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, ax = plt.subplots(1, 1, figsize=(12, 3))
ax.semilogy(costi, 'r', lw=2, label='ISTA')
ax.semilogy(costf, 'g', lw=2, label='FISTA')
ax.set_title('Cost functions', size=15, fontweight='bold')
ax.set_xlabel('Iteration')
ax.legend()
ax.grid(True)
plt.tight_layout()

• 
• 

As you can see, changing parametrization of the model and imposing sparsity in the Fourier domain has given an extra improvement to our ability of recovering the underlying densely sampled input signal. Moreover, FISTA converges much faster than ISTA as expected and should be preferred when using sparse solvers.

Finally we consider a slightly different cost function (note that in this case we try to solve a constrained problem):

\[J_1 = ||\mathbf{p}||_1 \quad subj.to \quad ||\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{p}||\]

A very popular solver to solve such kind of cost function is called spgl1 and can be accessed via pylops.optimization.sparsity.SPGL1.

xspgl1, pspgl1, info = \
    pylops.optimization.sparsity.SPGL1(Rop, y, FFTop, tau=3, iter_lim=200)

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle('Data reconstruction with SPGL1', fontsize=14,
             fontweight='bold', y=0.9)
axs[0].plot(f, np.abs(X), 'k', lw=3)
axs[0].plot(f, np.abs(pspgl1), '--m', lw=3)
axs[0].set_xlim(0, 30)
axs[0].set_title('Frequency domain')
axs[1].plot(t[iava], y, '.k', ms=20, label='available samples')
axs[1].plot(t, x, 'k', lw=3, label='original')
axs[1].plot(t, xspgl1, '--m', lw=3, label='SPGL1')
axs[1].set_title('Time domain')
axs[1].axis('tight')
axs[1].legend()
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, ax = plt.subplots(1, 1, figsize=(12, 3))
ax.semilogy(info['rnorm2'], 'k', lw=2, label='ISTA')
ax.set_title('Cost functions', size=15, fontweight='bold')
ax.set_xlabel('Iteration')
ax.legend()
ax.grid(True)
plt.tight_layout()

• 
• 

04. Bayesian Inversion

04. Bayesian Inversion

This tutorial focuses on Bayesian inversion, a special type of inverse problem that aims at incorporating prior information in terms of model and data probabilities in the inversion process.

In this case we will be dealing with the same problem that we discussed in 03. Solvers, but instead of defining ad-hoc regularization or preconditioning terms we parametrize and model our input signal in the frequency domain in a probabilistic fashion: the central frequency, amplitude and phase of the three sinusoids have gaussian distributions as follows:

\[X(f) = \sum_{i=1}^3 a_i e^{j \phi_i} \delta(f - f_i)\]

where \(f_i \sim N(f_{0,i}, \sigma_{f,i})\), \(a_i \sim N(a_{0,i}, \sigma_{a,i})\), and \(\phi_i \sim N(\phi_{0,i}, \sigma_{\phi,i})\).

Based on the above definition, we construct some prior models in the frequency domain, convert each of them to the time domain and use such an ensemble to estimate the prior mean \(\mu_\mathbf{x}\) and model covariance \(\mathbf{C_x}\).

We then create our data by sampling the true signal at certain locations and solve the resconstruction problem within a Bayesian framework. Since we are assuming gaussianity in our priors, the equation to obtain the posterion mean can be derived analytically:

\[\mathbf{x} = \mathbf{x_0} + \mathbf{C}_x \mathbf{R}^T (\mathbf{R} \mathbf{C}_x \mathbf{R}^T + \mathbf{C}_y)^{-1} (\mathbf{y} - \mathbf{R} \mathbf{x_0})\]

# sphinx_gallery_thumbnail_number = 2
import numpy as np
import matplotlib.pyplot as plt
import pylops

from scipy.sparse.linalg import lsqr


plt.close('all')
np.random.seed(10)

Let’s start by creating our true model and prior realizations

def prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft):
    """Create realization from prior mean and std for amplitude, frequency and
    phase
    """
    f = np.fft.rfftfreq(nfft, dt)
    df = f[1] - f[0]
    ifreqs = [int(np.random.normal(f, sigma)/df)
              for f, sigma in zip(f0, sigmaf)]
    amps = [np.random.normal(a, sigma) for a, sigma in zip(a0, sigmaa)]
    phis = [np.random.normal(phi, sigma) for phi, sigma in zip(phi0, sigmaphi)]

    # input signal in frequency domain
    X = np.zeros(nfft//2+1, dtype='complex128')
    X[ifreqs] = np.array(amps).squeeze() * \
                np.exp(1j * np.deg2rad(np.array(phis))).squeeze()

    # input signal in time domain
    FFTop = pylops.signalprocessing.FFT(nt, nfft=nfft, real=True)
    x = FFTop.H*X
    return x

# Priors
nreals = 100
f0 = [5, 3, 8]
sigmaf = [0.5, 1., 0.6]
a0 = [1., 1., 1.]
sigmaa = [0.1, 0.5, 0.6]
phi0 = [-90., 0., 0.]
sigmaphi = [0.1, 0.2, 0.4]
sigmad = 1e-2

# Prior models
nt = 200
nfft = 2**11
dt = 0.004
t = np.arange(nt)*dt
xs = \
    np.array([prior_realization(f0, a0, phi0, sigmaf,
                                sigmaa, sigmaphi, dt, nt, nfft)
              for _ in range(nreals)])

# True model (taken as one possible realization)
x = prior_realization(f0, a0, phi0, [0, 0, 0], [0, 0, 0],
                      [0, 0, 0], dt, nt, nfft)

We have now a set of prior models in time domain. We can easily use sample statistics to estimate the prior mean and covariance. For the covariance, we perform a second step where we average values around the main diagonal for each row and find a smooth, compact filter that we use to define a convolution linear operator that mimics the action of the covariance matrix on a vector

x0 = np.average(xs, axis=0)
Cm = ((xs - x0).T @ (xs - x0))/nreals

N = 30 # lenght of decorrelation
diags = np.array([Cm[i, i-N:i+N+1] for i in range(N, nt-N)])
diag_ave = np.average(diags, axis=0)
diag_ave *= np.hamming(2*N+1) # add a taper at the end to avoid edge effects

fig, ax = plt.subplots(1, 1, figsize=(12, 4))
ax.plot(t, xs.T, 'r', lw=1)
ax.plot(t, x0, 'g', lw=4)
ax.plot(t, x, 'k', lw=4)
ax.set_title('Prior realizations and mean')
ax.set_xlim(0, 0.8)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
im = ax1.imshow(Cm, interpolation='nearest', cmap='seismic',
                extent=(t[0], t[-1], t[-1], t[0]))
ax1.set_title(r"$\mathbf{C}_m^{prior}$")
ax1.axis('tight')
ax2.plot(np.arange(-N, N + 1) * dt, diags.T, '--r', lw=1)
ax2.plot(np.arange(-N, N + 1) * dt, diag_ave, 'k', lw=4)
ax2.set_title("Averaged covariance 'filter'")

• 
• 

Out:

Text(0.5, 1.0, "Averaged covariance 'filter'")

Let’s define now the sampling operator as well as create our covariance matrices in terms of linear operators. This may not be strictly necessary here but shows how even Bayesian-type of inversion can very easily scale to large model and data spaces.

# Sampling operator
perc_subsampling = 0.2
ntsub = int(np.round(nt*perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(nt))[:ntsub])
iava[-1] = nt-1 # assume we have the last sample to avoid instability
Rop = pylops.Restriction(nt, iava, dtype='float64')

# Covariance operators
Cm_op = \
    pylops.signalprocessing.Convolve1D(nt, diag_ave, offset=N)
Cd_op = sigmad**2 * pylops.Identity(ntsub)

We model now our data and add noise that respects our prior definition

n = np.random.normal(0, sigmad, nt)
y = Rop * x
yn = Rop * (x + n)
ymask = Rop.mask(x)
ynmask = Rop.mask(x + n)

First we apply the Bayesian inversion equation

xbayes = x0 + Cm_op * Rop.H * (lsqr(Rop * Cm_op * Rop.H + Cd_op,
                                    yn - Rop*x0, iter_lim=400)[0])

# Visualize
fig, ax = plt.subplots(1, 1, figsize=(12, 5))
ax.plot(t, x, 'k', lw=6, label='true')
ax.plot(t, ymask, '.k', ms=25, label='available samples')
ax.plot(t, ynmask, '.r', ms=25, label='available noisy samples')
ax.plot(t, xbayes, 'r', lw=3, label='bayesian inverse')
ax.legend()
ax.set_title('Signal')
ax.set_xlim(0, 0.8)

Out:

(0.0, 0.8)

So far we have been able to estimate our posterion mean. What about its uncertainties (i.e., posterion covariance)?

In real-life applications it is very difficult (if not impossible) to directly compute the posterior covariance matrix. It is much more useful to create a set of models that sample the posterion probability. We can do that by solving our problem several times using different prior realizations as starting guesses:

xpost = [x0 + Cm_op * Rop.H * (lsqr(Rop * Cm_op * Rop.H + Cd_op, yn - Rop * x0, iter_lim=400)[0])
         for x0 in xs[:30]]
xpost = np.array(xpost)

x0post = np.average(xpost, axis=0)
Cm_post = ((xpost - x0post).T @ (xpost - x0post)) / nreals

# Visualize
fig, ax = plt.subplots(1, 1, figsize=(12, 5))
ax.plot(t, x, 'k', lw=6, label='true')
ax.plot(t, xpost.T, '--r', lw=1)
ax.plot(t, x0post, 'r', lw=3, label='bayesian inverse')
ax.plot(t, ymask, '.k', ms=25, label='available samples')
ax.plot(t, ynmask, '.r', ms=25, label='available noisy samples')
ax.legend()
ax.set_title('Signal')
ax.set_xlim(0, 0.8)

fig, ax = plt.subplots(1, 1, figsize=(5, 4))
im = ax.imshow(Cm_post, interpolation='nearest', cmap='seismic',
               extent=(t[0], t[-1], t[-1], t[0]))
ax.set_title(r"$\mathbf{C}_m^{posterior}$")
ax.axis('tight')

• 
• 

Out:

(0.0, 0.796, 0.796, 0.0)

Note that here we have been able to compute a sample posterior covariance from its estimated samples. By displaying it we can see how both the overall variances and the correlation between different parameters have become narrower compared to their prior counterparts.

05. Image deblurring

05. Image deblurring

Deblurring is the process of removing blurring effects from images, caused for example by defocus aberration or motion blur.

In forward mode, such blurring effect is typically modelled as a 2-dimensional convolution between the so-called point spread function and a target sharp input image, where the sharp input image (which has to be recovered) is unknown and the point-spread function can be either known or unknown.

In this tutorial, an example of 2d blurring and deblurring will be shown using the pylops.signalprocessing.Convolve2D operator assuming knowledge of the point-spread function.

import numpy as np
import matplotlib.pyplot as plt
import pylops

Let’s start by importing a 2d image and defining the blurring operator

im = np.load('../testdata/python.npy')[::5, ::5, 0]

Nz, Nx = im.shape

# Blurring guassian operator
nh = [15, 25]
hz = np.exp(-0.1*np.linspace(-(nh[0]//2), nh[0]//2, nh[0])**2)
hx = np.exp(-0.03*np.linspace(-(nh[1]//2), nh[1]//2, nh[1])**2)
hz /= np.trapz(hz) # normalize the integral to 1
hx /= np.trapz(hx) # normalize the integral to 1
h = hz[:, np.newaxis] * hx[np.newaxis, :]

fig, ax = plt.subplots(1, 1, figsize=(5, 3))
him = ax.imshow(h)
ax.set_title('Blurring operator')
fig.colorbar(him, ax=ax)
ax.axis('tight')

Cop = pylops.signalprocessing.Convolve2D(Nz * Nx, h=h,
                                         offset=(nh[0] // 2,
                                                 nh[1] // 2),
                                         dims=(Nz, Nx), dtype='float32')

We first apply the blurring operator to the sharp image. We then try to recover the sharp input image by inverting the convolution operator from the blurred image. Note that when we perform inversion without any regularization, the deblurred image will show some ringing due to the instabilities of the inverse process. Using a L1 solver with a DWT preconditioner or TV regularization allows to recover sharper contrasts.

imblur = Cop * im.flatten()

imdeblur = \
    pylops.optimization.leastsquares.NormalEquationsInversion(Cop, None,
                                                              imblur,
                                                              maxiter=50)

Wop = pylops.signalprocessing.DWT2D((Nz, Nx), wavelet='haar', level=3)
Dop = [pylops.FirstDerivative(Nz * Nx, dims=(Nz, Nx), dir=0, edge=False),
       pylops.FirstDerivative(Nz * Nx, dims=(Nz, Nx), dir=1, edge=False)]
DWop = Dop + [Wop, ]

imdeblurfista = \
    pylops.optimization.sparsity.FISTA(Cop * Wop.H, imblur, eps=1e-1,
                                       niter=100)[0]
imdeblurfista = Wop.H * imdeblurfista

imdeblurtv = \
    pylops.optimization.sparsity.SplitBregman(Cop, Dop, imblur.flatten(),
                                              niter_outer=10, niter_inner=5,
                                              mu=1.5, epsRL1s=[2e0, 2e0],
                                              tol=1e-4, tau=1., show=False,
                                              ** dict(iter_lim=5, damp=1e-4))[0]

imdeblurtv1 = \
    pylops.optimization.sparsity.SplitBregman(Cop, DWop,
                                              imblur.flatten(),
                                              niter_outer=10, niter_inner=5,
                                              mu=1.5, epsRL1s=[1e0, 1e0, 1e0],
                                              tol=1e-4, tau=1., show=False,
                                              ** dict(iter_lim=5, damp=1e-4))[0]

# Reshape images
imblur = imblur.reshape((Nz, Nx))
imdeblur = imdeblur.reshape((Nz, Nx))
imdeblurfista = imdeblurfista.reshape((Nz, Nx))
imdeblurtv = imdeblurtv.reshape((Nz, Nx))
imdeblurtv1 = imdeblurtv1.reshape((Nz, Nx))

Finally we visualize the original, blurred, and recovered images.

# sphinx_gallery_thumbnail_number = 2
fig = plt.figure(figsize=(12, 6))
fig.suptitle('Deblurring', fontsize=14, fontweight='bold', y=0.95)
ax1 = plt.subplot2grid((2, 5), (0, 0))
ax2 = plt.subplot2grid((2, 5), (0, 1))
ax3 = plt.subplot2grid((2, 5), (0, 2))
ax4 = plt.subplot2grid((2, 5), (1, 0))
ax5 = plt.subplot2grid((2, 5), (1, 1))
ax6 = plt.subplot2grid((2, 5), (1, 2))
ax7 = plt.subplot2grid((2, 5), (0, 3), colspan=2)
ax8 = plt.subplot2grid((2, 5), (1, 3), colspan=2)
ax1.imshow(im, cmap='viridis', vmin=0, vmax=250)
ax1.axis('tight')
ax1.set_title('Original')
ax2.imshow(imblur, cmap='viridis', vmin=0, vmax=250)
ax2.axis('tight')
ax2.set_title('Blurred')
ax3.imshow(imdeblur, cmap='viridis', vmin=0, vmax=250)
ax3.axis('tight')
ax3.set_title('Deblurred')
ax4.imshow(imdeblurfista, cmap='viridis', vmin=0, vmax=250)
ax4.axis('tight')
ax4.set_title('FISTA deblurred')
ax5.imshow(imdeblurtv, cmap='viridis', vmin=0, vmax=250)
ax5.axis('tight')
ax5.set_title('TV deblurred')
ax6.imshow(imdeblurtv1, cmap='viridis', vmin=0, vmax=250)
ax6.axis('tight')
ax6.set_title('TV+Haar deblurred')
ax7.plot(im[Nz//2], 'k')
ax7.plot(imblur[Nz//2], '--r')
ax7.plot(imdeblur[Nz//2], '--b')
ax7.plot(imdeblurfista[Nz//2], '--g')
ax7.plot(imdeblurtv[Nz//2], '--m')
ax7.plot(imdeblurtv1[Nz//2], '--y')
ax7.axis('tight')
ax7.set_title('Horizontal section')
ax8.plot(im[:, Nx//2], 'k', label='Original')
ax8.plot(imblur[:, Nx//2], '--r', label='Blurred')
ax8.plot(imdeblur[:, Nx//2], '--b', label='Deblurred')
ax8.plot(imdeblurfista[:, Nx//2], '--g', label='FISTA deblurred')
ax8.plot(imdeblurtv[:, Nx//2], '--m', label='TV deblurred')
ax8.plot(imdeblurtv1[:, Nx//2], '--y', label='TV+Haar deblurred')
ax8.axis('tight')
ax8.set_title('Vertical section')
ax8.legend(loc=5, fontsize='small')
plt.tight_layout()
plt.subplots_adjust(top=0.8)

06. 2D Interpolation

06. 2D Interpolation

In the mathematical field of numerical analysis, interpolation is the problem of constructing new data points within the range of a discrete set of known data points. In signal and image processing, the data may be recorded at irregular locations and it is often required to regularize the data into a regular grid.

In this tutorial, an example of 2d interpolation of an image is carried out using a combination of PyLops operators (pylops.Restriction and pylops.Laplacian) and the pylops.optimization module.

Mathematically speaking, if we want to interpolate a signal using the theory of inverse problems, we can define the following forward problem:

\[\mathbf{y} = \mathbf{R} \mathbf{x}\]

where the restriction operator \(\mathbf{R}\) selects \(M\) elements from the regularly sampled signal \(\mathbf{x}\) at random locations. The input and output signals are:

\[\mathbf{y}= [y_1, y_2,...,y_N]^T, \qquad \mathbf{x}= [x_1, x_2,...,x_M]^T, \qquad\]

with \(M>>N\).

import numpy as np
import matplotlib.pyplot as plt
import pylops

plt.close('all')
np.random.seed(0)

To start we import a 2d image and define our restriction operator to irregularly and randomly sample the image for 30% of the entire grid

im = np.load('../testdata/python.npy')[:, :, 0]

Nz, Nx = im.shape
N = Nz * Nx

# Subsample signal
perc_subsampling = 0.2

Nsub2d = int(np.round(N*perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(N))[:Nsub2d])

# Create operators and data
Rop = pylops.Restriction(N, iava, dtype='float64')
D2op = pylops.Laplacian((Nz, Nx), weights=(1, 1), dtype='float64')

x = im.flatten()
y = Rop * x
y1 = Rop.mask(x)

We will now use two different routines from our optimization toolbox to estimate our original image in the regular grid.

xcg_reg_lop = \
    pylops.optimization.leastsquares.NormalEquationsInversion(Rop, [D2op], y,
                                                              epsRs=[np.sqrt(0.1)],
                                                              returninfo=False,
                                                              **dict(maxiter=200))

# Invert for interpolated signal, lsqrt
xlsqr_reg_lop, istop, itn, r1norm, r2norm = \
    pylops.optimization.leastsquares.RegularizedInversion(Rop, [D2op], y,
                                                          epsRs=[np.sqrt(0.1)],
                                                          returninfo=True,
                                                          **dict(damp=0,
                                                                 iter_lim=200,
                                                                 show=0))

# Reshape estimated images
im_sampled = y1.reshape((Nz, Nx))
im_rec_lap_cg = xcg_reg_lop.reshape((Nz, Nx))
im_rec_lap_lsqr = xlsqr_reg_lop.reshape((Nz, Nx))

Finally we visualize the original image, the reconstructed images and their error

fig, axs = plt.subplots(1, 4, figsize=(12, 4))
fig.suptitle('Data reconstruction - normal eqs', fontsize=14,
             fontweight='bold', y=0.95)
axs[0].imshow(im, cmap='viridis', vmin=0, vmax=250)
axs[0].axis('tight')
axs[0].set_title('Original')
axs[1].imshow(im_sampled, cmap='viridis', vmin=0, vmax=250)
axs[1].axis('tight')
axs[1].set_title('Sampled')
axs[2].imshow(im_rec_lap_cg, cmap='viridis', vmin=0, vmax=250)
axs[2].axis('tight')
axs[2].set_title('2D Regularization')
axs[3].imshow(im - im_rec_lap_cg, cmap='gray', vmin=-80, vmax=80)
axs[3].axis('tight')
axs[3].set_title('2D Regularization Error')
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, axs = plt.subplots(1, 4, figsize=(12, 4))
fig.suptitle('Data reconstruction - regularized eqs', fontsize=14,
             fontweight='bold', y=0.95)
axs[0].imshow(im, cmap='viridis', vmin=0, vmax=250)
axs[0].axis('tight')
axs[0].set_title('Original')
axs[1].imshow(im_sampled, cmap='viridis', vmin=0, vmax=250)
axs[1].axis('tight')
axs[1].set_title('Sampled')
axs[2].imshow(im_rec_lap_lsqr, cmap='viridis', vmin=0, vmax=250)
axs[2].axis('tight')
axs[2].set_title('2D Regularization')
axs[3].imshow(im - im_rec_lap_lsqr, cmap='gray', vmin=-80, vmax=80)
axs[3].axis('tight')
axs[3].set_title('2D Regularization Error')
plt.tight_layout()
plt.subplots_adjust(top=0.8)

• 
• 

07. Post-stack inversion

07. Post-stack inversion

Estimating subsurface properties from band-limited seismic data represents an important task for geophysical subsurface characterization.

In this tutorial, the pylops.avo.poststack.PoststackLinearModelling operator is used for modelling of both 1d and 2d synthetic post-stack seismic data from a profile or 2d model of the subsurface acoustic impedence.

\[d(t, \theta=0) = \frac{1}{2} w(t) * \frac{d\ln(AI(t))}{\quad dt}\]

where \(AI(t)\) is the acoustic impedance profile and \(w(t)\) is the time domain seismic wavelet. In compact form:

\[\mathbf{d}= \mathbf{W} \mathbf{D} \mathbf{ai}\]

where \(\mathbf{W}\) is a convolution operator, \(\mathbf{D}\) is a first derivative operator, and \(\mathbf{ai}\) is the input model. Subsequently the acoustic impedance model is estimated via the pylops.avo.poststack.PoststackInversion module. A two-steps inversion strategy is finally presented to deal with the case of noisy data.

# sphinx_gallery_thumbnail_number = 4
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import filtfilt

import pylops
from pylops.utils.wavelets import ricker

plt.close('all')
np.random.seed(10)

Let’s start with a 1d example. A synthetic profile of acoustic impedance is created and data is modelled using both the dense and linear operator version of pylops.avo.poststack.PoststackLinearModelling operator.

# model
nt0 = 301
dt0 = 0.004
t0 = np.arange(nt0)*dt0
vp = 1200 + np.arange(nt0) + \
     filtfilt(np.ones(5)/5., 1, np.random.normal(0, 80, nt0))
rho = 1000 + vp + \
      filtfilt(np.ones(5)/5., 1, np.random.normal(0, 30, nt0))
vp[131:] += 500
rho[131:] += 100
m = np.log(vp*rho)

# smooth model
nsmooth = 100
mback = filtfilt(np.ones(nsmooth)/float(nsmooth), 1, m)

# wavelet
ntwav = 41
wav, twav, wavc = ricker(t0[:ntwav//2+1], 20)

# dense operator
PPop_dense = \
    pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nt0, explicit=True)

# lop operator
PPop = pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nt0)

# data
d_dense = PPop_dense*m.flatten()
d = PPop*m

# add noise
dn_dense = d_dense + np.random.normal(0, 2e-2, d_dense.shape)

We can now estimate the acoustic profile from band-limited data using either the dense operator or linear operator.

# solve dense
minv_dense = \
    pylops.avo.poststack.PoststackInversion(d, wav / 2, m0=mback, explicit=True,
                                            simultaneous=False)[0]

# solve lop
minv = \
    pylops.avo.poststack.PoststackInversion(d_dense, wav / 2, m0=mback,
                                            explicit=False,
                                            simultaneous=False,
                                            **dict(iter_lim=2000))[0]

# solve noisy
mn = \
    pylops.avo.poststack.PoststackInversion(dn_dense, wav / 2, m0=mback,
                                            explicit=True,
                                            epsR=1e0, **dict(damp=1e-1))[0]

fig, axs = plt.subplots(1, 2, figsize=(6, 7), sharey=True)
axs[0].plot(d_dense, t0, 'k', lw=4, label='Dense')
axs[0].plot(d, t0, '--r', lw=2, label='Lop')
axs[0].plot(dn_dense, t0, '-.g', lw=2, label='Noisy')
axs[0].set_title('Data')
axs[0].invert_yaxis()
axs[0].axis('tight')
axs[0].legend(loc=1)
axs[1].plot(m, t0, 'k', lw=4, label='True')
axs[1].plot(mback, t0, '--b', lw=4, label='Back')
axs[1].plot(minv_dense, t0, '--m', lw=2, label='Inv Dense')
axs[1].plot(minv, t0, '--r', lw=2, label='Inv Lop')
axs[1].plot(mn, t0, '--g', lw=2, label='Inv Noisy')
axs[1].set_title('Model')
axs[1].axis('tight')
axs[1].legend(loc=1)

Out:

<matplotlib.legend.Legend object at 0x7fc20a9bd518>

We see how inverting a dense matrix is in this case faster than solving for the linear operator (a good estimate of the model is in fact obtained only after 2000 iterations of lsqr). Nevertheless, having a linear operator is useful when we deal with larger dimensions (2d or 3d) and we want to couple our modelling operator with different types of spatial regularizations or preconditioning.

Before we move onto a 2d example, let’s consider the case of non-stationary wavelet and see how we can easily use the same routines in this case

# wavelet
ntwav = 41
f0s = np.flip(np.arange(nt0) * 0.05 + 3)
wavs = np.array([ricker(t0[:ntwav], f0)[0] for f0 in f0s])
wavc = np.argmax(wavs[0])

plt.figure(figsize=(5, 4))
plt.imshow(wavs.T, cmap='gray', extent=(t0[0], t0[-1], t0[ntwav], -t0[ntwav]))
plt.xlabel('t')
plt.title('Wavelets')
plt.axis('tight')

# operator
PPop = \
    pylops.avo.poststack.PoststackLinearModelling(wavs / 2, nt0=nt0, explicit=True)

# data
d = PPop * m

# solve
minv = \
    pylops.avo.poststack.PoststackInversion(d, wavs / 2, m0=mback, explicit=True,
                                            **dict(cond=1e-10))[0]

fig, axs = plt.subplots(1, 2, figsize=(6, 7), sharey=True)
axs[0].plot(d, t0, 'k', lw=4)
axs[0].set_title('Data')
axs[0].invert_yaxis()
axs[0].axis('tight')
axs[1].plot(m, t0, 'k', lw=4, label='True')
axs[1].plot(mback, t0, '--b', lw=4, label='Back')
axs[1].plot(minv, t0, '--r', lw=2, label='Inv')
axs[1].set_title('Model')
axs[1].axis('tight')
axs[1].legend(loc=1)

• 
• 

Out:

<matplotlib.legend.Legend object at 0x7fc20ac5eb38>

We move now to a 2d example. First of all the model is loaded and data generated.

# model
inputfile = '../testdata/avo/poststack_model.npz'

model = np.load(inputfile)
m = np.log(model['model'][:, ::3])
x, z = model['x'][::3]/1000., model['z']/1000.
nx, nz = len(x), len(z)

# smooth model
nsmoothz, nsmoothx = 60, 50
mback = filtfilt(np.ones(nsmoothz)/float(nsmoothz), 1, m, axis=0)
mback = filtfilt(np.ones(nsmoothx)/float(nsmoothx), 1, mback, axis=1)

# dense operator
PPop_dense = \
    pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nz,
                                                  spatdims=nx, explicit=True)

# lop operator
PPop = pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nz,
                                                     spatdims=nx)

# data
d = (PPop_dense * m.flatten()).reshape(nz, nx)
n = np.random.normal(0, 1e-1, d.shape)
dn = d + n

Finally we perform 4 different inversions:

• trace-by-trace inversion with explicit solver and dense operator with noise-free data

• trace-by-trace inversion with explicit solver and dense operator with noisy data

• multi-trace regularized inversion with iterative solver and linear operator using the result of trace-by-trace inversion as starting guess

  \[J = ||\Delta \mathbf{d} - \mathbf{W} \Delta \mathbf{ai}||_2 + \epsilon_\nabla ^2 ||\nabla \mathbf{ai}||_2\]

  where \(\Delta \mathbf{d}=\mathbf{d}-\mathbf{W}\mathbf{AI_0}\) is the residual data

• multi-trace blocky inversion with iterative solver and linear operator

# dense inversion with noise-free data
minv_dense = \
    pylops.avo.poststack.PoststackInversion(d, wav / 2, m0=mback,
                                            explicit=True,
                                            simultaneous=False)[0]

# dense inversion with noisy data
minv_dense_noisy = \
    pylops.avo.poststack.PoststackInversion(dn, wav / 2, m0=mback,
                                            explicit=True, epsI=4e-2,
                                            simultaneous=False)[0]

# spatially regularized lop inversion with noisy data
minv_lop_reg = \
    pylops.avo.poststack.PoststackInversion(dn, wav / 2, m0=minv_dense_noisy,
                                            explicit=False, epsR=5e1,
                                            **dict(damp=np.sqrt(1e-4),
                                                   iter_lim=80))[0]

# blockiness promoting inversion with noisy data
minv_lop_blocky = \
    pylops.avo.poststack.PoststackInversion(dn, wav / 2, m0=mback,
                                            explicit=False,
                                            epsR=[0.4], epsRL1=[0.1],
                                            **dict(mu=0.1,
                                                   niter_outer=5,
                                                   niter_inner=10,
                                                   iter_lim=5, damp=1e-3))[0]

fig, axs = plt.subplots(2, 4, figsize=(15, 9))
axs[0][0].imshow(d, cmap='gray',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=-0.4, vmax=0.4)
axs[0][0].set_title('Data')
axs[0][0].axis('tight')
axs[0][1].imshow(dn, cmap='gray',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=-0.4, vmax=0.4)
axs[0][1].set_title('Noisy Data')
axs[0][1].axis('tight')
axs[0][2].imshow(m, cmap='gist_rainbow',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=m.min(), vmax=m.max())
axs[0][2].set_title('Model')
axs[0][2].axis('tight')
axs[0][3].imshow(mback, cmap='gist_rainbow',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=m.min(), vmax=m.max())
axs[0][3].set_title('Smooth Model')
axs[0][3].axis('tight')
axs[1][0].imshow(minv_dense, cmap='gist_rainbow',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=m.min(), vmax=m.max())
axs[1][0].set_title('Noise-free Inversion')
axs[1][0].axis('tight')
axs[1][1].imshow(minv_dense_noisy, cmap='gist_rainbow',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=m.min(), vmax=m.max())
axs[1][1].set_title('Trace-by-trace Noisy Inversion')
axs[1][1].axis('tight')
axs[1][2].imshow(minv_lop_reg, cmap='gist_rainbow',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=m.min(), vmax=m.max())
axs[1][2].set_title('Regularized Noisy Inversion - lop ')
axs[1][2].axis('tight')
axs[1][3].imshow(minv_lop_blocky, cmap='gist_rainbow',
                 extent=(x[0], x[-1], z[-1], z[0]),
                 vmin=m.min(), vmax=m.max())
axs[1][3].set_title('Blocky Noisy Inversion - lop ')
axs[1][3].axis('tight')

fig, ax = plt.subplots(1, 1, figsize=(3, 7))
ax.plot(m[:, nx//2], z, 'k', lw=4, label='True')
ax.plot(mback[:, nx//2], z, '--r', lw=4, label='Back')
ax.plot(minv_dense[:, nx//2], z, '--b', lw=2, label='Inv Dense')
ax.plot(minv_dense_noisy[:, nx//2], z, '--m', lw=2, label='Inv Dense noisy')
ax.plot(minv_lop_reg[:, nx//2], z, '--g', lw=2, label='Inv Lop regularized')
ax.plot(minv_lop_blocky[:, nx//2], z, '--y', lw=2, label='Inv Lop blocky')
ax.set_title('Model')
ax.invert_yaxis()
ax.axis('tight')
ax.legend()
plt.tight_layout()

• 
• 

That’s almost it. If you wonder how this can be applied to real data, head over to the following notebook where the open-source segyio library is used alongside pylops to create an end-to-end open-source seismic inversion workflow with SEG-Y input data.

08. Pre-stack (AVO) inversion

08. Pre-stack (AVO) inversion

Pre-stack inversion represents one step beyond post-stack inversion in that not only the profile of acoustic impedance can be inferred from seismic data, rather a set of elastic parameters is estimated from pre-stack data (i.e., angle gathers) using the information contained in the so-called AVO (amplitude versus offset) response. Such elastic parameters represent vital information for more sophisticated geophysical subsurface characterization than it would be possible to achieve working with post-stack seismic data.

In this tutorial, the pylops.avo.prestack.PrestackLinearModelling operator is used for modelling of both 1d and 2d synthetic pre-stack seismic data using 1d profiles or 2d models of different subsurface elastic parameters (P-wave velocity, S-wave velocity, and density) as input.

\[d(t, \theta) = w(t) * \sum_{i=1}^N G_i(t, \theta) \frac{d\ln(m_i(t))}{dt}\]

where \(\mathbf{m}(t)=[V_P(t), V_S(t), \rho(t)]\) is a vector containing three elastic parameters at time \(t\), \(G_i(t, \theta)\) are the coefficients of the AVO parametrization used to model pre-stack data and \(w(t)\) is the time domain seismic wavelet. In compact form:

\[\mathbf{d}= \mathbf{W} \mathbf{G} \mathbf{D} \mathbf{m}\]

where \(\mathbf{W}\) is a convolution operator, \(\mathbf{G}\) is the AVO modelling operator, \(\mathbf{D}\) is a block-diagonal derivative operator, and \(\mathbf{m}\) is the input model. Subsequently the elastic parameters are estimated via the pylops.avo.prestack.PrestackInversion module. Once again, a two-steps inversion strategy can also be used to deal with the case of noisy data.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import filtfilt

import pylops
from pylops.utils.wavelets import ricker

plt.close('all')
np.random.seed(0)

Let’s start with a 1d example. A synthetic profile of acoustic impedance is created and data is modelled using both the dense and linear operator version of pylops.avo.prestack.PrestackLinearModelling operator

# sphinx_gallery_thumbnail_number = 5

# model
nt0 = 301
dt0 = 0.004

t0 = np.arange(nt0)*dt0
vp = 1200 + np.arange(nt0) + \
     filtfilt(np.ones(5)/5., 1, np.random.normal(0, 80, nt0))
vs = 600  + vp/2 + \
     filtfilt(np.ones(5)/5., 1, np.random.normal(0, 20, nt0))
rho = 1000 + vp + \
      filtfilt(np.ones(5)/5., 1, np.random.normal(0, 30, nt0))
vp[131:] += 500
vs[131:] += 200
rho[131:] += 100
vsvp = 0.5
m = np.stack((np.log(vp), np.log(vs), np.log(rho)), axis=1)

# background model
nsmooth = 50
mback = filtfilt(np.ones(nsmooth)/float(nsmooth), 1, m, axis=0)

# angles
ntheta = 21
thetamin, thetamax = 0, 40
theta = np.linspace(thetamin, thetamax, ntheta)

# wavelet
ntwav = 41
wav = ricker(t0[:ntwav//2+1], 15)[0]

# lop
PPop = \
    pylops.avo.prestack.PrestackLinearModelling(wav, theta,
                                                vsvp=vsvp, nt0=nt0,
                                                linearization='akirich')

# dense
PPop_dense = \
    pylops.avo.prestack.PrestackLinearModelling(wav, theta,
                                                vsvp=vsvp, nt0=nt0,
                                                linearization='akirich',
                                                explicit=True)

# data lop
dPP = PPop * m.flatten()
dPP = dPP.reshape(nt0, ntheta)

# data dense
dPP_dense = PPop_dense * m.T.flatten()
dPP_dense = dPP_dense.reshape(ntheta, nt0).T

# noisy data
dPPn_dense = dPP_dense + np.random.normal(0, 1e-2, dPP_dense.shape)

We can now invert our data and retrieve elastic profiles for both noise-free and noisy data using pylops.avo.prestack.PrestackInversion.

# dense
minv_dense, dPP_dense_res = \
    pylops.avo.prestack.PrestackInversion(dPP_dense, theta, wav, m0=mback,
                                          linearization='akirich',
                                          explicit=True, returnres=True,
                                          **dict(cond=1e-10))

# lop
minv, dPP_res = \
    pylops.avo.prestack.PrestackInversion(dPP, theta, wav, m0=mback,
                                          linearization='akirich',
                                          explicit=False, returnres=True,
                                          **dict(damp=1e-10, iter_lim=2000))

# dense noisy
minv_dense_noise, dPPn_dense_res = \
    pylops.avo.prestack.PrestackInversion(dPPn_dense, theta, wav, m0=mback,
                                          linearization='akirich',
                                          explicit=True,
                                          returnres=True, **dict(cond=1e-1))

# lop noisy (with vertical smoothing)
minv_noise, dPPn_res = \
    pylops.avo.prestack.PrestackInversion(dPPn_dense, theta, wav, m0=mback,
                                          linearization='akirich',
                                          explicit=False,
                                          epsR=5e-1, returnres=True,
                                          **dict(damp=1e-1, iter_lim=100))

The data, inverted models and residuals are now displayed.

# data and model
fig, (axd, axdn, axvp, axvs, axrho) = \
    plt.subplots(1, 5, figsize=(8, 5), sharey=True)
axd.imshow(dPP_dense, cmap='gray',
           extent=(theta[0], theta[-1], t0[-1], t0[0]),
           vmin=-np.abs(dPP_dense).max(), vmax=np.abs(dPP_dense).max())
axd.set_title('Data')
axd.axis('tight')
axdn.imshow(dPPn_dense, cmap='gray',
            extent=(theta[0], theta[-1], t0[-1], t0[0]),
            vmin=-np.abs(dPP_dense).max(), vmax=np.abs(dPP_dense).max())
axdn.set_title('Noisy Data')
axdn.axis('tight')
axvp.plot(vp, t0, 'k', lw=4, label='True')
axvp.plot(np.exp(mback[:, 0]), t0, '--r', lw=4, label='Back')
axvp.plot(np.exp(minv_dense[:, 0]), t0, '--b', lw=2, label='Inv Dense')
axvp.plot(np.exp(minv[:, 0]), t0, '--m', lw=2, label='Inv Lop')
axvp.plot(np.exp(minv_dense_noise[:, 0]), t0, '--c', lw=2, label='Noisy Dense')
axvp.plot(np.exp(minv_noise[:, 0]), t0, '--g', lw=2, label='Noisy Lop')
axvp.set_title(r'$V_P$')
axvs.plot(vs, t0, 'k', lw=4, label='True')
axvs.plot(np.exp(mback[:, 1]), t0, '--r', lw=4, label='Back')
axvs.plot(np.exp(minv_dense[:, 1]), t0, '--b', lw=2, label='Inv Dense')
axvs.plot(np.exp(minv[:, 1]), t0, '--m', lw=2, label='Inv Lop')
axvs.plot(np.exp(minv_dense_noise[:, 1]), t0, '--c', lw=2, label='Noisy Dense')
axvs.plot(np.exp(minv_noise[:, 1]), t0, '--g', lw=2, label='Noisy Lop')
axvs.set_title(r'$V_S$')
axrho.plot(rho, t0, 'k', lw=4, label='True')
axrho.plot(np.exp(mback[:, 2]), t0, '--r', lw=4, label='Back')
axrho.plot(np.exp(minv_dense[:, 2]), t0, '--b', lw=2, label='Inv Dense')
axrho.plot(np.exp(minv[:, 2]), t0, '--m', lw=2, label='Inv Lop')
axrho.plot(np.exp(minv_dense_noise[:, 2]),
           t0, '--c', lw=2, label='Noisy Dense')
axrho.plot(np.exp(minv_noise[:, 2]), t0, '--g', lw=2, label='Noisy Lop')
axrho.set_title(r'$\rho$')
axrho.legend(loc='center left', bbox_to_anchor=(1, 0.5))
axd.axis('tight')
plt.tight_layout()

# residuals
fig, axs = plt.subplots(1, 4, figsize=(8, 5), sharey=True)
fig.suptitle('Residuals', fontsize=14,
             fontweight='bold', y=0.95)
im = axs[0].imshow(dPP_dense_res, cmap='gray',
                   extent=(theta[0], theta[-1], t0[-1], t0[0]),
                   vmin=-0.1, vmax=0.1)
axs[0].set_title('Dense')
axs[0].set_xlabel(r'$\theta$')
axs[0].set_ylabel('t[s]')
axs[0].axis('tight')
axs[1].imshow(dPP_res, cmap='gray',
              extent=(theta[0], theta[-1], t0[-1], t0[0]),
              vmin=-0.1, vmax=0.1)
axs[1].set_title('Lop')
axs[1].set_xlabel(r'$\theta$')
axs[1].axis('tight')
axs[2].imshow(dPPn_dense_res, cmap='gray',
              extent=(theta[0], theta[-1], t0[-1], t0[0]),
              vmin=-0.1, vmax=0.1)
axs[2].set_title('Noisy Dense')
axs[2].set_xlabel(r'$\theta$')
axs[2].axis('tight')
axs[3].imshow(dPPn_res, cmap='gray',
              extent=(theta[0], theta[-1], t0[-1], t0[0]),
              vmin=-0.1, vmax=0.1)
axs[3].set_title('Noisy Lop')
axs[3].set_xlabel(r'$\theta$')
axs[3].axis('tight')
plt.tight_layout()
plt.subplots_adjust(top=0.85)

• 
• 

Finally before moving to the 2d example, we consider the case when both PP and PS data are available. A joint PP-PS inversion can be easily solved as follows.

PSop = \
    pylops.avo.prestack.PrestackLinearModelling(2 * wav, theta,
                                                vsvp=vsvp, nt0=nt0,
                                                linearization='ps')
PPPSop = pylops.VStack((PPop, PSop))

# data
dPPPS = PPPSop * m.flatten()
dPPPS = dPPPS.reshape(2, nt0, ntheta)

dPPPSn = dPPPS + np.random.normal(0, 1e-2, dPPPS.shape)

# Invert
minvPPSP, dPPPS_res = \
    pylops.avo.prestack.PrestackInversion(dPPPS, theta,
                                          [wav, 2 * wav], m0=mback,
                                          linearization=['fatti', 'ps'],
                                          epsR=5e-1, returnres=True,
                                          **dict(damp=1e-1, iter_lim=100))

# Data and model
fig, (axd, axdn, axvp, axvs, axrho) = \
    plt.subplots(1, 5, figsize=(8, 5), sharey=True)
axd.imshow(dPPPSn[0], cmap='gray',
           extent=(theta[0], theta[-1], t0[-1], t0[0]),
           vmin=-np.abs(dPPPSn[0]).max(), vmax=np.abs(dPPPSn[0]).max())
axd.set_title('PP Data')
axd.axis('tight')
axdn.imshow(dPPPSn[1], cmap='gray',
            extent=(theta[0], theta[-1], t0[-1], t0[0]),
            vmin=-np.abs(dPPPSn[1]).max(), vmax=np.abs(dPPPSn[1]).max())
axdn.set_title('PS Data')
axdn.axis('tight')
axvp.plot(vp, t0, 'k', lw=4, label='True')
axvp.plot(np.exp(mback[:, 0]), t0, '--r', lw=4, label='Back')
axvp.plot(np.exp(minv_noise[:, 0]), t0, '--g', lw=2, label='PP')
axvp.plot(np.exp(minvPPSP[:, 0]), t0, '--b', lw=2, label='PP+PS')
axvp.set_title(r'$V_P$')
axvs.plot(vs, t0, 'k', lw=4, label='True')
axvs.plot(np.exp(mback[:, 1]), t0, '--r', lw=4, label='Back')
axvs.plot(np.exp(minv_noise[:, 1]), t0, '--g', lw=2, label='PP')
axvs.plot(np.exp(minvPPSP[:, 1]), t0, '--b', lw=2, label='PP+PS')
axvs.set_title(r'$V_S$')
axrho.plot(rho, t0, 'k', lw=4, label='True')
axrho.plot(np.exp(mback[:, 2]), t0, '--r', lw=4, label='Back')
axrho.plot(np.exp(minv_noise[:, 2]), t0, '--g', lw=2, label='PP')
axrho.plot(np.exp(minvPPSP[:, 2]), t0, '--b', lw=2, label='PP+PS')
axrho.set_title(r'$\rho$')
axrho.legend(loc='center left', bbox_to_anchor=(1, 0.5))
axd.axis('tight')
plt.tight_layout()

We move now to a 2d example. First of all the model is loaded and data generated.

# model
inputfile = '../testdata/avo/poststack_model.npz'

model = np.load(inputfile)
x, z = model['x'][::6]/1000., model['z'][:300]/1000.
nx, nz = len(x), len(z)
m = 1000*model['model'][:300, ::6]

mvp = m.copy()
mvs = m/2
mrho = m/3+300
m = np.log(np.stack((mvp, mvs, mrho), axis=1))

# smooth model
nsmoothz, nsmoothx = 30, 25
mback = filtfilt(np.ones(nsmoothz)/float(nsmoothz), 1, m, axis=0)
mback = filtfilt(np.ones(nsmoothx)/float(nsmoothx), 1, mback, axis=2)

# dense operator
PPop_dense = \
    pylops.avo.prestack.PrestackLinearModelling(wav, theta, vsvp=vsvp,
                                                nt0=nz, spatdims=(nx,),
                                                linearization='akirich',
                                                explicit=True)

# lop operator
PPop = pylops.avo.prestack.PrestackLinearModelling(wav, theta, vsvp=vsvp,
                                                   nt0=nz, spatdims=(nx,),
                                                   linearization='akirich')

# data
dPP = PPop_dense*m.swapaxes(0, 1).flatten()
dPP = dPP.reshape(ntheta, nz, nx).swapaxes(0, 1)
dPPn = dPP + np.random.normal(0, 5e-2, dPP.shape)

Finally we perform the same 4 different inversions as in the post-stack tutorial (see 07. Post-stack inversion for more details).

# dense inversion with noise-free data
minv_dense = \
    pylops.avo.prestack.PrestackInversion(dPP, theta, wav, m0=mback,
                                          explicit=True,
                                          simultaneous=False)

# dense inversion with noisy data
minv_dense_noisy = \
    pylops.avo.prestack.PrestackInversion(dPPn, theta, wav, m0=mback,
                                          explicit=True, epsI=4e-2,
                                          simultaneous=False)

# spatially regularized lop inversion with noisy data
minv_lop_reg = \
    pylops.avo.prestack.PrestackInversion(dPPn, theta, wav,
                                          m0=minv_dense_noisy,
                                          explicit=False, epsR=1e1,
                                          **dict(damp=np.sqrt(1e-4),
                                                 iter_lim=20))

# blockiness promoting inversion with noisy data
minv_blocky = \
    pylops.avo.prestack.PrestackInversion(dPPn, theta, wav,
                                          m0=mback,
                                          explicit=False,
                                          epsR=0.4, epsRL1=0.1,
                                          **dict(mu=0.1,
                                                 niter_outer=3,
                                                 niter_inner=3,
                                                 iter_lim=5, damp=1e-3))

Let’s now visualize the inverted elastic parameters for the different scenarios

def plotmodel(axs, m, x, z, vmin, vmax,
              params=('VP', 'VS', 'Rho'),
              cmap='gist_rainbow', title=None):
    """Quick visualization of model
    """
    for ip, param in enumerate(params):
        axs[ip].imshow(m[:, ip],
                       extent=(x[0], x[-1], z[-1], z[0]),
                       vmin=vmin, vmax=vmax, cmap=cmap)
        axs[ip].set_title('%s - %s' %(param, title))
        axs[ip].axis('tight')
    plt.setp(axs[1].get_yticklabels(), visible=False)
    plt.setp(axs[2].get_yticklabels(), visible=False)

# data
fig = plt.figure(figsize=(8, 9))
ax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)
ax2 = plt.subplot2grid((2, 3), (1, 0))
ax3 = plt.subplot2grid((2, 3), (1, 1), sharey=ax2)
ax4 = plt.subplot2grid((2, 3), (1, 2), sharey=ax2)
ax1.imshow(dPP[:, 0], cmap='gray',
           extent=(x[0], x[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax1.vlines([x[nx//5], x[nx//2], x[4*nx//5]], ymin=z[0], ymax=z[-1],
           colors='w', linestyles='--')
ax1.set_xlabel('x [km]')
ax1.set_ylabel('z [km]')
ax1.set_title(r'Stack ($\theta$=0)')
ax1.axis('tight')
ax2.imshow(dPP[:, :, nx//5], cmap='gray',
           extent=(theta[0], theta[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax2.set_xlabel(r'$\theta$')
ax2.set_ylabel('z [km]')
ax2.set_title(r'Gather (x=%.2f)' % x[nx//5])
ax2.axis('tight')
ax3.imshow(dPP[:, :, nx//2], cmap='gray',
           extent=(theta[0], theta[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax3.set_xlabel(r'$\theta$')
ax3.set_title(r'Gather (x=%.2f)' % x[nx//2])
ax3.axis('tight')
ax4.imshow(dPP[:, :, 4*nx//5], cmap='gray',
           extent=(theta[0], theta[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax4.set_xlabel(r'$\theta$')
ax4.set_title(r'Gather (x=%.2f)' % x[4*nx//5])
ax4.axis('tight')
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)

# noisy data
fig = plt.figure(figsize=(8, 9))
ax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)
ax2 = plt.subplot2grid((2, 3), (1, 0))
ax3 = plt.subplot2grid((2, 3), (1, 1), sharey=ax2)
ax4 = plt.subplot2grid((2, 3), (1, 2), sharey=ax2)
ax1.imshow(dPPn[:, 0], cmap='gray',
           extent=(x[0], x[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax1.vlines([x[nx//5], x[nx//2], x[4*nx//5]], ymin=z[0], ymax=z[-1],
           colors='w', linestyles='--')
ax1.set_xlabel('x [km]')
ax1.set_ylabel('z [km]')
ax1.set_title(r'Noisy Stack ($\theta$=0)')
ax1.axis('tight')
ax2.imshow(dPPn[:, :, nx//5], cmap='gray',
           extent=(theta[0], theta[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax2.set_xlabel(r'$\theta$')
ax2.set_ylabel('z [km]')
ax2.set_title(r'Gather (x=%.2f)' % x[nx//5])
ax2.axis('tight')
ax3.imshow(dPPn[:, :, nx//2], cmap='gray',
           extent=(theta[0], theta[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax3.set_title(r'Gather (x=%.2f)' % x[nx//2])
ax3.set_xlabel(r'$\theta$')
ax3.axis('tight')
ax4.imshow(dPPn[:, :, 4*nx//5], cmap='gray',
           extent=(theta[0], theta[-1], z[-1], z[0]),
           vmin=-0.4, vmax=0.4)
ax4.set_xlabel(r'$\theta$')
ax4.set_title(r'Gather (x=%.2f)' % x[4*nx//5])
ax4.axis('tight')
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)

# inverted models
fig, axs = plt.subplots(6, 3, figsize=(8, 19))
fig.suptitle('Model', fontsize=12, fontweight='bold', y=0.95)
plotmodel(axs[0], m, x, z, m.min(),
          m.max(), title='True')
plotmodel(axs[1], mback, x, z, m.min(),
          m.max(), title='Back')
plotmodel(axs[2], minv_dense, x, z,
          m.min(), m.max(), title='Dense')
plotmodel(axs[3], minv_dense_noisy, x, z,
          m.min(), m.max(), title='Dense noisy')
plotmodel(axs[4], minv_lop_reg, x, z,
          m.min(), m.max(), title='Lop regularized')
plotmodel(axs[5], minv_blocky, x, z,
          m.min(), m.max(), title='Lop blocky')
plt.tight_layout()
plt.subplots_adjust(top=0.92)

fig, axs = plt.subplots(1, 3, figsize=(8, 7))
for ip, param in enumerate(['VP', 'VS', 'Rho']):
    axs[ip].plot(m[:, ip, nx//2], z, 'k', lw=4, label='True')
    axs[ip].plot(mback[:, ip, nx//2], z, '--r', lw=4, label='Back')
    axs[ip].plot(minv_dense[:, ip, nx//2], z, '--b', lw=2, label='Inv Dense')
    axs[ip].plot(minv_dense_noisy[:, ip, nx//2], z, '--m', lw=2,
                 label='Inv Dense noisy')
    axs[ip].plot(minv_lop_reg[:, ip, nx//2], z, '--g', lw=2,
                 label='Inv Lop regularized')
    axs[ip].plot(minv_blocky[:, ip, nx // 2], z, '--y', lw=2,
                 label='Inv Lop blocky')
    axs[ip].set_title(param)
    axs[ip].invert_yaxis()
axs[2].legend(loc=8, fontsize='small')

• 
• 
• 
• 

Out:

<matplotlib.legend.Legend object at 0x7fc20a390390>

While the background model m0 has been provided in all the examples so far, it is worth showing that the module pylops.avo.prestack.PrestackInversion can also produce so-called relative elastic parameters (i.e., variations from an average medium property) when the background model m0 is not available.

dminv = \
    pylops.avo.prestack.PrestackInversion(dPP, theta, wav, m0=None,
                                          explicit=True,
                                          simultaneous=False)

fig, axs = plt.subplots(1, 3, figsize=(8, 3))
plotmodel(axs, dminv, x, z, -dminv.max(), dminv.max(),
          cmap='seismic', title='relative')

fig, axs = plt.subplots(1, 3, figsize=(8, 7))
for ip, param in enumerate(['VP', 'VS', 'Rho']):
    axs[ip].plot(dminv[:, ip, nx//2], z, 'k', lw=2)
    axs[ip].set_title(param)
    axs[ip].invert_yaxis()

• 
• 

09. Multi-Dimensional Deconvolution

09. Multi-Dimensional Deconvolution

This example shows how to set-up and run the pylops.waveeqprocessing.MDD inversion using synthetic data.

import warnings
import numpy as np
import matplotlib.pyplot as plt
import pylops

from pylops.utils.tapers import taper3d
from pylops.utils.wavelets import ricker
from pylops.utils.seismicevents import makeaxis, hyperbolic2d

warnings.filterwarnings('ignore')
plt.close('all')

# sphinx_gallery_thumbnail_number = 5

Let’s start by creating a set of hyperbolic events to be used as our MDC kernel

# Input parameters
par = {'ox':-150, 'dx':10, 'nx':31,
       'oy':-250, 'dy':10, 'ny':51,
       'ot':0, 'dt':0.004, 'nt':300,
       'f0': 20, 'nfmax': 200}

t0_m = [0.2]
vrms_m = [700.]
amp_m = [1.]

t0_G = [0.2, 0.5, 0.7]
vrms_G = [800., 1200., 1500.]
amp_G = [1., 0.6, 0.5]

# Taper
tap = taper3d(par['nt'], [par['ny'], par['nx']],
              (5, 5), tapertype='hanning')

# Create axis
t, t2, x, y = makeaxis(par)

# Create wavelet
wav = ricker(t[:41], f0=par['f0'])[0]

# Generate model
m, mwav = hyperbolic2d(x, t, t0_m, vrms_m, amp_m, wav)

# Generate operator
G, Gwav = np.zeros((par['ny'], par['nx'], par['nt'])), \
          np.zeros((par['ny'], par['nx'], par['nt']))
for iy, y0 in enumerate(y):
    G[iy], Gwav[iy] = hyperbolic2d(x-y0, t, t0_G, vrms_G, amp_G, wav)
G, Gwav = G*tap, Gwav*tap

# Add negative part to data and model
m = np.concatenate((np.zeros((par['nx'], par['nt']-1)), m), axis=-1)
mwav = np.concatenate((np.zeros((par['nx'], par['nt']-1)), mwav), axis=-1)
Gwav2 = np.concatenate((np.zeros((par['ny'], par['nx'], par['nt']-1)), Gwav),
                       axis=-1)

# Define MDC linear operator
Gwav_fft = np.fft.rfft(Gwav2, 2*par['nt']-1, axis=-1)
Gwav_fft = Gwav_fft[..., :par['nfmax']]

MDCop = pylops.waveeqprocessing.MDC(Gwav_fft, nt=2 * par['nt']-1, nv=1,
                                    dt=0.004, dr=1., dtype='float32')

# Create data
d = MDCop*m.flatten()
d = d.reshape(par['ny'], 2*par['nt']-1)

Let’s display what we have so far: operator, input model, and data

fig, axs = plt.subplots(1, 2, figsize=(8, 6))
axs[0].imshow(Gwav2[int(par['ny']/2)].T, aspect='auto',
              interpolation='nearest', cmap='gray',
              vmin=-np.abs(Gwav2.max()), vmax=np.abs(Gwav2.max()),
              extent=(x.min(), x.max(), t2.max(), t2.min()))
axs[0].set_title('G - inline view', fontsize=15)
axs[0].set_xlabel(r'$x_R$')
axs[1].set_ylabel(r'$t$')
axs[1].imshow(Gwav2[:, int(par['nx']/2)].T, aspect='auto',
              interpolation='nearest', cmap='gray',
              vmin=-np.abs(Gwav2.max()), vmax=np.abs(Gwav2.max()),
              extent=(y.min(), y.max(), t2.max(), t2.min()))
axs[1].set_title('G - inline view', fontsize=15)
axs[1].set_xlabel(r'$x_S$')
axs[1].set_ylabel(r'$t$')
fig.tight_layout()

fig, axs = plt.subplots(1, 2, figsize=(8, 6))
axs[0].imshow(mwav.T, aspect='auto', interpolation='nearest', cmap='gray',
              vmin=-np.abs(mwav.max()), vmax=np.abs(mwav.max()),
              extent=(x.min(), x.max(), t2.max(), t2.min()))
axs[0].set_title(r'$m$', fontsize=15)
axs[0].set_xlabel(r'$x_R$')
axs[1].set_ylabel(r'$t$')
axs[1].imshow(d.T, aspect='auto', interpolation='nearest', cmap='gray',
              vmin=-np.abs(d.max()), vmax=np.abs(d.max()),
              extent=(x.min(), x.max(), t2.max(), t2.min()))
axs[1].set_title(r'$d$', fontsize=15)
axs[1].set_xlabel(r'$x_S$')
axs[1].set_ylabel(r'$t$')
fig.tight_layout()

• 
• 

We are now ready to feed our operator to pylops.waveeqprocessing.MDD and invert back for our input model

minv, madj, psfinv, psfadj = \
    pylops.waveeqprocessing.MDD(Gwav, d[:, par['nt'] - 1:],
                                dt=par['dt'], dr=par['dx'],
                                nfmax=par['nfmax'], wav=wav,
                                twosided=True, add_negative=True,
                                adjoint=True, psf=True,
                                dtype='complex64', dottest=False,
                                **dict(damp=1e-4, iter_lim=20, show=0))

fig = plt.figure(figsize=(8, 6))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(madj.T, aspect='auto', interpolation='nearest', cmap='gray',
           vmin=-np.abs(madj.max()), vmax=np.abs(madj.max()),
           extent=(x.min(), x.max(), t2.max(), t2.min()))
ax1.set_title('Adjoint m', fontsize=15)
ax1.set_xlabel(r'$x_V$')
axs[1].set_ylabel(r'$t$')
ax2.imshow(minv.T, aspect='auto', interpolation='nearest', cmap='gray',
           vmin=-np.abs(minv.max()), vmax=np.abs(minv.max()),
           extent=(x.min(), x.max(), t2.max(), t2.min()))
ax2.set_title('Inverted m', fontsize=15)
ax2.set_xlabel(r'$x_V$')
axs[1].set_ylabel(r'$t$')
ax3.plot(madj[int(par['nx']/2)]/np.abs(madj[int(par['nx']/2)]).max(),
         t2, 'r', lw=5)
ax3.plot(minv[int(par['nx']/2)]/np.abs(minv[int(par['nx']/2)]).max(),
         t2, 'k', lw=3)
ax3.set_ylim([t2[-1], t2[0]])
fig.tight_layout()

fig, axs = plt.subplots(1, 2, figsize=(8, 6))
axs[0].imshow(psfinv[int(par['nx']/2)].T,
              aspect='auto', interpolation='nearest',
              vmin=-np.abs(psfinv.max()), vmax=np.abs(psfinv.max()),
              cmap='gray', extent=(x.min(), x.max(), t2.max(), t2.min()))
axs[0].set_title('Inverted psf - inline view', fontsize=15)
axs[0].set_xlabel(r'$x_V$')
axs[1].set_ylabel(r'$t$')
axs[1].imshow(psfinv[:, int(par['nx']/2)].T,
              aspect='auto', interpolation='nearest',
              vmin=-np.abs(psfinv.max()), vmax=np.abs(psfinv.max()),
              cmap='gray', extent=(y.min(), y.max(), t2.max(), t2.min()))
axs[1].set_title('Inverted psf - xline view', fontsize=15)
axs[1].set_xlabel(r'$x_V$')
axs[1].set_ylabel(r'$t$')
fig.tight_layout()

• 
• 

We repeat the same procedure but this time we will add a preconditioning by means of causality_precond parameter, which enforces the inverted model to be zero in the negative part of the time axis (as expected by theory). This preconditioning will have the effect of speeding up the convergence of the iterative solver and thus reduce the computation time of the deconvolution

minvprec = pylops.waveeqprocessing.MDD(Gwav, d[:, par['nt'] - 1:],
                                       dt=par['dt'], dr=par['dx'],
                                       nfmax=par['nfmax'], wav=wav,
                                       twosided=True, add_negative=True,
                                       adjoint=False, psf=False,
                                       causality_precond=True,
                                       dtype='complex64',
                                       dottest=False,
                                       **dict(damp=1e-4, iter_lim=50, show=0))

fig = plt.figure(figsize=(8, 6))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(madj.T, aspect='auto', interpolation='nearest', cmap='gray',
           vmin=-np.abs(madj.max()), vmax=np.abs(madj.max()),
           extent=(x.min(), x.max(), t2.max(), t2.min()))
ax1.set_title('Adjoint m', fontsize=15)
ax1.set_xlabel(r'$x_V$')
axs[1].set_ylabel(r'$t$')
ax2.imshow(minvprec.T, aspect='auto', interpolation='nearest', cmap='gray',
           vmin=-np.abs(minvprec.max()), vmax=np.abs(minvprec.max()),
           extent=(x.min(), x.max(), t2.max(), t2.min()))
ax2.set_title('Inverted m', fontsize=15)
ax2.set_xlabel(r'$x_V$')
axs[1].set_ylabel(r'$t$')
ax3.plot(madj[int(par['nx']/2)]/np.abs(madj[int(par['nx']/2)]).max(),
         t2, 'r', lw=5)
ax3.plot(minvprec[int(par['nx']/2)]/np.abs(minv[int(par['nx']/2)]).max(),
         t2, 'k', lw=3)
ax3.set_ylim([t2[-1], t2[0]])
fig.tight_layout()

10. Marchenko redatuming by inversion

10. Marchenko redatuming by inversion

This example shows how to set-up and run the pylops.waveeqprocessing.Marchenko inversion using synthetic data.

# sphinx_gallery_thumbnail_number = 5
# pylint: disable=C0103
import warnings
import numpy as np
import matplotlib.pyplot as plt

from scipy.signal import convolve
from pylops.waveeqprocessing import Marchenko

warnings.filterwarnings('ignore')
plt.close('all')

Let’s start by defining some input parameters and loading the test data

# Input parameters
inputfile = '../testdata/marchenko/input.npz'

vel = 2400.0         # velocity
toff = 0.045         # direct arrival time shift
nsmooth = 10         # time window smoothing
nfmax = 1000         # max frequency for MDC (#samples)
niter = 10           # iterations

inputdata = np.load(inputfile)

# Receivers
r = inputdata['r']
nr = r.shape[1]
dr = r[0, 1]-r[0, 0]

# Sources
s = inputdata['s']
ns = s.shape[1]
ds = s[0, 1]-s[0, 0]

# Virtual points
vs = inputdata['vs']

# Density model
rho = inputdata['rho']
z, x = inputdata['z'], inputdata['x']

# Reflection data (R[s, r, t]) and subsurface fields
R = inputdata['R'][:, :, :-100]
R = np.swapaxes(R, 0, 1) # just because of how the data was saved

Gsub = inputdata['Gsub'][:-100]
G0sub = inputdata['G0sub'][:-100]
wav = inputdata['wav']
wav_c = np.argmax(wav)

t = inputdata['t'][:-100]
ot, dt, nt = t[0], t[1]-t[0], len(t)

Gsub = np.apply_along_axis(convolve, 0, Gsub, wav, mode='full')
Gsub = Gsub[wav_c:][:nt]
G0sub = np.apply_along_axis(convolve, 0, G0sub, wav, mode='full')
G0sub = G0sub[wav_c:][:nt]

plt.figure(figsize=(10, 5))
plt.imshow(rho, cmap='gray', extent=(x[0], x[-1], z[-1], z[0]))
plt.scatter(s[0, 5::10], s[1, 5::10], marker='*', s=150, c='r', edgecolors='k')
plt.scatter(r[0, ::10], r[1, ::10], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(vs[0], vs[1], marker='.', s=250, c='m', edgecolors='k')
plt.axis('tight')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.title('Model and Geometry')
plt.xlim(x[0], x[-1])

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 7))
axs[0].imshow(R[0].T, cmap='gray', vmin=-1e-2, vmax=1e-2,
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[0].set_title('R shot=0')
axs[0].set_xlabel(r'$x_R$')
axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(R[ns//2].T, cmap='gray', vmin=-1e-2, vmax=1e-2,
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[1].set_title('R shot=%d' %(ns//2))
axs[1].set_xlabel(r'$x_R$')
axs[1].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)
axs[2].imshow(R[-1].T, cmap='gray', vmin=-1e-2, vmax=1e-2,
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[2].set_title('R shot=%d' %ns)
axs[2].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1.5, 0)
fig.tight_layout()

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(8, 6))
axs[0].imshow(Gsub, cmap='gray', vmin=-1e6, vmax=1e6,
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[0].set_title('G')
axs[0].set_xlabel(r'$x_R$')
axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(G0sub, cmap='gray', vmin=-1e6, vmax=1e6,
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[1].set_title('G0')
axs[1].set_xlabel(r'$x_R$')
axs[1].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)
fig.tight_layout()

• 
• 
• 

Let’s now create an object of the pylops.waveeqprocessing.Marchenko class and apply redatuming for a single subsurface point vs.

# direct arrival window
trav = np.sqrt((vs[0]-r[0])**2+(vs[1]-r[1])**2)/vel

MarchenkoWM = Marchenko(R, dt=dt, dr=dr, nfmax=nfmax, wav=wav,
                        toff=toff, nsmooth=nsmooth)

f1_inv_minus, f1_inv_plus, p0_minus, g_inv_minus, g_inv_plus = \
    MarchenkoWM.apply_onepoint(trav, G0=G0sub.T, rtm=True, greens=True,
                               dottest=True, **dict(iter_lim=niter, show=True))
g_inv_tot = g_inv_minus + g_inv_plus

Out:

Dot test passed, v^T(Opu)=405.165096 - u^T(Op^Tv)=405.165096
Dot test passed, v^T(Opu)=172.065076 - u^T(Op^Tv)=172.065076

LSQR            Least-squares solution of  Ax = b
The matrix A has 282598 rows and 282598 columns
damp = 0.00000000000000e+00   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   3.134e+07  3.134e+07    1.0e+00  3.3e-08
     1  0.00000e+00   1.374e+07  1.374e+07    4.4e-01  9.3e-01   1.1e+00  1.0e+00
     2  0.00000e+00   7.770e+06  7.770e+06    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   5.750e+06  5.750e+06    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   3.930e+06  3.930e+06    1.3e-01  3.4e-01   2.5e+00  5.1e+00
     5  0.00000e+00   3.042e+06  3.042e+06    9.7e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.423e+06  2.423e+06    7.7e-02  2.2e-01   3.3e+00  8.6e+00
     7  0.00000e+00   1.675e+06  1.675e+06    5.3e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.248e+06  1.248e+06    4.0e-02  2.0e-01   3.9e+00  1.3e+01
     9  0.00000e+00   1.004e+06  1.004e+06    3.2e-02  1.5e-01   4.2e+00  1.4e+01
    10  0.00000e+00   7.762e+05  7.762e+05    2.5e-02  1.8e-01   4.4e+00  1.6e+01

LSQR finished
The iteration limit has been reached

istop =       7   r1norm = 7.8e+05   anorm = 4.4e+00   arnorm = 6.1e+05
itn   =      10   r2norm = 7.8e+05   acond = 1.6e+01   xnorm  = 3.6e+07

We can now compare the result of Marchenko redatuming via LSQR with standard redatuming

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 7))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-5e5, vmax=5e5,
              extent=(r[0, 0], r[0, -1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$')
axs[0].set_xlabel(r'$x_R$')
axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1.2, 0)
axs[1].imshow(g_inv_minus.T, cmap='gray', vmin=-5e5, vmax=5e5,
              extent=(r[0, 0], r[0, -1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$')
axs[1].set_xlabel(r'$x_R$')
axs[1].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1.2, 0)
axs[2].imshow(g_inv_plus.T, cmap='gray', vmin=-5e5, vmax=5e5,
              extent=(r[0, 0], r[0, -1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$')
axs[2].set_xlabel(r'$x_R$')
axs[2].set_ylabel(r'$t$')
axs[2].axis('tight')
axs[2].set_ylim(1.2, 0)
fig.tight_layout()

fig = plt.figure(figsize=(12, 7))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5,
           extent=(r[0, 0], r[0, -1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$')
axs[0].set_xlabel(r'$x_R$')
axs[0].set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot.T, cmap='gray', vmin=-5e5, vmax=5e5,
           extent=(r[0, 0], r[0, -1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$')
axs[1].set_xlabel(r'$x_R$')
axs[1].set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)
ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot[nr//2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0)
fig.tight_layout()

• 
• 

Note that Marchenko redatuming can also be applied simultaneously to multiple subsurface points. Use pylops.waveeqprocessing.Marchenko.apply_multiplepoints instead of pylops.waveeqprocessing.Marchenko.apply_onepoint.

11. Radon filtering

11. Radon filtering

In this example we will be taking advantage of the pylops.signalprocessing.Radon2D operator to perform filtering of unwanted events from a seismic data. For those of you not familiar with seismic data, let’s imagine that we have a data composed of a certain number of flat events and a parabolic event , we are after a transform that allows us to separate such an event from the others and filter it out. Those of you with a geophysics background may immediately realize this is the case of seismic angle (or offset) gathers after migration and those events with parabolic moveout are generally residual multiples that we would like to suppress prior to performing further analysis of our data.

The Radon transform is actually a very good transform to perform such a separation. We can thus devise a simple workflow that takes our data as input, applies a Radon transform, filters some of the events out and goes back to the original domain.

import numpy as np
import matplotlib.pyplot as plt

import pylops
from pylops.utils.wavelets import ricker

plt.close('all')
np.random.seed(0)

Let’s first create a data composed on 3 linear events and a parabolic event.

par = {'ox':0, 'dx':2, 'nx':121,
       'ot':0, 'dt':0.004, 'nt':100,
       'f0': 30}

# linear events
v = 1500
t0 = [0.1, 0.2, 0.3]
theta = [0, 0, 0]
amp = [1., -2, 0.5]

# parabolic event
tp0 = [0.13]
px = [0]
pxx = [5e-7]
ampp = [0.7]

# create axis
taxis, taxis2, xaxis, yaxis = pylops.utils.seismicevents.makeaxis(par)

# create wavelet
wav = ricker(taxis[:41], f0=par['f0'])[0]

# generate model
y = pylops.utils.seismicevents.linear2d(xaxis, taxis, v, t0,
                                        theta, amp, wav)[1] + \
    pylops.utils.seismicevents.parabolic2d(xaxis, taxis, tp0,
                                           px, pxx, ampp, wav)[1]

We can now create the pylops.signalprocessing.Radon2D operator. We also apply its adjoint to the data to obtain a representation of those 3 linear events overlapping to a parabolic event in the Radon domain. Similarly, we feed the operator to a sparse solver like pylops.optimization.sparsity.FISTA to obtain a sparse represention of the data in the Radon domain. At this point we try to filter out the unwanted event. We can see how this is much easier for the sparse transform as each event has a much more compact representation in the Radon domain than for the adjoint transform.

# radon operator
npx = 61
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)

Rop = pylops.signalprocessing.Radon2D(taxis, xaxis, px, kind='linear',
                                      interp='nearest',  centeredh=False,
                                      dtype='float64')

# adjoint Radon transform
xadj = Rop.H * y.flatten()
xadj = xadj.reshape(npx, par['nt'])

# sparse Radon transform
xinv, niter, cost = \
    pylops.optimization.sparsity.FISTA(Rop, y.flatten(), 15,
                                       eps=1e1, returninfo=True)
xinv = xinv.reshape(npx, par['nt'])

# filtering
xfilt = np.zeros_like(xadj)
xfilt[npx//2-3:npx//2+4] = xadj[npx//2-3:npx//2+4]

yfilt = Rop * xfilt.flatten()
yfilt = yfilt.reshape(par['nx'], par['nt'])

# filtering on sparse transform
xinvfilt = np.zeros_like(xinv)
xinvfilt[npx//2-3:npx//2+4] = xinv[npx//2-3:npx//2+4]

yinvfilt = Rop * xinvfilt.flatten()
yinvfilt = yinvfilt.reshape(par['nx'], par['nt'])

Finally we visualize our results.

fig, axs = plt.subplots(1, 5, sharey=True, figsize=(12, 5))
axs[0].imshow(y.T, cmap='gray',
              vmin=-np.abs(y).max(), vmax=np.abs(y).max(),
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xadj.T, cmap='gray',
              vmin=-np.abs(xadj).max(), vmax=np.abs(xadj).max(),
              extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1].axvline(px[npx//2-3], color='r', linestyle='--')
axs[1].axvline(px[npx//2+3], color='r', linestyle='--')
axs[1].set_title('Radon')
axs[1].axis('tight')
axs[2].imshow(yfilt.T, cmap='gray',
              vmin=-np.abs(yfilt).max(), vmax=np.abs(yfilt).max(),
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[2].set_title('Filtered data')
axs[2].axis('tight')
axs[3].imshow(xinv.T, cmap='gray',
              vmin=-np.abs(xinv).max(), vmax=np.abs(xinv).max(),
              extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[3].axvline(px[npx//2-3], color='r', linestyle='--')
axs[3].axvline(px[npx//2+3], color='r', linestyle='--')
axs[3].set_title('Sparse Radon')
axs[3].axis('tight')
axs[4].imshow(yinvfilt.T, cmap='gray',
              vmin=-np.abs(y).max(), vmax=np.abs(y).max(),
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[4].set_title('Sparse filtered data')
axs[4].axis('tight')

Out:

(0.0, 240.0, 0.396, 0.0)

As expected, the Radon domain is a suitable domain for this type of filtering and the sparse transform improves the ability to filter out parabolic events with small curvature.

On the other hand, it is important to note that we have not been able to correctly preserve the amplitudes of each event. This is because the sparse Radon transform can only identify a sparsest response that explain the data within a certain threshold. For this reason a more suitable approach for preserving amplitudes could be to apply a parabolic Raodn transform with the aim of reconstructing only the unwanted event and apply an adaptive subtraction between the input data and the reconstructed unwanted event.

12. Seismic regularization

12. Seismic regularization

The problem of seismic data regularization (or interpolation) is a very simple one to write, yet ill-posed and very hard to solve.

The forward modelling operator is a simple pylops.Restriction operator which is applied along the spatial direction(s).

\[\mathbf{y} = \mathbf{R} \mathbf{x}\]

Here \(\mathbf{y} = [\mathbf{y}_{R1}^T, \mathbf{y}_{R2}^T,..., \mathbf{y}_{RN^T}]^T\) where each vector \(\mathbf{y}_{Ri}\) contains all time samples recorded in the seismic data at the specific receiver \(R_i\). Similarly, \(\mathbf{x} = [\mathbf{x}_{r1}^T, \mathbf{x}_{r2}^T,..., \mathbf{x}_{rM}^T]\), contains all traces at the regularly and finely sampled receiver locations \(r_i\).

By inverting such an equation we can create a regularized data with densely and regularly spatial direction(s).

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve

import pylops
from pylops.utils.wavelets import ricker
from pylops.utils.seismicevents import makeaxis, linear2d

np.random.seed(0)
plt.close('all')

Let’s start by creating a very simple 2d data composed of 3 linear events input parameters

par = {'ox':0, 'dx':2, 'nx':70,
       'ot':0, 'dt':0.004, 'nt':80,
       'f0': 20}

v = 1500
t0_m = [0.1, 0.2, 0.28]
theta_m = [0, 30, -80]
phi_m = [0]
amp_m = [1., -2, 0.5]

# axis
taxis, t2, xaxis, y = makeaxis(par)

# wavelet
wav = ricker(taxis[:41], f0=par['f0'])[0]

# model
_, x = linear2d(xaxis, taxis, v, t0_m, theta_m, amp_m, wav)

We can now define the spatial locations along which the data has been sampled. In this specific example we will assume that we have access only to 40% of the ‘original’ locations.

perc_subsampling = 0.6
nxsub = int(np.round(par['nx']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['nx']))[:nxsub])

# restriction operator
Rop = pylops.Restriction(par['nx']*par['nt'],
                         iava, dims=(par['nx'], par['nt']),
                         dir=0, dtype='float64')

# data
y = Rop*x.ravel()
y = y.reshape(nxsub, par['nt'])

# mask
ymask = Rop.mask(x.flatten())

# inverse
xinv = Rop / y.ravel()
xinv = xinv.reshape(par['nx'], par['nt'])

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(5, 4))
axs[0].imshow(x.T, cmap='gray', vmin=-2, vmax=2,
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2,
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[1].set_title('Masked model')
axs[1].axis('tight')

Out:

(0.0, 138.0, 0.316, 0.0)

As we can see, inverting the restriction operator is not possible without adding any prior information into the inverse problem. In the following we will consider two possible routes:

• regularized inversion with second derivative along the spatial axis

  \[J = ||\mathbf{y} - \mathbf{R} \mathbf{x}||_2 + \epsilon_\nabla ^2 ||\nabla \mathbf{x}||_2\]

• sparsity-promoting inversion with pylops.FFT2 operator used as sparsyfing transform

  \[J = ||\mathbf{y} - \mathbf{R} \mathbf{F}^H \mathbf{x}||_2 + \epsilon ||\mathbf{F}^H \mathbf{x}||_1\]

# smooth inversion
D2op = pylops.SecondDerivative(par['nx']*par['nt'],
                               dims=(par['nx'], par['nt']),
                               dir=0, dtype='float64')

xsmooth, _, _ = \
    pylops.waveeqprocessing.SeismicInterpolation(y, par['nx'], iava,
                                                 kind='spatial',
                                                 **dict(epsRs=[np.sqrt(0.1)],
                                                        damp=np.sqrt(1e-4),
                                                        iter_lim=50, show=0))

# sparse inversion with FFT2
nfft = 2**8
FFTop = pylops.signalprocessing.FFT2D(dims=[par['nx'], par['nt']],
                                      nffts=[nfft, nfft],
                                      sampling=[par['dx'], par['dt']])
X = FFTop*x.flatten()
X = np.reshape(X, (nfft, nfft))

xl1, Xl1, cost = \
    pylops.waveeqprocessing.SeismicInterpolation(y, par['nx'], iava, kind='fk',
                                                 nffts=(nfft, nfft),
                                                 sampling=(par['dx'],
                                                           par['dt']),
                                                 **dict(niter=50, eps=1e-1,
                                                        returninfo=True))

fig, axs = plt.subplots(1, 4, sharey=True, figsize=(13, 4))
axs[0].imshow(x.T, cmap='gray', vmin=-2, vmax=2,
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2,
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[1].set_title('Masked model')
axs[1].axis('tight')
axs[2].imshow(xsmooth.T, cmap='gray', vmin=-2, vmax=2,
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[2].set_title('Smoothed model')
axs[2].axis('tight')
axs[3].imshow(xl1.T, cmap='gray', vmin=-2, vmax=2,
              extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[3].set_title('L1 model')
axs[3].axis('tight')

fig, axs = plt.subplots(1, 3, figsize=(10, 2))
axs[0].imshow(np.fft.fftshift(np.abs(X[:, :nfft//2-1]), axes=0).T,
              extent=(np.fft.fftshift(FFTop.f1)[0],
                      np.fft.fftshift(FFTop.f1)[-1],
                      FFTop.f2[nfft // 2 - 1], FFTop.f2[0]))
axs[0].set_title('Model in f-k domain')
axs[0].axis('tight')
axs[0].set_xlim(-0.1, 0.1)
axs[0].set_ylim(50, 0)
axs[1].imshow(np.fft.fftshift(np.abs(Xl1[:, :nfft // 2 - 1]), axes=0).T,
              extent=(np.fft.fftshift(FFTop.f1)[0],
                      np.fft.fftshift(FFTop.f1)[-1],
                      FFTop.f2[nfft // 2 - 1], FFTop.f2[0],))
axs[1].set_title('Reconstructed model in f-k domain')
axs[1].axis('tight')
axs[1].set_xlim(-0.1, 0.1)
axs[1].set_ylim(50, 0)
axs[2].plot(cost, 'k', lw=3)
axs[2].set_title('FISTA convergence')

• 
• 

Out:

Text(0.5, 1.0, 'FISTA convergence')

We see how adding prior information to the inversion can help improving the estimate of the regularized seismic data. Nevertheless, in both cases the reconstructed data is not perfect. A better sparsyfing transform could in fact be chosen here to be the linear pylops.signalprocessing.Radon2D transform in spite of the pylops.FFT2 transform.

npx = 40
pxmax = 1e-3
px = np.linspace(-pxmax, pxmax, npx)
Radop = pylops.signalprocessing.Radon2D(taxis, xaxis, px, engine='numba')

RRop = Rop*Radop

# adjoint
Xadj_fromx = Radop.H*x.flatten()
Xadj_fromx = Xadj_fromx.reshape(npx, par['nt'])

Xadj = RRop.H*y.flatten()
Xadj = Xadj.reshape(npx, par['nt'])

# L1 inverse
xl1, Xl1, cost = \
    pylops.waveeqprocessing.SeismicInterpolation(y, par['nx'], iava,
                                                 kind='radon-linear',
                                                 spataxis=xaxis,
                                                 taxis=taxis, paxis=px,
                                                 centeredh=True,
                                                 **dict(niter=50, eps=1e-1,
                                                        returninfo=True))

fig, axs = plt.subplots(2, 3, sharey=True, figsize=(12, 7))
axs[0][0].imshow(x.T, cmap='gray', vmin=-2, vmax=2,
                 extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0][0].set_title('Data', fontsize=12)
axs[0][0].axis('tight')
axs[0][1].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2,
                 extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0][1].set_title('Masked data', fontsize=12)
axs[0][1].axis('tight')
axs[0][2].imshow(xl1.T, cmap='gray', vmin=-2, vmax=2,
                 extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0][2].set_title('Reconstructed data', fontsize=12)
axs[0][2].axis('tight')
axs[1][0].imshow(Xadj_fromx.T, cmap='gray', vmin=-70, vmax=70,
                 extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1][0].set_title('Adj. Radon on data', fontsize=12)
axs[1][0].axis('tight')
axs[1][1].imshow(Xadj.T, cmap='gray', vmin=-50, vmax=50,
                 extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1][1].set_title('Adj. Radon on subsampled data', fontsize=12)
axs[1][1].axis('tight')
axs[1][2].imshow(Xl1.T, cmap='gray', vmin=-0.2, vmax=0.2,
                 extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1][2].set_title('Inverse Radon on subsampled data', fontsize=12)
axs[1][2].axis('tight')

Out:

(-0.001, 0.001, 0.316, 0.0)

Finally, let’s take now a more realistic dataset. We will use once again the linear pylops.signalprocessing.Radon2D transform but we will take advantnge of the pylops.signalprocessing.Sliding2D operator to perform such a transform locally instead of globally to the entire dataset.

inputfile = '../testdata/marchenko/input.npz'
inputdata = np.load(inputfile)

x = inputdata['R'][50, :, ::2]
x = x/np.abs(x).max()
taxis, xaxis = inputdata['t'][::2], inputdata['r'][0]

par = {}
par['nx'], par['nt'] = x.shape
par['dx'] = inputdata['r'][0, 1] - inputdata['r'][0, 0]
par['dt'] = inputdata['t'][1] - inputdata['t'][0]

# add wavelet
wav = inputdata['wav'][::2]
wav_c = np.argmax(wav)
x = np.apply_along_axis(convolve, 1, x, wav, mode='full')
x = x[:, wav_c:][:, :par['nt']]

# gain
gain = np.tile((taxis**2)[:, np.newaxis], (1, par['nx'])).T
x = x*gain

# subsampling locations
perc_subsampling = 0.5
Nsub = int(np.round(par['nx']*perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(par['nx']))[:Nsub])

# restriction operator
Rop = pylops.Restriction(par['nx']*par['nt'], iava,
                         dims=(par['nx'], par['nt']),
                         dir=0, dtype='float64')

y = Rop*x.flatten()
xadj = Rop.H*y.flatten()

y = y.reshape(Nsub, par['nt'])
xadj = xadj.reshape(par['nx'], par['nt'])

# apply mask
ymask = Rop.mask(x.flatten())

# sliding windows with radon transform
dx = par['dx']
nwins = 4
nwin = 27
nover = 3
npx = 31
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)
dimsd = x.shape
dims = (nwins*npx, dimsd[1])

Op = \
    pylops.signalprocessing.Radon2D(taxis, np.linspace(-par['dx']*nwin//2,
                                                       par['dx']*nwin//2,
                                                       nwin),
                                    px, centeredh=True, kind='linear',
                                    engine='numba')
Slidop = pylops.signalprocessing.Sliding2D(Op, dims, dimsd, nwin, nover,
                                           tapertype='cosine', design=True)

# adjoint
RSop = Rop*Slidop

Xadj_fromx = Slidop.H*x.flatten()
Xadj_fromx = Xadj_fromx.reshape(npx*nwins, par['nt'])

Xadj = RSop.H*y.flatten()
Xadj = Xadj.reshape(npx*nwins, par['nt'])

# inverse
xl1, Xl1, _ = \
    pylops.waveeqprocessing.SeismicInterpolation(y, par['nx'], iava,
                                                 kind='sliding',
                                                 spataxis=xaxis,
                                                 taxis=taxis, paxis=px,
                                                 nwins=nwins, nwin=nwin,
                                                 nover=nover,
                                                 **dict(niter=50, eps=1e-2))

fig, axs = plt.subplots(2, 3, sharey=True, figsize=(12, 14))
axs[0][0].imshow(x.T, cmap='gray', vmin=-0.1, vmax=0.1,
                 extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0][0].set_title('Data')
axs[0][0].axis('tight')
axs[0][1].imshow(ymask.T, cmap='gray', vmin=-0.1, vmax=0.1,
                 extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0][1].set_title('Masked data')
axs[0][1].axis('tight')
axs[0][2].imshow(xl1.T, cmap='gray', vmin=-0.1, vmax=0.1,
                 extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]))
axs[0][2].set_title('Reconstructed data')
axs[0][2].axis('tight')
axs[1][0].imshow(Xadj_fromx.T, cmap='gray', vmin=-1, vmax=1,
                 extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1][0].set_title('Adjoint Radon on data')
axs[1][0].axis('tight')
axs[1][1].imshow(Xadj.T, cmap='gray', vmin=-0.6, vmax=0.6,
                 extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1][1].set_title('Adjoint Radon on subsampled data')
axs[1][1].axis('tight')
axs[1][2].imshow(Xl1.T, cmap='gray', vmin=-0.03, vmax=0.03,
                 extent=(px[0], px[-1], taxis[-1], taxis[0]))
axs[1][2].set_title('Inverse Radon on subsampled data')
axs[1][2].axis('tight')

Out:

(-0.0005, 0.0005, 1.995, 0.0)

As expected the linear pylops.signalprocessing.Radon2D is able to locally explain events in the input data and leads to a satisfactory recovery. Note that increasing the number of iterations and sliding windows can further refine the result, especially the accuracy of weak events, as shown in this companion notebook.

13. Deghosting

13. Deghosting

Single-component seismic data can be decomposed in their up- and down-going constituents in a model driven fashion. This task can be achieved by defining an f-k propagator (or ghost model) and solving an inverse problem as described in pylops.waveeqprocessing.Deghosting.

# sphinx_gallery_thumbnail_number = 3
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import lsqr

import pylops

np.random.seed(0)
plt.close('all')

Let’s start by loading the input dataset and geometry

inputfile = '../testdata/updown/input.npz'
inputdata = np.load(inputfile)

vel_sep = 2400.0 # velocity at separation level
clip = 1e-1 # plotting clip


# Receivers
r = inputdata['r']
nr = r.shape[1]
dr = r[0, 1]-r[0, 0]

# Sources
s = inputdata['s']

# Model
rho = inputdata['rho']

# Axes
t = inputdata['t']
nt, dt = len(t), t[1]-t[0]
x, z = inputdata['x'], inputdata['z']
dx, dz = x[1] - x[0], z[1] - z[0]

# Data
p = inputdata['p'].T
p /= p.max()

fig = plt.figure(figsize=(9, 4))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=4)
ax2 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(rho, cmap='gray', extent=(x[0], x[-1], z[-1], z[0]))
ax1.scatter(r[0, ::5], r[1, ::5], marker='v', s=150, c='b', edgecolors='k')
ax1.scatter(s[0], s[1], marker='*', s=250, c='r', edgecolors='k')
ax1.axis('tight')
ax1.set_xlabel('x [m]')
ax1.set_ylabel('y [m]')
ax1.set_title('Model and Geometry')
ax1.set_xlim(x[0], x[-1])
ax1.set_ylim(z[-1], z[0])
ax2.plot(rho[:, len(x)//2], z, 'k', lw=2)
ax2.set_ylim(z[-1], z[0])
ax2.set_yticks([])

Out:

[]

To be able to deghost the input dataset, we need to remove its direct arrival. In this example we will create a mask based on the analytical traveltime of the direct arrival.

direct = np.sqrt(np.sum((s[:, np.newaxis]-r)**2, axis=0))/vel_sep

# Window
off = 0.035
direct_off = direct + off
win = np.zeros((nt, nr))
iwin = np.round(direct_off/dt).astype(np.int)
for i in range(nr):
    win[iwin[i]:, i] = 1


fig, axs = plt.subplots(1, 2, sharey=True, figsize=(8, 7))
axs[0].imshow(p.T, cmap='gray', vmin=-clip*np.abs(p).max(),
              vmax=clip*np.abs(p).max(),
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[0].plot(r[0], direct_off, 'r', lw=2)
axs[0].set_title(r'$P$')
axs[0].axis('tight')
axs[1].imshow(win * p.T, cmap='gray', vmin=-clip*np.abs(p).max(),
              vmax=clip*np.abs(p).max(),
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[1].set_title(r'Windowed $P$')
axs[1].axis('tight')
axs[1].set_ylim(1, 0)

Out:

(1.0, 0.0)

We can now perform deghosting

pup, pdown = \
    pylops.waveeqprocessing.Deghosting(p.T, nt, nr, dt, dr, vel_sep,
                                       r[1, 0] + dz, win=win,
                                       npad=5, ntaper=11, solver=lsqr,
                                       dottest=False, dtype='complex128',
                                       **dict(damp=1e-10, iter_lim=60))

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 7))
axs[0].imshow(p.T, cmap='gray', vmin=-clip * np.abs(p).max(),
              vmax=clip * np.abs(p).max(),
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[0].set_title(r'$P$')
axs[0].axis('tight')
axs[1].imshow(pup, cmap='gray', vmin=-clip * np.abs(p).max(),
              vmax=clip * np.abs(p).max(),
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[1].set_title(r'$P^-$')
axs[1].axis('tight')
axs[2].imshow(pdown, cmap='gray', vmin=-clip * np.abs(p).max(),
              vmax=clip * np.abs(p).max(),
              extent=(r[0, 0], r[0, -1], t[-1], t[0]))
axs[2].set_title(r'$P^+$')
axs[2].axis('tight')
axs[2].set_ylim(1, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr // 2], 'k', lw=2, label=r'$p$')
plt.plot(t, pup[:, nr // 2], 'r', lw=2, label=r'$p^-$')
plt.xlim(0, t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pdown[:, nr // 2], 'b', lw=2, label=r'$p^+$')
plt.plot(t, pup[:, nr // 2], 'r', lw=2, label=r'$p^-$')
plt.xlim(0, t[200])
plt.ylim(-0.2, 0.2)
plt.legend()

• 
• 
• 

Out:

<matplotlib.legend.Legend object at 0x7fc20ab9cef0>

To see more examples head over to the following notebook: notebook1.

14. Seismic wavefield decomposition

14. Seismic wavefield decomposition

Multi-component seismic data can be decomposed in their up- and down-going constituents in a purely data driven fashion. This task can be accurately achieved by linearly combining the input pressure and particle velocity data in the frequency-wavenumber described in details in pylops.waveeqprocessing.UpDownComposition2D and pylops.waveeqprocessing.WavefieldDecomposition.

In this tutorial we will consider a simple synthetic data composed of six events (three up-going and three down-going). We will first combine them to create pressure and particle velocity data and then show how we can retrieve their directional constituents both by directly combining the input data as well as by setting an inverse problem. The latter approach results vital in case of spatial aliasing, as applying simple scaled summation in the frequency-wavenumber would result in sub-optimal decomposition due to the superposition of different frequency-wavenumber pairs at some (aliased) locations.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import filtfilt

import pylops
from pylops.utils.wavelets import ricker
from pylops.utils.seismicevents import makeaxis, hyperbolic2d

np.random.seed(0)
plt.close('all')

Let’s first the input up- and down-going wavefields

par = {'ox':-220, 'dx':5, 'nx':89,
       'ot':0, 'dt':0.004, 'nt':200,
       'f0': 40}

t0_plus = np.array([0.2, 0.5, 0.7])
t0_minus = t0_plus + 0.04
vrms = np.array([1400., 1500., 2000.])
amp = np.array([1., -0.6, 0.5])
vel_sep = 1000.0 # velocity at separation level
rho_sep = 1000.0 # density at separation level

# Create axis
t, t2, x, y = makeaxis(par)

# Create wavelet
wav = ricker(t[:41], f0=par['f0'])[0]

# Create data
_, p_minus = hyperbolic2d(x, t, t0_minus, vrms, amp, wav)
_, p_plus = hyperbolic2d(x, t, t0_plus, vrms, amp, wav)

We can now combine them to create pressure and particle velocity data

critical = 1.1
ntaper = 51
nfft = 2**10

# 2d fft operator
FFTop = pylops.signalprocessing.FFT2D(dims=[par['nx'], par['nt']],
                                      nffts=[nfft, nfft],
                                      sampling=[par['dx'], par['dt']])

#obliquity factor
[Kx, F] = np.meshgrid(FFTop.f1, FFTop.f2, indexing='ij')
k = F/vel_sep
Kz = np.sqrt((k**2-Kx**2).astype(np.complex))
Kz[np.isnan(Kz)] = 0
OBL = rho_sep*(np.abs(F)/Kz)
OBL[Kz == 0] = 0

mask = np.abs(Kx) < critical*np.abs(F)/vel_sep
OBL *= mask
OBL = filtfilt(np.ones(ntaper)/float(ntaper), 1, OBL, axis=0)
OBL = filtfilt(np.ones(ntaper)/float(ntaper), 1, OBL, axis=1)

# composition operator
UPop = \
    pylops.waveeqprocessing.UpDownComposition2D(par['nt'], par['nx'],
                                                par['dt'], par['dx'],
                                                rho_sep, vel_sep,
                                                nffts=(nfft, nfft),
                                                critical=critical*100.,
                                                ntaper=ntaper,
                                                dtype='complex128')

# wavefield modelling
d = UPop * np.concatenate((p_plus.flatten(), p_minus.flatten())).flatten()
d = np.real(d.reshape(2*par['nx'], par['nt']))
p, vz = d[:par['nx']], d[par['nx']:]

# obliquity scaled vz
VZ = FFTop * vz.flatten()
VZ = VZ.reshape(nfft, nfft)

VZ_obl = OBL * VZ
vz_obl = FFTop.H*VZ_obl.flatten()
vz_obl = np.real(vz_obl.reshape(par['nx'], par['nt']))

fig, axs = plt.subplots(1, 4, figsize=(10, 5))
axs[0].imshow(p.T, aspect='auto', vmin=-1, vmax=1,
              interpolation='nearest', cmap='gray',
              extent=(x.min(), x.max(), t.max(), t.min()))
axs[0].set_title(r'$p$', fontsize=15)
axs[0].set_xlabel('x')
axs[0].set_ylabel('t')
axs[1].imshow(vz_obl.T, aspect='auto', vmin=-1, vmax=1,
              interpolation='nearest', cmap='gray',
              extent=(x.min(), x.max(), t.max(), t.min()))
axs[1].set_title(r'$v_z^{obl}$', fontsize=15)
axs[1].set_xlabel('x')
axs[1].set_ylabel('t')
axs[2].imshow(p_plus.T, aspect='auto', vmin=-1, vmax=1,
              interpolation='nearest', cmap='gray',
              extent=(x.min(), x.max(), t.max(), t.min()))
axs[2].set_title(r'$p^+$', fontsize=15)
axs[2].set_xlabel('x')
axs[2].set_ylabel('t')
axs[3].imshow(p_minus.T, aspect='auto',
              interpolation='nearest', cmap='gray',
              extent=(x.min(), x.max(), t.max(), t.min()),
              vmin=-1, vmax=1)
axs[3].set_title(r'$p^-$', fontsize=15)
axs[3].set_xlabel('x')
axs[3].set_ylabel('t')
plt.tight_layout()

Wavefield separation is first performed using the analytical expression for combining pressure and particle velocity data in the wavenumber-frequency domain

pup_sep, pdown_sep = \
    pylops.waveeqprocessing.WavefieldDecomposition(p, vz, par['nt'], par['nx'],
                                                   par['dt'], par['dx'],
                                                   rho_sep, vel_sep,
                                                   nffts=(nfft, nfft),
                                                   kind='analytical',
                                                   critical=critical*100,
                                                   ntaper=ntaper,
                                                   dtype='complex128')
fig = plt.figure(figsize=(12, 5))
axs0 = plt.subplot2grid((2, 5), (0, 0), rowspan=2)
axs1 = plt.subplot2grid((2, 5), (0, 1), rowspan=2)
axs2 = plt.subplot2grid((2, 5), (0, 2), colspan=3)
axs3 = plt.subplot2grid((2, 5), (1, 2), colspan=3)
axs0.imshow(pup_sep.T, cmap='gray', vmin=-1, vmax=1,
            extent=(x.min(), x.max(), t.max(), t.min()))
axs0.set_title(r'$p^-$ analytical')
axs0.axis('tight')
axs1.imshow(pdown_sep.T, cmap='gray', vmin=-1, vmax=1,
            extent=(x.min(), x.max(), t.max(), t.min()))
axs1.set_title(r'$p^+$ analytical')
axs1.axis('tight')
axs2.plot(t, p[par['nx']//2], 'r', lw=2, label=r'$p$')
axs2.plot(t, vz_obl[par['nx']//2], '--b', lw=2, label=r'$v_z^{obl}$')
axs2.set_ylim(-1, 1)
axs2.set_title('Data at x=%.2f' % x[par['nx']//2])
axs2.set_xlabel('t [s]')
axs2.legend()
axs3.plot(t, pup_sep[par['nx']//2], 'r', lw=2, label=r'$p^-$ ana')
axs3.plot(t, pdown_sep[par['nx']//2], '--b', lw=2, label=r'$p^+$ ana')
axs3.set_title('Separated wavefields at x=%.2f' % x[par['nx']//2])
axs3.set_xlabel('t [s]')
axs3.set_ylim(-1, 1)
axs3.legend()
plt.tight_layout()

We repeat the same exercise but this time we invert the composition operator pylops.waveeqprocessing.UpDownComposition2D

pup_inv, pdown_inv = \
    pylops.waveeqprocessing.WavefieldDecomposition(p, vz, par['nt'], par['nx'],
                                                   par['dt'], par['dx'],
                                                   rho_sep, vel_sep,
                                                   nffts=(nfft, nfft),
                                                   kind='inverse',
                                                   critical=critical*100,
                                                   ntaper=ntaper,
                                                   dtype='complex128',
                                                   **dict(damp=1e-10,
                                                          iter_lim=20))

fig = plt.figure(figsize=(12, 5))
axs0 = plt.subplot2grid((2, 5), (0, 0), rowspan=2)
axs1 = plt.subplot2grid((2, 5), (0, 1), rowspan=2)
axs2 = plt.subplot2grid((2, 5), (0, 2), colspan=3)
axs3 = plt.subplot2grid((2, 5), (1, 2), colspan=3)
axs0.imshow(pup_inv.T, cmap='gray', vmin=-1, vmax=1,
            extent=(x.min(), x.max(), t.max(), t.min()))
axs0.set_title(r'$p^-$ inverse')
axs0.axis('tight')
axs1.imshow(pdown_inv.T, cmap='gray', vmin=-1, vmax=1,
            extent=(x.min(), x.max(), t.max(), t.min()))
axs1.set_title(r'$p^+$ inverse')
axs1.axis('tight')
axs2.plot(t, p[par['nx']//2], 'r', lw=2, label=r'$p$')
axs2.plot(t, vz_obl[par['nx']//2], '--b', lw=2, label=r'$v_z^{obl}$')
axs2.set_ylim(-1, 1)
axs2.set_title('Data at x=%.2f' % x[par['nx']//2])
axs2.set_xlabel('t [s]')
axs2.legend()
axs3.plot(t, pup_inv[par['nx']//2], 'r', lw=2, label=r'$p^-$ inv')
axs3.plot(t, pdown_inv[par['nx']//2], '--b', lw=2, label=r'$p^+$ inv')
axs3.set_title('Separated wavefields at x=%.2f' % x[par['nx']//2])
axs3.set_xlabel('t [s]')
axs3.set_ylim(-1, 1)
axs3.legend()
plt.tight_layout()

The up- and down-going constituents have been succesfully separated in both cases. Finally, we use the pylops.waveeqprocessing.UpDownComposition2D operator to reconstruct the particle velocity wavefield from its up- and down-going pressure constituents

PtoVop = \
    pylops.waveeqprocessing.PressureToVelocity(par['nt'], par['nx'],
                                               par['dt'], par['dx'],
                                               rho_sep, vel_sep,
                                               nffts=(nfft, nfft),
                                               critical=critical * 100.,
                                               ntaper=ntaper,
                                               topressure=False)

vdown_rec = (PtoVop * pdown_inv.ravel()).reshape(par['nx'], par['nt'])
vup_rec = (PtoVop * pup_inv.ravel()).reshape(par['nx'], par['nt'])
vz_rec = np.real(vdown_rec - vup_rec)

fig, axs = plt.subplots(1, 3, figsize=(13, 6))
axs[0].imshow(vz.T, cmap='gray', vmin=-1e-6, vmax=1e-6,
              extent=(x.min(), x.max(), t.max(), t.min()))
axs[0].set_title(r'$vz$')
axs[0].axis('tight')
axs[1].imshow(vz_rec.T, cmap='gray', vmin=-1e-6, vmax=1e-6,
              extent=(x.min(), x.max(), t[-1], t[0]))
axs[1].set_title(r'$vz rec$')
axs[1].axis('tight')
axs[2].imshow(vz.T - vz_rec.T, cmap='gray', vmin=-1e-6, vmax=1e-6,
              extent=(x.min(), x.max(), t[-1], t[0]))
axs[2].set_title(r'$error$')
axs[2].axis('tight')

Out:

(-220.0, 220.0, 0.796, 0.0)

To see more examples, including applying wavefield separation and regularization simultaneously, as well as 3D examples, head over to the following notebooks: notebook1 and notebook2

15. Least-squares migration

15. Least-squares migration

Seismic migration is the process by which seismic data are manipulated to create an image of the subsurface reflectivity.

While traditionally solved as the adjont of the demigration operator, it is becoming more and more common to solve the underlying inverse problem in the quest for more accurate and detailed subsurface images.

Indipendently of the choice of the modelling operator (i.e., ray-based or full wavefield-based), the demigration/migration process can be expressed as a linear operator of such a kind:

\[d(\mathbf{x_r}, \mathbf{x_s}, t) = w(t) * \int_V G(\mathbf{x}, \mathbf{x_s}, t) G(\mathbf{x_r}, \mathbf{x}, t) m(\mathbf{x}) d\mathbf{x}\]

where \(m(\mathbf{x})\) is the reflectivity at every location in the subsurface, \(G(\mathbf{x}, \mathbf{x_s}, t)\) and \(G(\mathbf{x_r}, \mathbf{x}, t)\) are the Green’s functions from source-to-subsurface-to-receiver and finally \(w(t)\) is the wavelet. Ultimately, while the Green’s functions can be computed in many different ways, solving this system of equations for the reflectivity model is what we generally refer to as Least-squares migration (LSM).

In this tutorial we will consider the most simple scenario where we use an eikonal solver to compute the Green’s functions and show how we can use the pylops.waveeqprocessing.LSM operator to perform LSM.

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import lsqr

import pylops

plt.close('all')
np.random.seed(0)

To start we create a simple model with 2 interfaces

# 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. # gradient
vel = np.outer(np.ones(nx), v0 +kv*z)

# Reflectivity Model
refl = np.zeros((nx, nz))
refl[:, 30] = -1
refl[:, 50] = 0.5

# Receivers
nr = 11
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 = 10
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]

plt.figure(figsize=(10, 5))
im = plt.imshow(vel.T, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(recs[0],  recs[1], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(sources[0], sources[1], marker='*', s=150, c='r', edgecolors='k')
plt.colorbar(im)
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]')
plt.title('Velocity')
plt.xlim(x[0], x[-1])

plt.figure(figsize=(10, 5))
im = plt.imshow(refl.T, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(recs[0],  recs[1], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(sources[0], sources[1], marker='*', s=150, c='r', edgecolors='k')
plt.colorbar(im)
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]')
plt.title('Reflectivity')
plt.xlim(x[0], x[-1])

• 
• 

Out:

(0.0, 320.0)

We can now create our LSM object and invert for the reflectivity using two different solvers: scipy.sparse.linalg.lsqr (LS solution) and pylops.optimization.sparsity.FISTA (LS solution with sparse model).

nt = 651
dt = 0.004
t = np.arange(nt)*dt
wav, wavt, wavc = pylops.utils.wavelets.ricker(t[:41], f0=20)


lsm = pylops.waveeqprocessing.LSM(z, x, t, sources, recs, v0, wav, wavc,
                                  mode='analytic')

d = lsm.Demop * refl.ravel()
d = d.reshape(ns, nr, nt)

madj = lsm.Demop.H * d.ravel()
madj = madj.reshape(nx, nz)

minv = lsm.solve(d.ravel(), solver=lsqr, **dict(iter_lim=100))
minv = minv.reshape(nx, nz)

minv_sparse = lsm.solve(d.ravel(),
                        solver=pylops.optimization.sparsity.FISTA,
                        **dict(eps=1e2, niter=100))
minv_sparse = minv_sparse.reshape(nx, nz)

# demigration
dadj = lsm.Demop * madj.ravel()
dadj = dadj.reshape(ns, nr, nt)

dinv = lsm.Demop * minv.ravel()
dinv = dinv.reshape(ns, nr, nt)

dinv_sparse = lsm.Demop * minv_sparse.ravel()
dinv_sparse = dinv_sparse.reshape(ns, nr, nt)

# sphinx_gallery_thumbnail_number = 2
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
axs[0][0].imshow(refl.T, cmap='gray', vmin=-1, vmax=1)
axs[0][0].axis('tight')
axs[0][0].set_title(r'$m$')
axs[0][1].imshow(madj.T, cmap='gray', vmin=-madj.max(), vmax=madj.max())
axs[0][1].set_title(r'$m_{adj}$')
axs[0][1].axis('tight')
axs[1][0].imshow(minv.T, cmap='gray', vmin=-1, vmax=1)
axs[1][0].axis('tight')
axs[1][0].set_title(r'$m_{inv}$')
axs[1][1].imshow(minv_sparse.T, cmap='gray', vmin=-1, vmax=1)
axs[1][1].axis('tight')
axs[1][1].set_title(r'$m_{FISTA}$')

fig, axs = plt.subplots(1, 4, figsize=(10, 4))
axs[0].imshow(d[0, :, :300].T, cmap='gray',
              vmin=-d.max(), vmax=d.max())
axs[0].set_title(r'$d$')
axs[0].axis('tight')
axs[1].imshow(dadj[0, :, :300].T, cmap='gray',
              vmin=-dadj.max(), vmax=dadj.max())
axs[1].set_title(r'$d_{adj}$')
axs[1].axis('tight')
axs[2].imshow(dinv[0, :, :300].T, cmap='gray',
              vmin=-d.max(), vmax=d.max())
axs[2].set_title(r'$d_{inv}$')
axs[2].axis('tight')
axs[3].imshow(dinv_sparse[0, :, :300].T, cmap='gray',
              vmin=-d.max(), vmax=d.max())
axs[3].set_title(r'$d_{fista}$')
axs[3].axis('tight')

fig, axs = plt.subplots(1, 4, figsize=(10, 4))
axs[0].imshow(d[ns//2, :, :300].T, cmap='gray',
              vmin=-d.max(), vmax=d.max())
axs[0].set_title(r'$d$')
axs[0].axis('tight')
axs[1].imshow(dadj[ns//2, :, :300].T, cmap='gray',
              vmin=-dadj.max(), vmax=dadj.max())
axs[1].set_title(r'$d_{adj}$')
axs[1].axis('tight')
axs[2].imshow(dinv[ns//2, :, :300].T, cmap='gray',
              vmin=-d.max(), vmax=d.max())
axs[2].set_title(r'$d_{inv}$')
axs[2].axis('tight')
axs[3].imshow(dinv_sparse[ns//2, :, :300].T, cmap='gray',
              vmin=-d.max(), vmax=d.max())
axs[3].set_title(r'$d_{fista}$')
axs[3].axis('tight')

• 
• 
• 

Out:

(-0.5, 10.5, 299.5, -0.5)

This was just a short teaser, for a more advanced set of examples of 2D and 3D traveltime-based LSM head over to this notebook.

16. CT Scan Imaging

16. CT Scan Imaging

This tutorial considers a very well-known inverse problem from the field of medical imaging.

We will be using the pylops.signalprocessing.Radon2D operator to model a sinogram, which is a graphic representation of the raw data obtained from a CT scan. The sinogram is further inverted using both a L2 solver and a TV-regularized solver like Split-Bregman.

# sphinx_gallery_thumbnail_number = 2
import numpy as np
import matplotlib.pyplot as plt
import pylops

from numba import jit

plt.close('all')
np.random.seed(10)

Let’s start by loading the Shepp-Logan phantom model. We can then construct the sinogram by providing a custom-made function to the pylops.signalprocessing.Radon2D that samples parametric curves of such a type:

\[t(r,\theta; x) = tan(90-\theta)*x + r/sin(\theta)\]

where \(\theta\) is the angle between the x-axis (\(x\)) and the perpendicular to the summation line and \(r\) is the distance from the origin of the summation line.

@jit(nopython=True)
def radoncurve(x, r, theta):
    return (r - ny//2)/(np.sin(np.deg2rad(theta))+1e-15) + np.tan(np.deg2rad(90 - theta))*x  + ny//2

x = np.load('../testdata/optimization/shepp_logan_phantom.npy')
x = x / x.max()
ny, nx = x.shape

ntheta = 150
theta = np.linspace(0., 180., ntheta, endpoint=False)

RLop = \
    pylops.signalprocessing.Radon2D(np.arange(ny), np.arange(nx),
                                    theta, kind=radoncurve,
                                    centeredh=True, interp=False,
                                    engine='numba', dtype='float64')

y = RLop.H * x.T.ravel()
y = y.reshape(ntheta, ny).T

We can now first perform the adjoint, which in the medical imaging literature is also referred to as back-projection.

This is the first step of a common reconstruction technique, named filtered back-projection, which simply applies a correction filter in the frequency domain to the adjoint model.

xrec = RLop*y.T.ravel()
xrec = xrec.reshape(nx, ny).T

fig, axs = plt.subplots(1, 3, figsize=(10, 4))
axs[0].imshow(x, vmin=0, vmax=1, cmap='gray')
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(y, cmap='gray')
axs[1].set_title('Data')
axs[1].axis('tight')
axs[2].imshow(xrec, cmap='gray')
axs[2].set_title('Adjoint model')
axs[2].axis('tight')
fig.tight_layout()

Finally we take advantage of our different solvers and try to invert the modelling operator both in a least-squares sense and using TV-reg.

Dop = [
    pylops.FirstDerivative(ny * nx, dims=(nx, ny), dir=0, edge=True,
                           kind='backward', dtype=np.float),
    pylops.FirstDerivative(ny * nx, dims=(nx, ny), dir=1, edge=True,
                           kind='backward', dtype=np.float)]
D2op = pylops.Laplacian(dims=(nx, ny), edge=True, dtype=np.float)

# L2
xinv_sm = \
    pylops.optimization.leastsquares.RegularizedInversion(RLop.H,
                                                          [D2op],
                                                          y.T.flatten(),
                                                          epsRs=[1e1],
                                                          **dict(iter_lim=20))
xinv_sm = np.real(xinv_sm.reshape(nx, ny)).T

# TV
mu = 1.5
lamda = [1., 1.]
niter = 3
niterinner = 4

xinv, niter = pylops.optimization.sparsity.SplitBregman(RLop.H, Dop, y.T.flatten(), niter, niterinner,
                           mu=mu, epsRL1s=lamda, tol=1e-4, tau=1., show=False,
                           **dict(iter_lim=20, damp=1e-2))
xinv = np.real(xinv.reshape(nx, ny)).T

fig, axs = plt.subplots(1, 3, figsize=(10, 4))
axs[0].imshow(x, vmin=0, vmax=1, cmap='gray')
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(xinv_sm, vmin=0, vmax=1, cmap='gray')
axs[1].set_title('L2 Inversion')
axs[1].axis('tight')
axs[2].imshow(xinv, vmin=0, vmax=1, cmap='gray')
axs[2].set_title('TV-Reg Inversion')
axs[2].axis('tight')
fig.tight_layout()

17. Real/Complex Inversion

17. Real/Complex Inversion

In this tutorial we will discuss two equivalent approaches to the solution of inverse problems with real-valued model vector and complex-valued data vector. In other words, we consider a modelling operator \(\mathbf{A}:\mathbb{F}^m \to \mathbb{C}^n\) (which could be the case for example for the real FFT).

Mathematically speaking, this problem can be solved equivalently by inverting the complex-valued problem:

\[\mathbf{y} = \mathbf{A} \mathbf{x}\]

or the real-valued augmented system

\[\begin{split}\begin{bmatrix} Re(\mathbf{y}) \\ Im(\mathbf{y}) \end{bmatrix} = \begin{bmatrix} Re(\mathbf{A}) \\ Im(\mathbf{A}) \end{bmatrix} \mathbf{x}\end{split}\]

Whilst we already know how to solve the first problem, let’s see how we can solve the second one by taking advantage of the real method of the pylops.LinearOperator object. We will also wrap our linear operator into a pylops.MemoizeOperator which remembers the last N model and data vectors and by-passes the computation of the forward and/or adjoint pass whenever the same pair reappears. This is very useful in our case when we want to compute the real and the imag components of

import numpy as np
import matplotlib.pyplot as plt
import pylops

plt.close('all')
np.random.seed(0)

To start we create the forward problem

n = 5
x = np.arange(n) + 1.

# make A
Ar = np.random.normal(0, 1, (n, n))
Ai = np.random.normal(0, 1, (n, n))
A = Ar + 1j * Ai
Aop = pylops.MatrixMult(A, dtype=np.complex)
y = Aop @ x

Let’s check we can solve this problem using the first formulation

A1op = Aop.toreal(forw=False, adj=True)
xinv = A1op.div(y)

print('xinv=%s\n' % xinv)

Out:

xinv=[1. 2. 3. 4. 5.]

Let’s now see how we formulate the second problem

Amop = pylops.MemoizeOperator(Aop, max_neval=10)
Arop = Amop.toreal()
Aiop = Amop.toimag()

A1op = pylops.VStack([Arop, Aiop])
y1 = np.concatenate([np.real(y), np.imag(y)])
xinv1 = np.real(A1op.div(y1))

print('xinv1=%s\n' % xinv1)

Out:

xinv1=[1. 2. 3. 4. 5.]

Frequenty Asked Questions

1. Can I visualize my operator?

Yes, you can. Every operator has a method called todense that will return the dense matrix equivalent of the operaotor. Note, however, that in order to do so we need to allocate a numpy array of the size of your operator and apply the operator N times, where N is the number of columns of the operator. The allocation can be very heavy on your memory and the computation may take long time, so use it with care only for small toy examples to understand what your operator looks like. This method should however not be abused, as the reason of working with linear operators is indeed that you don’t really need to access the explicit matrix representation of an operator.

2. Can I have an older version of cupy installed in my system (**``cupy-cudaXX<8.1.0``)?** Yes. Nevertheless you need to tell PyLops that you don’t want to use its cupy backend by setting the environment variable CUPY_PYLOPS=0. Failing to do so will lead to an error when you import pylops because some of the cupyx routines that we use are not available in earlier version of cupy.

PyLops API

The Application Programming Interface (API) of PyLops can be loosely seen as composed of a stack of three main layers:

• Linear operators: building blocks for the setting up of inverse problems
• Solvers: interfaces to a variety of solvers, providing an easy way to augment an inverse problem with additional regularization and/or preconditioning term
• Applications: high-level interfaces allowing users to easily setup and solve specific problems (while hiding the non-needed details - i.e., creation and setup of linear operators and solvers).

Linear operators

Templates

LinearOperator([Op, explicit, clinear])	Common interface for performing matrix-vector products.
FunctionOperator(f, *args, **kwargs)	Function Operator.
MemoizeOperator(Op[, max_neval])	Memoize Operator.

Basic operators

MatrixMult(A[, dims, dtype])	Matrix multiplication.
Identity(N[, M, dtype, inplace])	Identity operator.
Zero(N[, M, dtype])	Zero operator.
Diagonal(diag[, dims, dir, dtype])	Diagonal operator.
Transpose(dims, axes[, dtype])	Transpose operator.
Flip(N[, dims, dir, dtype])	Flip along an axis.
Roll(N[, dims, dir, shift, dtype])	Roll along an axis.
Pad(dims, pad[, dtype])	Pad operator.
Sum(dims, dir[, dtype])	Sum operator.
Symmetrize(N[, dims, dir, dtype])	Symmetrize along an axis.
Restriction(M, iava[, dims, dir, dtype, inplace])	Restriction (or sampling) operator.
Regression(taxis, order[, dtype])	Polynomial regression.
LinearRegression(taxis[, dtype])	Linear regression.
CausalIntegration(N[, dims, dir, sampling, …])	Causal integration.
Spread(dims, dimsd[, table, dtable, fh, …])	Spread operator.
VStack(ops[, nproc, dtype])	Vertical stacking.
HStack(ops[, nproc, dtype])	Horizontal stacking.
Block(ops[, nproc, dtype])	Block operator.
BlockDiag(ops[, nproc, dtype])	Block-diagonal operator.
Kronecker(Op1, Op2[, dtype])	Kronecker operator.
Real(dims[, dtype])	Real operator.
Imag(dims[, dtype])	Imag operator.
Conj(dims[, dtype])	Complex conjugate operator.

Smoothing and derivatives

Smoothing1D(nsmooth, dims[, dir, dtype])	1D Smoothing.
Smoothing2D(nsmooth, dims[, nodir, dtype])	2D Smoothing.
FirstDerivative(N[, dims, dir, sampling, …])	First derivative.
SecondDerivative(N[, dims, dir, sampling, …])	Second derivative.
Laplacian(dims[, dirs, weights, sampling, …])	Laplacian.
Gradient(dims[, sampling, edge, dtype, kind])	Gradient.
FirstDirectionalDerivative(dims, v[, …])	First Directional derivative.
SecondDirectionalDerivative(dims, v[, …])	Second Directional derivative.

Signal processing

Convolve1D(N, h[, offset, dims, dir, dtype, …])	1D convolution operator.
Convolve2D(N, h, dims[, offset, nodir, …])	2D convolution operator.
ConvolveND(N, h, dims[, offset, dirs, …])	ND convolution operator.
Interp(M, iava[, dims, dir, kind, dtype])	Interpolation operator.
Bilinear(iava, dims[, dtype])	Bilinear interpolation operator.
FFT(dims[, dir, nfft, sampling, real, …])	One dimensional Fast-Fourier Transform.
FFT2D(dims[, dirs, nffts, sampling, dtype])	Two dimensional Fast-Fourier Transform.
FFTND(dims[, dirs, nffts, sampling, dtype])	N-dimensional Fast-Fourier Transform.
Shift(dims, shift[, dir, nfft, sampling, …])	Shift operator
DWT(dims[, dir, wavelet, level, dtype])	One dimensional Wavelet operator.
DWT2D(dims[, dirs, wavelet, level, dtype])	Two dimensional Wavelet operator.
Seislet(slopes[, sampling, level, kind, …])	Two dimensional Seislet operator.
Radon2D(taxis, haxis, pxaxis[, kind, …])	Two dimensional Radon transform.
Radon3D(taxis, hyaxis, hxaxis, pyaxis, pxaxis)	Three dimensional Radon transform.
ChirpRadon2D(taxis, haxis, pmax[, dtype])	2D Chirp Radon transform
ChirpRadon3D(taxis, hyaxis, hxaxis, pmax[, …])	3D Chirp Radon transform
Sliding1D(Op, dim, dimd, nwin, nover[, …])	1D Sliding transform operator.
Sliding2D(Op, dims, dimsd, nwin, nover[, …])	2D Sliding transform operator.
Sliding3D(Op, dims, dimsd, nwin, nover, nop)	3D Sliding transform operator.
Patch2D(Op, dims, dimsd, nwin, nover, nop[, …])	2D Patch transform operator.
Fredholm1(G[, nz, saveGt, usematmul, dtype])	Fredholm integral of first kind.

Wave-Equation processing

PressureToVelocity(nt, nr, dt, dr, rho, vel)	Pressure to Vertical velocity conversion.
UpDownComposition2D(nt, nr, dt, dr, rho, vel)	2D Up-down wavefield composition.
UpDownComposition3D(nt, nr, dt, dr, rho, vel)	3D Up-down wavefield composition.
MDC(G, nt, nv[, dt, dr, twosided, fast, …])	Multi-dimensional convolution.
PhaseShift(vel, dz, nt, freq, kx[, ky, dtype])	Phase shift operator
Demigration(z, x, t, srcs, recs, vel, wav, …)	Kirchoff Demigration operator.

Geophysical subsurface characterization

avo.AVOLinearModelling(theta[, vsvp, nt0, …])	AVO Linearized modelling.
poststack.PoststackLinearModelling(wav, nt0)	Post-stack linearized seismic modelling operator.
prestack.PrestackLinearModelling(wav, theta)	Pre-stack linearized seismic modelling operator.
prestack.PrestackWaveletModelling(m, theta, nwav)	Pre-stack linearized seismic modelling operator for wavelet.

Solvers

Basic

solver.cg(Op, y, x0[, niter, damp, tol, …])	Conjugate gradient
solver.cgls(Op, y, x0[, niter, damp, tol, …])	Conjugate gradient least squares
solver.lsqr(Op, y, x0[, damp, atol, btol, …])	LSQR

Least-squares

leastsquares.NormalEquationsInversion(Op, …)	Inversion of normal equations.
leastsquares.RegularizedInversion(Op, Regs, data)	Regularized inversion.
leastsquares.PreconditionedInversion(Op, P, data)	Preconditioned inversion.

Sparsity

sparsity.IRLS(Op, data, nouter[, threshR, …])	Iteratively reweighted least squares.
sparsity.OMP(Op, data[, niter_outer, …])	Orthogonal Matching Pursuit (OMP).
sparsity.ISTA(Op, data, niter[, eps, alpha, …])	Iterative Shrinkage-Thresholding Algorithm (ISTA).
sparsity.FISTA(Op, data, niter[, eps, …])	Fast Iterative Shrinkage-Thresholding Algorithm (FISTA).
sparsity.SPGL1(Op, data[, SOp, tau, sigma, x0])	Spectral Projected-Gradient for L1 norm.
sparsity.SplitBregman(Op, RegsL1, data[, …])	Split Bregman for mixed L2-L1 norms.

Applications

Wave-Equation processing

SeismicInterpolation(data, nrec, iava[, …])	Seismic interpolation (or regularization).
Deghosting(p, nt, nr, dt, dr, vel, zrec[, …])	Wavefield deghosting.
WavefieldDecomposition(p, vz, nt, nr, dt, …)	Up-down wavefield decomposition.
MDD(G, d[, dt, dr, nfmax, wav, twosided, …])	Multi-dimensional deconvolution.
Marchenko(R[, R1, dt, nt, dr, nfmax, wav, …])	Marchenko redatuming
LSM(z, x, t, srcs, recs, vel, wav, wavcenter)	Least-squares Migration (LSM).

Geophysical subsurface characterization

poststack.PoststackInversion(data, wav[, …])	Post-stack linearized seismic inversion.
prestack.PrestackInversion(data, theta, wav)	Pre-stack linearized seismic inversion.

PyLops Utilities

Alongside with its Linear Operators and Solvers, PyLops contains also a number of auxiliary routines performing universal tasks that are used by several operators or simply within one or more Tutorials for the preparation of input data and subsequent visualization of results.

Shared

Dot-test

dottest(Op[, nr, nc, tol, complexflag, …])

Dot test.

Estimators

trace_hutchinson(Op[, neval, batch_size, …])	Trace of linear operator using the Hutchinson method.
trace_hutchpp(Op[, neval, sampler, backend])	Trace of linear operator using the Hutch++ method.
trace_nahutchpp(Op[, neval, sampler, c1, …])	Trace of linear operator using the NA-Hutch++ method.

Scalability test

scalability_test(Op, x[, workers, forward])

Scalability test.

Others

Synthetics

seismicevents.makeaxis(par)	Create axes t, x, and y axes
seismicevents.linear2d(x, t, v, t0, theta, …)	Linear 2D events
seismicevents.parabolic2d(x, t, t0, px, pxx, …)	Parabolic 2D events
seismicevents.hyperbolic2d(x, t, t0, vrms, …)	Hyperbolic 2D events
seismicevents.linear3d(x, y, t, v, t0, …)	Linear 3D events
seismicevents.hyperbolic3d(x, y, t, t0, …)	Hyperbolic 3D events

marchenko.directwave(wav, trav, nt, dt[, …])

Analytical direct wave in acoustic media

Signal-processing

signalprocessing.convmtx(h, n)	Convolution matrix
signalprocessing.nonstationary_convmtx(H, n)	Convolution matrix from a bank of filters
signalprocessing.slope_estimate(d, dz, dx[, …])	Local slope estimation

Tapers

tapers.taper2d(nt, nmask, ntap[, tapertype])	2D taper
tapers.taper3d(nt, nmask, ntap[, tapertype])	3D taper

Wavelets

wavelets.ricker(t[, f0])	Ricker wavelet
wavelets.gaussian(t[, std])	Gaussian wavelet

Geophysicical Reservoir characterization

avo.zoeppritz_scattering(vp1, vs1, rho1, …)	Zoeppritz solution.
avo.zoeppritz_element(vp1, vs1, rho1, vp0, …)	Single element of Zoeppritz solution.
avo.zoeppritz_pp(vp1, vs1, rho1, vp0, vs0, …)	PP reflection coefficient from the Zoeppritz scattering matrix.
avo.approx_zoeppritz_pp(vp1, vs1, rho1, vp0, …)	PP reflection coefficient from the approximate Zoeppritz equation.
avo.akirichards(theta, vsvp[, n])	Three terms Aki-Richards approximation.
avo.fatti(theta, vsvp[, n])	Three terms Fatti approximation.
avo.ps(theta, vsvp[, n])	PS reflection coefficient

Implementing new operators

Users are welcome to create new operators and add them to the PyLops library.

In this tutorial, we will go through the key steps in the definition of an operator, using the pylops.Diagonal as an example. This is a very simple operator that applies a diagonal matrix to the model in forward mode and to the data in adjoint mode.

Creating the operator

The first thing we need to do is to create a new file with the name of the operator we would like to implement. Note that as the operator will be a class, we need to follow the UpperCaseCamelCase convention both for the class itself and for the filename.

Once we have created the file, we will start by importing the modules that will be needed by the operator. While this varies from operator to operator, you will always need to import the pylops.LinearOperator class, which will be used as parent class for any of our operators:

from pylops import LinearOperator

This class is a child of the scipy.sparse.linalg.LinearOperator class itself which implements the same methods of its parent class as well as an additional method for quick inversion: such method can be easily accessed by using \ between the operator and the data (e.g., A\y).

After that we define our new object:

class Diagonal(LinearOperator):

followed by a numpydoc docstring (starting with r""" and ending with """) containing the documentation of the operator. Such docstring should contain at least a short description of the operator, a Parameters section with a detailed description of the input parameters and a Notes section providing a mathematical explanation of the operator. Take a look at some of the core operators of PyLops to get a feeling of the level of details of the mathematical explanation.

We then need to create the __init__ where the input parameters are passed and saved as members of our class. While the input parameters change from operator to operator, it is always required to create three members, the first called shape with a tuple containing the dimensions of the operator in the data and model space, the second called dtype with the data type object (np.dtype) of the model and data, and the third called explicit with a boolean (True or False) identifying if the operator can be inverted by a direct solver or requires an iterative solver. This member is True if the operator has also a member A that contains the matrix to be inverted like for example in the pylops.MatrixMult operator, and it will be False otherwise. In this case we have another member called d which is equal to the input vector containing the diagonal elements of the matrix we want to multiply to the model and data.

def __init__(self, d, dtype=None):
    self.d = d.flatten()
    self.shape = (len(self.d), len(self.d))
    self.dtype = np.dtype(dtype)
    self.explicit = False

We can then move onto writing the forward mode in the method _matvec. In other words, we will need to write the piece of code that will implement the following operation \(\mathbf{y} = \mathbf{A}\mathbf{x}\). Such method is always composed of two inputs (the object itself self and the input model x). In our case the code to be added to the forward is very simple, we will just need to apply element-wise multiplication between the model \(\mathbf{x}\) and the elements along the diagonal contained in the array \(\mathbf{d}\). We will finally need to return the result of this operation:

def _matvec(self, x):
    return self.d*x

Finally we need to implement the adjoint mode in the method _rmatvec. In other words, we will need to write the piece of code that will implement the following operation \(\mathbf{x} = \mathbf{A}^H\mathbf{y}\). Such method is also composed of two inputs (the object itself self and the input data y). In our case the code to be added to the forward is the same as the one from the forward (but this will obviously be different from operator to operator):

def _rmatvec(self, x):
    return self.d*x

And that’s it, we have implemented our first linear operator!

Testing the operator

Being able to write an operator is not yet a guarantee of the fact that the operator is correct, or in other words that the adjoint code is actually the adjoint of the forward code. Luckily for us, a simple test can be performed to check the validity of forward and adjoint operators, the so called dot-test.

We can generate random vectors \(\mathbf{u}\) and \(\mathbf{v}\) and verify the the following equality within a numerical tolerance:

\[(\mathbf{A}*\mathbf{u})^H*\mathbf{v} = \mathbf{u}^H*(\mathbf{A}^H*\mathbf{v})\]

The method pylops.utils.dottest implements such a test for you, all you need to do is create a new test within an existing test_*.py file in the pytests folder (or in a new file).

Generally a test file will start with a number of dictionaries containing different parameters we would like to use in the testing of one or more operators. The test itself starts with a decorator that contains a list of all (or some) of dictionaries that will would like to use for our specific operator, followed by the definition of the test

@pytest.mark.parametrize("par", [(par1),(par2)])
def test_Diagonal(par):

At this point we can first of all create the operator and run the pylops.utils.dottest preceded by the assert command. Moreover, the forward and adjoint methods should tested towards expected outputs or even better, when the operator allows it (i.e., operator is invertible), a small inversion should be run and the inverted model tested towards the input model.

"""Dot-test and inversion for diagonal operator
"""
d = np.arange(par['nx']) + 1.

Dop = Diagonal(d)
assert dottest(Dop, par['nx'], par['nx'],
               complexflag=0 if par['imag'] == 1 else 3)

x = np.ones(par['nx'])
xlsqr = lsqr(Dop, Dop * x, damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=4)

Documenting the operator

Once the operator has been created, we can add it to the documentation of PyLops. To do so, simply add the name of the operator within the index.rst file in docs/source/api directory.

Moreover, in order to facilitate the user of your operator by other users, a simple example should be provided as part of the Sphinx-gallery of the documentation of the PyLops library. The directory examples containes several scripts that can be used as template.

Final checklist

Before submitting your new operator for review, use the following checklist to ensure that your code adheres to the guidelines of PyLops:

• you have created a new file containing a single class (or a function when the new operator is a simple combination of existing operators - see pylops.Laplacian for an example of such operator) and added to a new or existing directory within the pylops package.
• the new class contains at least __init__, _matvec and _matvec methods.
• the new class (or function) has a numpydoc docstring documenting at least the input Parameters and with a Notes section providing a mathematical explanation of the operator
• a new test has been added to an existing test_*.py file within the pytests folder. The test should verify that the new operator passes the pylops.utils.dottest. Moreover it is advisable to create a small toy example where the operator is applied in forward mode and the resulting data is inverted using \ from pylops.LinearOperator.
• the new operator is used within at least one example (in examples directory) or one tutorial (in tutorials directory).

Contributing

Contributions are welcome and greatly appreciated!

The best way to get in touch with the core developers and mantainers is to join the PyLops slack channel as well as open new Issues directly from the github repo.

Moreover, take a look at the Roadmap page for a list of current ideas for improvements and additions to the PyLops library.

Types of Contributions

Report Bugs

Report bugs at https://github.com/PyLops/pylops/issues.

If you are playing with the PyLops library and find a bug, please report it including:

• Your operating system name and version.
• Any details about your Python environment.
• Detailed steps to reproduce the bug.

Propose New Operators or Features

Open an issue at https://github.com/PyLops/pylops/issues with tag enhancement.

If you are proposing a new operator or a new feature:

• Explain in detail how it should work.
• Keep the scope as narrow as possible, to make it easier to implement.

Implement Operators or Features

Look through the Git issues for operator or feature requests. Anything tagged with enhancement is open to whoever wants to implement it.

Add Examples or improve Documentation

Writing new operators is not the only way to get involved and contribute. Create examples with existing operators as well as improving the documentation of existing operators is as important as making new operators and very much encouraged.

Getting Started to contribute

Ready to contribute?

1. Fork the PyLops repo.
2. Clone your fork locally:

>>  git clone https://github.com/your_name_here/pylops.git

3. Follow the installation instructions for developers that you find in Installation page. Ensure that you are able to pass all the tests before moving forward.
4. Add the main repository to the list of your remotes (this will be important to ensure you pull the latest changes before tyring to merge your local changes):

>>  git remote add upstream https://github.com/equinor/pylops

5. Create a branch for local development:

>>  git checkout -b name-of-your-branch

Now you can make your changes locally.

6. When you’re done making changes, check that your code follows the guidelines for Implementing new operators and that the both old and new tests pass successfully:

>>  make tests

7. Commit your changes and push your branch to GitLab:

>>  git add .
>> git commit -m "Your detailed description of your changes."
>> git push origin name-of-your-branch

Remember to add -u when pushing the branch for the first time.

8. Submit a pull request through the GitHub website.

Pull Request Guidelines

Before you submit a pull request, check that it meets these guidelines:

1. The pull request should include new tests for all the core routines that have been developed.
2. If the pull request adds functionality, the docs should be updated accordingly.
3. Ensure that the updated code passes all tests.

Changelog

Version 1.16.0

Released on: 11/12/2021

• Added pylops.utils.estimators module for trace estimation
• Added x0 in pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA to handle non-zero initial guess
• Modified pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA to handle multiple right hand sides
• Modified creation of haxis in pylops.signalprocessing.Radon2D and pylops.signalprocessing.Radon3D to allow for uncentered spatial axes
• Fixed _rmatvec for explicit in pylops.LinearOperator._ColumnLinearOperator

Version 1.15.0

Released on: 23/10/2021

• Added pylops.signalprocessing.Shift operator.
• Added option to choose derivative kind in pylops.avo.poststack.PoststackInversion and pylops.avo.prestack.PrestackInversion.
• Improved efficiency of adjoint of pylops.signalprocessing.Fredholm1 by applying complex conjugation to the vectors.
• Added vsvp to pylops.avo.prestack.PrestackInversion allowing to use user defined VS/VP ratio.
• Added kind to pylops.basicoperators.CausalIntegration allowing full, half, or trapezoidal integration.
• Fixed _hardthreshold_percentile in pylops.optimization.sparsity - Issue #249.
• Fixed r2norm in pylops.optimization.solver.cgls.

Version 1.14.0

Released on: 09/07/2021

• Added pylops.optimization.solver.lsqr solver
• Added utility routine pylops.utils.scalability_test for scalability tests when using multiprocessing
• Added pylops.avo.avo.ps AVO modelling option and restructured pylops.avo.prestack.PrestackLinearModelling to allow passing any function handle that can perform AVO modelling apart from those directly available
• Added R-linear operators (when setting the property clinear=False of a linear operator). pylops.basicoperators.Real, pylops.basicoperators.Imag, and pylops.basicoperators.Conj
• Added possibility to run operators pylops.basicoperators.HStack, pylops.basicoperators.VStack, pylops.basicoperators.Block pylops.basicoperators.BlockDiag, and pylops.signalprocessing.Sliding3D using multiprocessing
• Added dtype to vector X when using scipy.sparse.linalg.lobpcg in eigs method of pylops.LinearOperator
• Use kind=forward fot FirstDerivative in pylops.avo.poststack.PoststackInversion inversion when dealing with L1 regularized inversion as it makes the inverse problem more stable (no ringing in solution)
• Changed cost in pylops.optimization.solver.cg and pylops.optimization.solver.cgls to be L2 norms of residuals
• Fixed pylops.utils.dottest.dottest for imaginary vectors and to ensure u and v vectors are of same dtype of the operator

Version 1.13.0

Released on: 26/03/2021

• Added pylops.signalprocessing.Sliding1D and pylops.signalprocessing.Patch2D operators
• Added pylops.basicoperators.MemoizeOperator operator
• Added decay and analysis option in pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA solvers
• Added toreal and toimag methods to pylops.LinearOperator
• Make nr and nc optional in pylops.utils.dottest.dottest
• Fixed complex check in pylops.basicoperators.MatrixMult when working with complex-valued cupy arrays
• Fixed bug in data reshaping in check in pylops.avo.prestack.PrestackInversion
• Fixed loading error when using old cupy and/or cusignal (see Issue #201)

Version 1.12.0

Released on: 22/11/2020

• Modified all operators and solvers to work with cupy arrays
• Added eigs and solver submodules to pylops.optimization
• Added deps and backend submodules to pylops.utils
• Fixed bug in pylops.signalprocessing.Convolve2D. and pylops.signalprocessing.ConvolveND. when dealing with filters that have less dimensions than the input vector.

Version 1.11.1

Released on: 24/10/2020

• Fixed import of pyfttw when not available in :py:class:``pylops.signalprocessing.ChirpRadon3D`

Version 1.11.0

Released on: 24/10/2020

• Added pylops.signalprocessing.ChirpRadon2D and pylops.signalprocessing.ChirpRadon3D operators.
• Fixed bug in the inferred dimensions for regularization data creation in pylops.optimization.leastsquares.NormalEquationsInversion, pylops.optimization.leastsquares.RegularizedInversion, and pylops.optimization.sparsity.SplitBregman.
• Changed dtype of pylops.HStack to allow automatic inference from dtypes of input operator.
• Modified dtype of pylops.waveeqprocessing.Marchenko operator to ensure that outputs of forward and adjoint are real arrays.
• Reverted to previous complex-friendly implementation of pylops.optimization.sparsity._softthreshold to avoid division by 0.

Version 1.10.0

Released on: 13/08/2020

• Added tosparse method to pylops.LinearOperator.
• Added kind=linear in pylops.signalprocessing.Seislet operator.
• Added kind to pylops.FirstDerivative. operator to perform forward and backward (as well as centered) derivatives.
• Added kind to pylops.optimization.sparsity.IRLS solver to choose between data or model sparsity.
• Added possibility to use scipy.sparse.linalg.lobpcg in pylops.LinearOperator.eigs and pylops.LinearOperator.cond
• Added possibility to use scipy.signal.oaconvolve in pylops.signalprocessing.Convolve1D.
• Added NRegs to pylops.optimization.leastsquares.NormalEquationsInversion to allow providing regularization terms directly in the form of H^T H.

Version 1.9.1

Released on: 25/05/2020

• Changed internal behaviour of pylops.sparsity.OMP when niter_inner=0. Automatically reverts to Matching Pursuit algorithm.
• Changed handling of dtype in pylops.signalprocessing.FFT and pylops.signalprocessing.FFT2D to ensure that the type of the input vector is retained when applying forward and adjoint.
• Added dtype parameter to the FFT calls in the definition of the pylops.waveeqprocessing.MDD operation. This ensure that the type of the real part of G input is enforced to the output vectors of the forward and adjoint operations.

Version 1.9.0

Released on: 13/04/2020

• Added pylops.waveeqprocessing.Deghosting and pylops.signalprocessing.Seislet operators
• Added hard and half thresholds in pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA solvers
• Added prescaled input parameter to pylops.waveeqprocessing.MDC and pylops.waveeqprocessing.Marchenko
• Added sinc interpolation to pylops.signalprocessing.Interp (kind == 'sinc')
• Modified pylops.waveeqprocessing.marchenko.directwave to to model analytical responses from both sources of volume injection (derivative=False) and source of volume injection rate (derivative=True)
• Added pylops.LinearOperator.asoperator method to pylops.LinearOperator
• Added pylops.utils.signalprocessing.slope_estimate function
• Fix bug in pylops.signalprocessing.Radon2D and pylops.signalprocessing.Radon3D when onthefly=True returning the same result as when onthefly=False

Version 1.8.0

Released on: 12/01/2020

• Added pylops.LinearOperator.todense method to pylops.LinearOperator
• Added pylops.signalprocessing.Bilinear, pylops.signalprocessing.DWT, and pylops.signalprocessing.DWT2 operators
• Added pylops.waveeqprocessing.PressureToVelocity, pylops.waveeqprocessing.UpDownComposition3Doperator, and pylops.waveeqprocessing.PhaseShift operators
• Fix bug in pylops.basicoperators.Kronecker (see Issue #125)

Version 1.7.0

Released on: 10/11/2019

• Added pylops.Gradient, pylops.Sum, pylops.FirstDirectionalDerivative, and pylops.SecondDirectionalDerivative operators
• Added pylops.LinearOperator._ColumnLinearOperator private operator
• Added possibility to directly mix Linear operators and numpy/scipy 2d arrays in pylops.VStack and pylops.HStack and pylops.BlockDiag operators
• Added pylops.optimization.sparsity.OMP solver

Version 1.6.0

Released on: 10/08/2019

• Added pylops.signalprocessing.ConvolveND operator
• Added pylops.utils.signalprocessing.nonstationary_convmtx to create matrix for non-stationary convolution
• Added possibility to perform seismic modelling (and inversion) with non-stationary wavelet in pylops.avo.poststack.PoststackLinearModelling
• Create private methods for pylops.Block, pylops.avo.poststack.PoststackLinearModelling, pylops.waveeqprocessing.MDC to allow calling different operators (e.g., from pylops-distributed or pylops-gpu) within the method

Version 1.5.0

Released on: 30/06/2019

• Added conj method to pylops.LinearOperator
• Added pylops.Kronecker, pylops.Roll, and pylops.Transpose operators
• Added pylops.signalprocessing.Fredholm1 operator
• Added pylops.optimization.sparsity.SPGL1 and pylops.optimization.sparsity.SplitBregman solvers
• Sped up pylops.signalprocessing.Convolve1D using scipy.signal.fftconvolve for multi-dimensional signals
• Changes in implementation of pylops.waveeqprocessing.MDC and pylops.waveeqprocessing.Marchenko to take advantage of primitives operators
• Added epsRL1 option to pylops.avo.poststack.PoststackInversion and pylops.avo.prestack.PrestackInversion to include TV-regularization terms by means of pylops.optimization.sparsity.SplitBregman solver

Version 1.4.0

Released on: 01/05/2019

• Added numba engine to pylops.Spread and pylops.signalprocessing.Radon2D operators
• Added pylops.signalprocessing.Radon3D operator
• Added pylops.signalprocessing.Sliding2D and pylops.signalprocessing.Sliding3D operators
• Added pylops.signalprocessing.FFTND operator
• Added pylops.signalprocessing.Radon3D operator
• Added niter option to pylops.LinearOperator.eigs method
• Added show option to pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA solvers
• Added pylops.waveeqprocessing.seismicinterpolation, pylops.waveeqprocessing.waveeqdecomposition and pylops.waveeqprocessing.lsm submodules
• Added tests for engine in various operators
• Added documentation regarding usage of pylops Docker container

Version 1.3.0

Released on: 24/02/2019

• Added fftw engine to pylops.signalprocessing.FFT operator
• Added pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA sparse solvers
• Added possibility to broadcast (handle multi-dimensional arrays) to pylops.Diagonal and pylops..Restriction operators
• Added pylops.signalprocessing.Interp operator
• Added pylops.Spread operator
• Added pylops.signalprocessing.Radon2D operator

Version 1.2.0

Released on: 13/01/2019

• Added pylops.LinearOperator.eigs and pylops.LinearOperator.cond methods to estimate estimate eigenvalues and conditioning number using scipy wrapping of ARPACK
• Modified default dtype for all operators to be float64 (or complex128) to be consistent with default dtypes used by numpy (and scipy) for real and complex floating point numbers.
• Added pylops.Flip operator
• Added pylops.Symmetrize operator
• Added pylops.Block operator
• Added pylops.Regression operator performing polynomial regression and modified pylops.LinearRegression to be a simple wrapper of pylops.Regression when order=1
• Modified pylops.MatrixMult operator to work with both numpy ndarrays and scipy sparse matrices
• Added pylops.avo.prestack.PrestackInversion routine
• Added possibility to have a data weight via Weight input parameter to pylops.optimization.leastsquares.NormalEquationsInversion and pylops.optimization.leastsquares.RegularizedInversion solvers
• Added pylops.optimization.sparsity.IRLS solver

Version 1.1.0

Released on: 13/12/2018

• Added pylops.CausalIntegration operator

Version 1.0.1

Released on: 09/12/2018

• Changed module from lops to pylops for consistency with library name (and pip install).
• Removed quickplots from utilities and matplotlib from requirements of PyLops.

Version 1.0.0

Released on: 04/12/2018

• First official release.

Roadmap

This roadmap is aimed at providing an high-level overview on the bug fixes, improvements and new functionality that are planned for the PyLops library.

Any of the fixes/additions mentioned in the roadmap are directly linked to a Github Issue that provides more details onto the reason and initial thoughts for the implementation of such a fix/addition.

Striked tasks have been completed and related github issue closed with more details regarding how this task has been carried out.

Library structure

• Create a child repository and python library called geolops (just a suggestion) where geoscience-related operators and examples are moved across, keeping the core pylops library very generic and multi-purpose - Issue #22.

Code cleaning

• Change all np.flatten() into np.ravel() - Issue #24.
• Fix all if: return ... else: ... statements to enforce a single return with the same number of outputs - Issue #26.
• Protected attributes and @property attributes in linear operator classes? - Issue #27.

Code optimization

• Investigate speed-up given by decorating _matvec and _rmatvec methods with numba @jit and @stencil decorators - Issue #23.
• Replace np.fft.* routines used in several submodules with double engine, numpy and pyFFTW - Issue #20.

Modules

avo

• Add possibility to choose different damping factors for each elastic parameter to invert for in pylops.avo.prestack.PrestackInversion - Issue #25.

basicoperators

• Create Kronecker operator - Issue #28.
• Deal with edges in FirstDerivative and SecondDerivative operators - Issue #34.

optimization

• Sparse solvers - Issue #44.

signalprocessing

• Compare performance in FTT operator of performing np.swap+np.fft.fft(…, axis=-1) versus np.fft.fft(…, axis=chosen) - Issue #33.
• Add Wavelet operator performing the wavelet transform - Issue #21.
• Fredholm1 operator applying Fredholm integrals of first kind - Issue #31.
• Fredholm2 operators applying Fredholm integrals of second kind - Issue #31.

utils

Nothing so far

waveeqprocessing

• numpy.matmul as a way to speed up integral computation (i.e., inner for loop) in MDC operator - Issue #32.
• NMO operator performing NMO modelling - Issue #29.
• WavefieldDecomposition operator performing acoustic wavefield separation by inversion - Issue #30.

Papers

This section lists various conference abstracts and papers using the PyLops framework:

• Vakalis, S.Chen, D., Yan, M., and Nanzer, J. A., Image enhancement in active incoherent millimeter-wave imaging. Passive and Active Millimeter-Wave Imaging XXIV 11745 link (2021).
• Li, X., Becker, T., Ravasi, M., Robertsson, J., and van Manen D.J., Closed-aperture unbounded acoustics experimentation using multidimensional deconvolution. The Journal of the Acoustical Society of America 149 (3), 1813-1828 link (2021).
• Kuijpers, D., Vasconcelos, I., and Putzky, P., Reconstructing missing seismic data using Deep Learning. arXiv, link (2021).
• Ravasi, M., and Birnie, C., A Joint Inversion-Segmentation approach to Assisted Seismic Interpretation. arXiv, link - code (2021).
• Haindl, C., Ravasi, M., and Broggini, F., Handling gaps in acquisition geometries — Improving Marchenko-based imaging using sparsity-promoting inversion and joint inversion of time-lapse data. Geophysics, 86 (2), S143-S154 link - code (2021).
• Ulrich, I. E., Zunino, A., Boehm, C., and Fichtner, A., Sparsifying regularizations for stochastic sample average minimization in ultrasound computed tomography. Medical Imaging 2021: Ultrasonic Imaging and Tomography. International Society for Optics and Photonics, link (2021).
• R Feng, R., Mejer Hansen, T., Grana, D., and Balling, N., An unsupervised deep-learning method for porosity estimation based on poststack seismic data. Geophysics, 85 (6), M97-M105. link (2020).
• Nightingale J. W., Hayes, R.G, et al., PyAutoLens: Open-Source Strong Gravitational Lensing, JOSS, link - code (2020).
• Zhang, M., Marchenko Green’s functions from compressive sensing acquisition SEG Technical Program Expanded Abstracts, link (2020).
• Vargas, D., and Vasconcelos I., Rayleigh-Marchenko Redatuming Using Scattered Fields in Highly Complex Media, EAGE Technical Program Expanded Abstracts, link (2020).
• Ravasi, M., and Vasconcelos I., Implementation of Large-Scale Integral Operators with Modern HPC Solutions, EAGE Technical Program Expanded Abstracts, link - code (2020).
• Ruan, J., and Vasconcelos I., Data-and prior-driven sampling and wavefield reconstruction for sparse, irregularly-sampled, higher-order gradient data, SEG Technical Program Expanded Abstracts, link - code (2019).

Citing

When using pylops in scientific publications, please cite the following paper:

• Ravasi, M., and Vasconcelos I., PyLops–A linear-operator Python library for scalable algebra and optimization, Software X, (2020). link.

Contributors

• Matteo Ravasi, mrava87
• Carlos da Costa, cako
• Dieter Werthmüller, prisae
• Tristan van Leeuwen, TristanvanLeeuwen
• Leonardo Uieda, leouieda
• Filippo Broggini, filippo82
• Tyler Hughes, twhughes
• Lyubov Skopintseva, lskopintseva
• Francesco Picetti, fpicetti
• Alan Richardson, ar4
• BurningKarl, BurningKarl
• Nick Luiken, NickLuiken