Under the stationary phase approximation, the equation is in the standard linear form:
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):
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:
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:
In this case, the gradient can be derived as
Neglecting second-derivative terms, the Hessian can be derived in the Gauss-Newton sense as:
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:
where the elements along the diagonal are given by:
where the term has been added for stability. The inverse diagonal Hessian is then simply the inverse of each diagonal element of :
where is another diagonal operator with elements:
Applying the approximate inverse diagonal Hessian operator (85) to the gradient (80) results in an approximate l2 solution for the geometric reflection coefficient:
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 l2 amplitude estimate.
It should be noted that the form and weight structure of the l2 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.