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)

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 the physical property distributions of each geologic units are distinct

approx_evalbool

If True, use the L2-approximation evaluation of the smallness term by assuming the physical property distributions of each geologic units are distinct

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.
gmm	Gaussian mixture model.
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.
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.
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.