pylops.avo.poststack.PoststackInversion(data, wav, m0=None, explicit=False, simultaneous=False, epsI=None, epsR=None, dottest=False, epsRL1=None, **kwargs_solver)[source]

Post-stack linearized seismic inversion.

Invert post-stack seismic operator to retrieve an elastic parameter of choice from band-limited seismic post-stack data. Depending on the choice of input parameters, inversion can be trace-by-trace with explicit operator or global with either explicit or linear operator.

data : np.ndarray

Band-limited seismic post-stack data of size \([n_{t_0}\,(\times n_x \times n_y)]\)

wav : np.ndarray

Wavelet in time domain (must have odd number of elements and centered to zero). If 1d, assume stationary wavelet for the entire time axis. If 2d of size \([n_{t_0} \times n_h]\) use as non-stationary wavelet

m0 : np.ndarray, optional

Background model of size \([n_{t_0}\,(\times n_x \times n_y)]\)

explicit : bool, optional

Create a chained linear operator (False, preferred for large data) or a MatrixMult linear operator with dense matrix (True, preferred for small data)

simultaneous : bool, optional

Simultaneously invert entire data (True) or invert trace-by-trace (False) when using explicit operator (note that the entire data is always inverted when working with linear operator)

epsI : float, optional

Damping factor for Tikhonov regularization term

epsR : float, optional

Damping factor for additional Laplacian regularization term

dottest : bool, optional

Apply dot-test

epsRL1 : float, optional

Damping factor for additional blockiness regularization term


Arbitrary keyword arguments for scipy.linalg.lstsq solver (if explicit=True and epsR=None) or scipy.sparse.linalg.lsqr solver (if explicit=False and/or epsR is not None)

minv : np.ndarray

Inverted model of size \([n_{t_0}\,(\times n_x \times n_y)]\)

datar : np.ndarray

Residual data (i.e., data - background data) of size \([n_{t_0}\,(\times n_x \times n_y)]\)


The cost function and solver used in the seismic post-stack inversion module depends on the choice of explicit, simultaneous, epsI, and epsR parameters:

  • explicit=True, epsI=None and epsR=None: the explicit solver scipy.linalg.lstsq is used if simultaneous=False (or the iterative solver scipy.sparse.linalg.lsqr is used if simultaneous=True)
  • explicit=True with epsI and epsR=None: the regularized normal equations \(\mathbf{W}^T\mathbf{d} = (\mathbf{W}^T \mathbf{W} + \epsilon_\mathbf{I}^2 \mathbf{I}) \mathbf{AI}\) are instead fed into the scipy.linalg.lstsq solver if simultaneous=False (or the iterative solver scipy.sparse.linalg.lsqr if simultaneous=True)
  • explicit=False and epsR=None: the iterative solver scipy.sparse.linalg.lsqr is used
  • explicit=False with epsR and epsRL1=None: the iterative solver pylops.optimization.leastsquares.regularized_inversion is used to solve the spatially regularized problem.
  • explicit=False with epsR and epsRL1: the iterative solver pylops.optimization.sparsity.SplitBregman is used to solve the blockiness-promoting (in vertical direction) and spatially regularized (in additional horizontal directions) problem.

Note that the convergence of iterative solvers such as scipy.sparse.linalg.lsqr can be very slow for this type of operator. It is suggested to take a two steps approach with first a trace-by-trace inversion using the explicit operator, followed by a regularized global inversion using the outcome of the previous inversion as initial guess.