SimPEG.regularization.PGI#
- class SimPEG.regularization.PGI(mesh, gmmref, alpha_x=None, alpha_y=None, alpha_z=None, alpha_xx=0.0, alpha_yy=0.0, alpha_zz=0.0, gmm=None, wiresmap=None, maplist=None, alpha_pgi=1.0, approx_hessian=True, approx_gradient=True, approx_eval=True, weights_list=None, non_linear_relationships: bool = False, reference_model_in_smooth: bool = False, **kwargs)[source]#
Bases:
ComboObjectiveFunction
Regularization function for petrophysically guided inversion (PGI).
PGI
is used to recover models in which 1) the physical property values are consistent with petrophysical information and 2) structures in the recovered model are geologically plausible.PGI
regularization is a weighted sum ofPGIsmallness
,SmoothnessFirstOrder
andSmoothnessSecondOrder
(optional) regularization functions. The PGI smallness term assumes the statistical distribution of physical property values defining the model is characterized by a Gaussian mixture model (GMM). And the smoothness terms penalize large spatial derivatives in the recovered model.PGI
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:
- mesh
SimPEG.regularization.RegularizationMesh
,discretize.base.BaseMesh
Mesh on which the regularization is discretized. Implemented for tensor, QuadTree or Octree meshes.
- gmmref
SimPEG.utils.WeightedGaussianMixture
Reference Gaussian mixture model.
- gmm
None
,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. IfNone
, thedirectives.PGI_UpdateParameters
directive must be used to ensure there is a Gaussian mixture model for the inversion.- alpha_pgi
float
Scaling constant for the PGI smallness term.
- alpha_x, alpha_y, alpha_z
float
orNone
,optional
Scaling constants for the first order smoothness along x, y and z, respectively. If set to
None
, the scaling constant is set automatically according to the value of the length_scale parameter.- alpha_xx, alpha_yy, alpha_zz0,
float
Scaling constants for the second order smoothness along x, y and z, respectively. If set to
None
, the scaling constant is set automatically according to the length scales; seeregularization.WeightedLeastSquares
.- wiresmap
None
,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.- maplist
None
,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 anmaps.IdentityMap
from model values to physical property values.- 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.
- reference_model_in_smoothbool,
optional
Whether to include the reference model in the smoothness terms.
- mesh
Notes
For one or more physical property types (e.g. conductivity, density, susceptibility), the
PGI
regularization function (objective function) is derived by using a Gaussian mixture model (GMM) to construct 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) forPGI
is given by:\[\begin{split}\phi (\mathbf{m}) &= \alpha_\text{pgi} \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 ] \\ &+ \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,
\(\mathbf{G_x, \, G_y, \; G_z}\) are partial cell gradients operators along x, y and z,
\(\mathbf{L_x, \, L_y, \; L_z}\) are second-order derivative operators with respect to x, y and z,
\(\mathbf{W}(\Theta , \mathbf{z}^\ast )\) is the weighting matrix for PGI smallness, and
\(\mathbf{W_x, \, W_y, \; W_z}\) are weighting matrices for smoothness terms.
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 isTrue
, 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:\[\begin{split}\phi (\mathbf{m}) &= \alpha_\text{pgi} \Big \| \mathbf{W}_{\! 1/2}(\Theta, \mathbf{z}^\ast ) \, \big [ \mathbf{m} - \mathbf{m_{ref}}(\Theta, \mathbf{z}^\ast ) \big ] \, \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}\]When the
approx_eval
property isTrue
, you may also set theapprox_gradient
property toTrue
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 isTrue
, 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 thegmm
property set asNone
, the GMM is learned and updated after every model update. That is, we assume the GMM defined using thegmmref
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
Full weighting matrix for the combo objective function.
Scaling constant for the PGI smallness term.
Gaussian mixture model.
Ordered list of mappings from model values to physical property values.
Mapping from the model to the quantity evaluated in the object function.
reference_model.mref has been deprecated.
Multipliers for the objective functions.
Number of model parameters.
Reference model.
Whether to include the reference model in the smoothness objective functions.
Regularization mesh.
Mapping from the model to the model parameters of each type.
Methods
__call__
(m[, f])Evaluate the objective functions for a given model.
Compute and return quasi geology model.
deriv
(m[, f])Gradient of the objective function evaluated for the model provided.
deriv2
(m[, v, f])Hessian of the objective function evaluated for the model provided.
get_functions_of_type
(fun_class)Return objective functions of a given type(s).
map_class
alias of
IdentityMap
membership
(m)Compute and return membership array for the model provided.
test
([x, num])Run a convergence test on both the first and second derivatives.
Galleries and Tutorials using SimPEG.regularization.PGI
#
Petrophysically guided inversion (PGI): Linear example
Petrophysically guided inversion: Joint linear example with nonlinear relationships
Joint PGI of Gravity + Magnetic on an Octree mesh using full petrophysical information
Joint PGI of Gravity + Magnetic on an Octree mesh without petrophysical information