Model Roughness

Spatial Gradient Roughness Operator

The gradient roughness operator output by Mamba2D in the sect_penalty_file approximates the integral of the model spatial gradient over each parameter using the formula

\[\| \mathbf{R} \mathbf{m} \|^{2} = \sum_{i=1}^{m} \left ( \sum_{j=1}^{N(i)} \left [ w_j \left( m_i - m_j \right) \right ]^2 \right )\]

where \(N(i)\) is the number of neighboring parameters around parameter \(\mathbf{m}_i\) and \(w_j\) are weights set so that the sum of differences approximates the integral of the local gradient. Note that the parameter indices correspond to the parameter numbers in the Resistivity File (.resistivity), not the region numbers since each region can have fixed or free parameters, and because each region will have more than one parameter for anisotropic models. See [Key16] for more details.

Minimum Gradient Support (MGS) Roughness Operator

Minimum gradient support regularization can used to find models that have relatively sharper gradients between blocky regions of nearly constant conductivity (e.g., [PoZh99]). This is accomplished by using roughness penalty relaxation weights that are inversely proportional to gradients in the previous model iteration \(\mathbf{m}_{i-1}\)

\[\begin{split}\begin{eqnarray} \mathbf{q} &=& \mathbf{R}\mathbf{m}_{i-1} \\ \mathbf{W}_{mgs} &=& \left[ \mathrm{diag}(\mathbf{q}^T \mathbf{q}) + \mathbf{I}\beta_{mgs}^2 \right ]^{-1}. \end{eqnarray}\end{split}\]

The complete model roughness norm is updated to use this diagonal weighting matrix via

\[\| \mathbf{W}_{mgs} \mathbf{R} \mathbf{m} \|^{2}.\]

This is applied iteratively with a new MGS weighting matrix computed for each iteration using the model vector from the previous iteration. Since the elements of \(|\mathbf{q}|\) will be large where the model has large gradients, the MGS weights (which are the inverse of the square of this quantity) will be small at these locations and hence relax the model roughness penalty. This creates a self-reinforcing relaxation scheme that iteratively allows the model to become rougher and rougher at these locations as the inversion iterations proceed.

The \(\beta_{mgs}\) term is used to avoid division by zero, but also serves as a threshold that that small values of \(|\mathbf{q}|\) will have MGS weights equal to unity, and locations with larger \(|\mathbf{q}|\) values will have proportionally smaller MGS weights. Trial and error is need to find a good value for this parameter.

Because the MGS weights can greatly reduce the model regularization constraint, they can greatly destabilize the inversion and their application must be done with care. One strategy is to perform a normal smooth inversion to find an acceptable model, then turn on the MGS regularization for a few final model polishing iterations.

See [BKS12] for more details and a variant of this that uses gradients in seismic models for the term in the denominator.

Anisotropy Roughness Operator

For anisotropic models, the roughness is augmented by splitting the model vector into anisotropic subsets

\[\begin{split}\mathbf{m} = \left [ \begin{matrix} \mathbf{m}_x \\ \mathbf{m}_y \\ \mathbf{m}_z \end{matrix} \right]\end{split}\]

so that

\[\| \mathbf{R} \mathbf{m} \|^{2} \equiv \| \mathbf{R} \mathbf{m}_x \|^{2} + \| \mathbf{R} \mathbf{m}_y \|^{2} + \| \mathbf{R}\mathbf{m}_z \|^{2} + \alpha \| \mathbf{m}-\mathbf{m}' \|^{2}\]

where the last term on the right is used to penalize anisotropy and can be arbitrarily dialed up or down with the scalar parameter \(\alpha\) (also referred to as the anisotropy penalty weight), and the permutation of model vector

\[\begin{split}\mathbf{m}' = \left [ \begin{matrix} \mathbf{m}_y \\ \mathbf{m}_z \\ \mathbf{m}_x \end{matrix} \right].\end{split}\]

For the case of transversely isotropic models, there are only two anisotropic components instead of the three components shown above.