Least-squares imaging and deconvolution using the hybrid norm conjugate-direction solver

Application - Least-Squares Inverse Imaging

This application originated from the work of target-oriented wave-equation LSI imaging, as described in Tang (2008); M. Clapp and Biondi (2005); Valenciano (2006). The concept of LSI imaging starts with a simple inversion problem:

$$\bf F (\bold m) = \vert\vert\bf {L} \bold m - \bold d_{obs}\vert\vert _2 ,$$ (1)

where L is a linearized wave-equation modeling operator, the adjoint of which is the imaging operator, m is the subsurface reflectivity model, and ${\bf d_{obs}}$ is the observed surface seismic data. In theory, the solution to this inversion problem can be written as follows:

$${\bf m} = {\bf H}^{+}{\bf L}'{\bf d}_{\rm obs},$$ (2)

Where $\bf {H=L'L}$ is called Hessian operator, and $\bf {H^{+}}$ is the pseudo inverse of H. In practice, it is usually impossible to invert H directly; thus a gradient-based optimization method is often used to find the solution.

One disadvantage of this data-space inversion scheme is that it can not be computed in a target-oriented way, since theoretically even a local perturbation in the model space will affect the entire data space and vice versa. To overcome this difficulty, Valenciano (2006) transformed (1) to a model space inversion based on (2):

$${\bf J}(\bold m) = \vert\vert\bold{H} \bold m - \bold{L'}\bold d_{obs}\vert\vert _2 .$$

Valenciano (2008) and Tang (2008) showed that unlike L, matrix H is usually very sparse (i.e., most of the non-zero elements are centered around the diagonal); thus despite the huge size of H, it is feasible to store an approximation of H matrix by keeping only a few off-diagonal elements without losing much accuracy.

If we write $\bold m_{mig}=\bold{L'} \bold d_{obs}$ , and add a model regularization term (since most likely H has a null space). Then the inversion formula is

$$\boldsymbol{J}(\bold m) = \vert\vert\boldsymbol{H} \bold m - \bold m_{mig}\vert\vert _2 + \epsilon \vert\vert \bold m \vert\vert _{norm} ,$$

in which we applied the hybrid norm to the regularization term.

Tang (2009) provided a way to efficiently compute the Hessian matrix using the phase-encoding technique, and this Hessian matrix is computed only once and stored for all iterations.

Least-squares imaging and deconvolution using the hybrid norm conjugate-direction solver