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.

Attributes

W

Weighting matrix.

active_cells

Active cells defined on the regularization mesh.

cell_weights

Deprecated property for 'volume' and user defined weights.

gmm

Gaussian mixture model.

indActive

active_cells.indActive has been deprecated.

maplist

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

mapping

Mapping from the inversion model parameters to the regularization mesh.

model

The model parameters.

mref

reference_model.mref has been deprecated.

nP

Number of model parameters.

non_linear_relationships

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

parent

The parent objective function

reference_model

Reference model.

regmesh

regularization_mesh.regmesh has been deprecated.

regularization_mesh

Regularization mesh.

shape

Number of model parameters.

units

Units for the model parameters.

weights_keys

Return the keys for the existing cell weights

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.

f_m(m)

Evaluate the regularization kernel function.

f_m_deriv(m)

Derivative of the regularization kernel function.

get_weights(key)

Cell weights for a given key.

map_class

alias of IdentityMap

membership(m)

Compute and return membership array for the model provided.

remove_weights(key)

Removes the weights for the key provided.

set_weights(**weights)

Adds (or updates) the specified weights.

test([x, num, random_seed])

Run a convergence test on both the first and second derivatives.

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.

When the approx_eval property is True, we assume the physical property distributions of each geologic units are distinct (no significant overlap of their respective physical properties distribution). The GMM probability density value at any each point of the physical property space can then be approximated by the locally dominant Gaussian distribution. In this case, the PGI regularization function (objective function) can be expressed as a least-square:

\[\begin{split}\phi (\mathbf{m}) &= \alpha_\text{pgi} \big | \mathbf{W} ( \Theta , \mathbf{z}^\ast ) \, (\mathbf{m} - \mathbf{m_{ref}}(\Theta, \mathbf{z}^\ast ) \, \Big \|^2 &+ \sum_{j=x,y,z} \alpha_j \Big \| \mathbf{W_j G_j \, m} \, \Big \|^2 \\ &+ \sum_{j=x,y,z} \alpha_{jj} \Big \| \mathbf{W_{jj} L_j \, m} \, \Big \|^2 \;\;\;\;\;\;\;\; \big ( \textrm{optional} \big )\end{split}\]

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.

For the full, non-approximated PGI regularization, please refer to (Astic, et al 2019; Astic et al 2020).

When the approx_eval property is True, you may also set the approx_gradient and approx_hessian properties to True so that the least-squares approximation is used to compute the gradient, as it is making the same assumptions about the GMM.

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}\]

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}}^{\frac{-1}{2}} \, diag \big ( \mathbf{w} \big )\]

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.