Linear least-squares inversion

Tarantola (1987) formalizes the geophysical inverse problem by giving a theoretical approach to compensate for the experiment's deficiencies (acquisition geometry, obstacles, etc.), while being consistent with the acquired data. His approach can be summarized as follows: given a linear modeling operator ${\bf L}$ to compute synthetic data d,

$${\bf d}={\bf L}{\bf m},$$ (1)

where m is a reflectivity model, and given the recorded data ${\bf d}_{obs}$, a quadratic cost function

$$S({\bf m})=\Vert {\bf d} - {\bf d}_{obs} \Vert^2 =\Vert {\bf L}{\bf m} - {\bf d}_{obs} \Vert^2$$ (2)

is formed. The model of the earth $\hat{{\bf m}}$ that minimize $S({\bf m})$ is given by  

$$\hat{{\bf m}}=({\bf L}'{\bf L})^{-1}{\bf L}' {\bf d}_{obs} = {\bf H}^{-1}{\bf L}' {\bf d}_{obs},$$ (3)

where ${\bf L}'$ is the adjoint of the linear operator ${\bf L}$, and ${\bf H}={\bf L}'{\bf L}$ is the Hessian of $S({\bf m})$.

The main difficulty with this approach is that the explicit calculation of inverse of the Hessian for the entire model space is practically unfeasible. That is why iterative algorithms like conjugate-gradient have been used to implicitly calculate the inverse of the Hessian Duquet and Marfurt (1999); Kuehl and Sacchi (2001); Nemeth et al. (1999); Prucha et al. (2000); Ronen and Liner (2000).

In the case of wave-equation migration or inversion, the operator ${\bf L}$ is expensive to apply. Thus, applying this operator and its transpose iteratively is sometimes prohibitive. Among other factors, the computational cost is proportional to the number of depth steps the wavefields need to be propagated Audebert (1994), and the number of iterations.