Shaping regularization in geophysical estimation problems

Review of Tikhonov's regularization

If the data are represented by vector $\mathbf{d}$, model parameters by vector $\mathbf{m}$, and their functional relationship is defined by the forward modeling operator $\mathbf{L}$, the least-squares optimization approach amounts to minimizing the least-squares norm of the residual difference $\mathbf{L\,m - d}$. In Tikhonov's regularization approach, one additionally attempts to minimize the norm of $\mathbf{D\,m}$, where $\mathbf {D}$ is the regularization operator. Thus, we are looking for the model $\mathbf{m}$ that minimizes the least-squares norm of the compound vector $\left[\begin{array}{cc} \mathbf{L m - d} & \epsilon \mathbf{D m} \end{array}\right]^T$, where $\epsilon$ is a scalar scaling parameter. The formal solution has the well-known form

$$\widehat{\mathbf{m}} = \left(\mathbf{L}^T\,\mathbf{L} + ... ...bf{D}^T\,\mathbf{D}\right)^{-1}\,\mathbf{L}^T\,\mathbf{d}\;,$$ (1)

where $\widehat\mathbf{m}$ denotes the least-squares estimate of $\mathbf{m}$, and $\mathbf{L}^T$ denotes the adjoint operator. One can carry out the optimization iteratively with the help of the conjugate-gradient method (Hestenes and Steifel, 1952) or its analogs. Iterative methods have computational advantages in large-scale problems when forward and adjoint operators are represented by sparse matrices and can be computed efficiently (Saad, 2003; van der Vorst, 2003).

In an alternative approach, one obtains the regularized estimate by minimizing the least-squares norm of the compound vector $\left[\begin{array}{cc} \mathbf{p} & \mathbf{r} \end{array}\right]^T$ under the constraint

$$\epsilon \mathbf{r = d - L m = d - L P p}\;.$$ (2)

Here $\mathbf{P}$ is the model reparameterization operator that translates vector $\mathbf{p}$ into the model vector $\mathbf{m}$, $\mathbf{r}$ is the scaled residual vector, and $\epsilon$ has the same meaning as before. The formal solution of the preconditioned problem is given by

$$\widehat{\mathbf{m}} = \mathbf{P}\,\widehat{\mathbf{p}} =... ...f{L}^T + \epsilon^2\,\mathbf{I}\right)^{-1}\, \mathbf{d}\;,$$ (3)

where $\mathbf{I}$ is the identity operator in the data space. Estimate 3 is mathematically equivalent to estimate 1 if $\mathbf{D}^T\,\mathbf{D}$ is invertible and

$$\left(\mathbf{D}^T\,\mathbf{D}\right)^{-1} = \mathbf{P}\,\mathbf{P}^T = \mathbf{C}\;.$$ (4)

Statistical theory of least-squares estimation connects $\mathbf{C}$ with the model covariance operator (Tarantola, 2004). In a more general case of reparameterization, the size of $\mathbf{p}$ may be different from the size of $\mathbf{m}$, and $\mathbf{C}$ may not have the full rank. In iterative methods, the preconditioned formulation often leads to faster convergence. Fomel and Claerbout (2003) suggest constructing preconditioning operators in multi-dimensional problems by recursive helical filtering.

Shaping regularization in geophysical estimation problems