Introduction

The generic geophysical inverse problem Claerbout (1998a); Tarantola (1987) can be summarized as follows - given a linear forward modeling operator ${\bf A}$, and some recorded data ${\bf d}$,estimate a model ${\bf m}$ such that ${\bf A} \, {\bf m} \approx {\bf d}$. If the system is over-determined, the model that minimizes the expected (L2) error in predicted data is given by the solution to the normal equations ():  

$${\bf m}_{L2} = ({\bf A}'\, {\bf A})^{-1} \; {\bf A}' \, {\bf d}.$$ (92)

Typically the matrices involved in industrial-scale geophysical inverse problems are too large to invert directly, and we depend on iterative gradient-based linear solvers to estimate solutions. However, operators such as prestack depth migration are so expensive to apply that we can only afford to iterate a handful of times, at best.

To attempt to speed convergence, we can always change model-space variables from ${\bf m}$ to ${\bf x}$ through a linear operator ${\bf P}$, and solve the following new system for ${\bf x}$,

$${\bf d}={\bf A \, P \, x} = {\bf B \, x}.$$ (93)

When we find a solution, we can then recover the model estimate, ${\bf m}_{L2}={\bf P} \, {\bf x}$.

If we choose the operator ${\bf P}$ such that ${\bf B}' {\bf B} \approx {\bf I}$, then even simply applying the adjoint (${\bf B}'$)will yield a good model estimate; furthermore, gradient-based solvers should converge to a solution of the new system rapidly in only a few iterations.