SimPEG.regularization.PGIsmallness#

class SimPEG.regularization.PGIsmallness(gmmref, gmm=None, wiresmap=None, maplist=None, mesh=None, approx_gradient=True, approx_eval=True, approx_hessian=True, non_linear_relationships=False, **kwargs)[source]#

Bases: Smallness

Smallness regularization function for petrophysically guided inversion (PGI).

PGIsmallness is used to recover models in which the physical property values are consistent with petrophysical information. PGIsmallness regularization assumes that the statistical distribution of physical property values defining the model is characterized by a Gaussian mixture model (GMM). That is, the physical property values for each specified geological unit are characterized by a separate multivariate Gaussian distribution, which are summed to define the GMM. PGIsmallness is generally combined with other regularization classes to form a complete regularization for the inverse problem; see PGI.

PGIsmallness can be implemented to invert for a single physical property or multiple physical properties, each of which are defined on a linear scale (e.g. density) or a log-scale (e.g. electrical conductivity). If the statistical distribution(s) of physical property values for each property type are known, the GMM can be constructed and left static throughout the inversion. Otherwise, the recovered model at each iteration is used to update the GMM. And the updated GMM is used to constrain the recovered model for the following iteration.

Parameters:
gmmrefSimPEG.utils.WeightedGaussianMixture

Reference Gaussian mixture model.

gmmNone, SimPEG.utils.WeightedGaussianMixture

Set the Gaussian mixture model used to constrain the recovered physical property model. Can be left static throughout the inversion or updated using the directives.PGI_UpdateParameters directive. If None, the directives.PGI_UpdateParameters directive must be used to ensure there is a Gaussian mixture model for the inversion.

wiresmapNone, SimPEG.maps.Wires

Mapping from the model to the model parameters of each type. If None, we assume only a single physical property type in the inversion.

maplistNone, list of SimPEG.maps

Ordered list of mappings from model values to physical property values; one for each physical property. If None, we assume a single physical property type in the regularization and an maps.IdentityMap from model values to physical property values.

meshSimPEG.regularization.RegularizationMesh, discretize.base.BaseMesh

Mesh on which the regularization is discretized. Implemented for tensor, QuadTree or Octree meshes.

approx_gradientbool

If True, use the L2-approximation of the gradient by assuming physical property values of different types are uncorrelated.

approx_evalbool

If True, use the L2-approximation evaluation of the smallness term by assuming physical property values of different types are uncorrelated.

approx_hessianbool

Approximate the Hessian of the regularization function.

non_linear_relationshipbool

Whether relationships in the Gaussian mixture model are non-linear.

Notes

For one or more physical property types (e.g. conductivity, density, susceptibility), the PGIsmallness regularization function (objective function) is derived by setting a Gaussian mixture model (GMM) as the prior within a Baysian inversion scheme. For a comprehensive description, see (Astic, et al 2019; Astic et al 2020).

We let \(\Theta\) store all of the means (\(\boldsymbol{\mu}\)), covariances (\(\boldsymbol{\Sigma}\)) and proportion constants (\(\boldsymbol{\gamma}\)) defining the GMM. And let \(\mathbf{z}^\ast\) define an membership array that extracts the GMM parameters for the most representative rock unit within each active cell in the RegularizationMesh. The regularization function (objective function) for PGIsmallness is given by:

\[\phi (\mathbf{m}) = \frac{1}{2} \big [ \mathbf{m} - \mathbf{m_{ref}}(\Theta, \mathbf{z}^\ast ) \big ]^T \mathbf{W} ( \Theta , \mathbf{z}^\ast ) \, \big [ \mathbf{m} - \mathbf{m_{ref}}(\Theta, \mathbf{z}^\ast ) \big ]\]

where

  • \(\mathbf{m}\) is the model,

  • \(\mathbf{m_{ref}}(\Theta, \mathbf{z}^\ast )\) is the reference model, and

  • \(\mathbf{W}(\Theta , \mathbf{z}^\ast )\) is a weighting matrix.

PGIsmallness regularization can be used for models consisting of one or more physical property types. The ordering of the physical property types within the model is defined using the wiresmap. And the mapping from model parameter values to physical property values is specified with maplist. For \(K\) physical property types, the model is an array vector of the form:

\[\begin{split}\mathbf{m} = \begin{bmatrix} \mathbf{m}_1 \\ \mathbf{m}_2 \\ \vdots \\ \mathbf{m}_K \end{bmatrix}\end{split}\]

When the approx_eval property is True, we assume the physical property types have values that are uncorrelated. In this case, the weighting matrix is diagonal and the regularization function (objective function) can be expressed as:

\[\phi (\mathbf{m}) = \frac{1}{2} \Big \| \mathbf{W}_{\! 1/2}(\Theta, \mathbf{z}^\ast ) \, \big [ \mathbf{m} - \mathbf{m_{ref}}(\Theta, \mathbf{z}^\ast ) \big ] \, \Big \|^2\]

When the approx_eval property is True, you may also set the approx_gradient property to True so that the least-squares approximation is used to compute the gradient.

Constructing the Reference Model and Weighting Matrix:

The reference model used in the regularization function is constructed by extracting the means \(\boldsymbol{\mu}\) from the GMM using the membership array \(\mathbf{z}^\ast\). We represent this vector as:

\[\mathbf{m_{ref}} (\Theta ,{\mathbf{z}^\ast}) = \boldsymbol{\mu}_{\mathbf{z}^\ast}\]

To construct the weighting matrix, \(\mathbf{z}^\ast\) is used to extract the covariances \(\boldsymbol{\Sigma}\) for each cell. And the weighting matrix is given by:

\[\mathbf{W}(\Theta ,{\mathbf{z}^\ast } ) = \boldsymbol{\Sigma}_{\mathbf{z^\ast}}^{-1} \, diag \big ( \mathbf{v \odot w} \big )\]

where \(\mathbf{v}\) are the volumes of the active cells, and \(\mathbf{w}\) are custom cell weights. When the approx_eval property is True, the off-diagonal covariances are zero and we can use a weighting matrix of the form:

\[\mathbf{W}_{\! 1/2}(\Theta ,{\mathbf{z}^\ast } ) = diag \Big ( \big [ \mathbf{v \odot w} \odot \boldsymbol{\sigma}_{\mathbf{z}^\ast}^{-2} \big ]^{1/2} \Big )\]

where \(\boldsymbol{\sigma}_{\mathbf{z}^\ast}^2\) are the variances extracted using the membership array \(\mathbf{z}^\ast\).

Updating the Gaussian Mixture Model:

When the GMM is set using the gmm property, the GMM remains static throughout the inversion. When the gmm property set as None, the GMM is learned and updated after every model update. That is, we assume the GMM defined using the gmmref property is not completely representative of the physical property distributions for each rock unit, and we update the all of the means (\(\boldsymbol{\mu}\)), covariances (\(\boldsymbol{\Sigma}\)) and proportion constants (\(\boldsymbol{\gamma}\)) defining the GMM \(\Theta\). This is done by solving:

\[\max_\Theta \; \mathcal{P}(\Theta | \mathbf{m})\]

using a MAP variation of the expectation-maximization clustering algorithm introduced in Dempster (et al. 1977).

Updating the Membership Array:

As the model (and GMM) are updated throughout the inversion, the rock unit considered most indicative of the geology within each cell is updated; which is represented by the membership array \(\mathbf{z}^\ast\). W. For the current GMM with means (\(\boldsymbol{\mu}\)), covariances (\(\boldsymbol{\Sigma}\)) and proportion constants (\(\boldsymbol{\gamma}\)), we solve the following for each cell:

\[z_i^\ast = \max_n \; \gamma_{i,n} \, \mathcal{N} (\mathbf{m}_i | \boldsymbol{\mu}_n , \boldsymbol{\Sigma}_n)\]

where

  • \(\mathbf{m_i}\) are the model values for cell \(i\),

  • \(\gamma_{i,n}\) is the proportion for cell \(i\) and rock unit \(n\)

  • \(\boldsymbol{\mu}_n\) are the mean property values for unit \(n\),

  • \(\boldsymbol{\Sigma}_n\) are the covariances for unit \(n\), and

  • \(\mathcal{N}\) represent the multivariate Gaussian distribution.

Attributes

gmm

Gaussian mixture model.

maplist

Ordered list of mappings from model values to physical property values.

non_linear_relationships

Whether relationships in the Gaussian mixture model are non-linear.

shape

Number of model parameters.

wiresmap

Mapping from the model to the model parameters of each type.

Methods

__call__(m[, external_weights])

Evaluate the regularization function for the model provided.

compute_quasi_geology_model()

Compute and return quasi geology model.

deriv(m)

Gradient of the regularization function evaluated for the model provided.

deriv2(m[, v])

Hessian of the regularization function evaluated for the model provided.

membership(m)

Compute and return membership array for the model provided.

set_weights(**weights)

Adds (or updates) the specified weights.