Under the stationary phase approximation, the equation is in the standard linear form:

(76) |

where , , and *F* is the
forward theoretical mapping of to *U*, given (66).
This linear form has a well-known damped least-squares
Gauss-Newton solution (e.g., Gill et al., 1981):

(77) |

where and are the gradient and Hessian of *E* respectively, and
is the adjoint operator to . The adjoint operator is
*defined* by the inner product:

(78) |

Writing down the integral equation implied in (78), interchanging the order of integration, and taking the complex conjugate, results in an integral expression for the adjoint operator:

(79) |

In this case, the gradient can be derived as

(80) |

Neglecting second-derivative terms, the Hessian can be derived in the Gauss-Newton sense as:

(81) |

Evaluation of one ``column'' of the Hessian operator requires a complete
forward modeling response to an impulse located at a *single* subsurface
location , followed by a complete prestack depth migration of that
impulse response back to *all* subsurface positions .
Therefore, the complete evaluation of the full Hessian operator requires
forward modeling runs plus prestack depth migrations,
where is the total number of discrete surface points under
consideration. This would be extremely expensive, if not impossible to
compute!

In this form, the Hessian is spatially non-diagonal in the sense that it
tries to incorporate and compensate for all point-to-point correlations
of reflectivity in the subsurface, at due to .
This non-diagonality should not be
confused with *multi-parameter* crosscorrelations which occur in the
case of acoustic or elastic parameter inversions
(e.g., Beydoun and Mendes, 1989). There is only a single parameter being
estimated here (), and it is likely to be spatially correlated due
to finite data bandwidth.
The Hessian in (81) is trying to *deconvolve*
that spatial correlation. The effect of applying the full Hessian is to
perform a *least-squares acquisition geometry and aperture compensation*
to the reflectivity estimate.

In the interest of efficiency, the Hessian operator can be diagonalized (i.e., ignoring spatial correlations of ), by assuming the ``'' in the volume integrand of (81) is a spatially uncorrelated delta function (e.g., Le Bras and Clayton, 1988). In this case, we find a diagonal operator which is the diagonalized Hessian such that:

(82) |

where the elements along the diagonal are given by:

(83) |

where the term has been added for stability. The inverse diagonal Hessian is then simply the inverse of each diagonal element of :

(84) |

where is another diagonal operator with elements:

(85) |

Applying the approximate inverse diagonal Hessian operator (85)
to the gradient (80) results in an approximate *l _{2}*
solution for the geometric reflection coefficient:

(86) |

The approximate
diagonal Hessian solution (86) costs little more than a single
conventional prestack migration in terms of floating point operations.
However, one drawback is that it requires twice the memory, since the
and images must be stored separately during the estimation
(migration) process. They may only be combined as in (86) after
both have been completely imaged. This seems to be a small price to pay
for a reasonably good *l _{2}* amplitude estimate.

It should be noted that the form and weight structure
of the *l _{2}* inverse solution (86) is somewhat similar, but different
in detail, to the trial solution of Bleistein (1987). In
particular, (86) involves two separate but simultaneous images
to be evaluated, the gradient (numerator) being a weighted Kirchhoff
prestack depth migration, and the approximate inverse Hessian (denominator)
being an accumulation of the squared migration weights.
Bleistein's solution involves only
one migration integral since it is not a least-squares estimate, and
his different migration kernel amplitude weights will give a different
amplitude result in general.

11/16/1997