Linear least-squares inversion

Tarantola (1987) formalizes the geophysical inverse problem by giving a theoretical approach to compensate for experimental deficiency (e.g., acquisition geometry, complex overburden), while being consistent with the acquired data. His approach can be summarized as follows: given a linear modeling operator ${\bf L}$, compute synthetic data d, using,

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

where m is a reflectivity model. 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 reflectivity model $\hat{{\bf m}}$ that minimizes $S({\bf m})$ is given by

$$\begin{eqnarray} \hat{{\bf m}}&=&({\bf L}'{\bf L})^{-1}{\bf L}' {\bf d}_{obs} \\ \hat{{\bf m}}&=& {\bf H}^{-1} {\bf m}_{mig}, \end{eqnarray}$$ (3) (4)

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

The main difficulty with this approach is the explicit calculation of the Hessian inverse. In practice, it is more feasible to compute the least-squares inverse image as the solution of the linear system of equations,  

$${\bf H} \hat{{\bf m}}={\bf m}_{mig},$$ (5)

by using an iterative conjugate gradient algorithm.

The inversion inherent in equation 5 needs regularization. Prucha et al. (2000) and Kuehl and Sacchi (2001) propose smoothing the image in the offset ray parameter dimension, which is equivalent to the same procedure in the reflection angle dimension. This idea can be generalize to include the azimuth dimension.

The least squares solution of equation 5 is obtained using the fitting goals,

$$\begin{eqnarray} {\bf H}({\bf x,\Theta};{\bf x',\Theta'}) \hat{{\bf m}}({\bf x},... ... \\ {\bf D}(\Theta) \hat{{\bf m}}({\bf x},\Theta)&\approx& 0, \end{eqnarray}$$ (6)

where ${\Theta}=(\theta,\alpha)$ are the reflection and the azimuth angles, and ${\bf D}(\Theta)$ is a smoothing operator in the reflection and azimuth angle dimensions.

The next sections show how to include the subsurface offset dimension in the Hessian computation and how to go from subsurface offset to reflection and azimuth angle dimensions following the Sava and Fomel (2003) approach.