Reflectivity

Under the stationary phase approximation, the $\grave{P}\!\acute{P}$ equation is in the standard linear form: ${\bf d}= {\bf F}{\bf m}$. This has a well-known damped least-squares Gauss-Newton solution (e.g., Gill et al., 1981): ${\bf m}\approx -\H^{-1}{\bf g}\approx ({\bf F}'{\bf F}+ \epsilon^2)^{-1} {\bf F}'{\bf d}$, where ${\bf g}$ and $\H$ are the gradient and Hessian of E respectively, and ${\bf F}'$ is the adjoint operator to ${\bf F}$. The adjoint operator is defined by the scalar product:

 

$$< {\bf d}, {\bf F}{\bf m}\gt \; = \; < {\bf F}'{\bf d}, {\bf m}\gt \;.$$ (11)

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

 

$${\bf F}'\{\cdot \} = - \int_w \int_{{\bf x}_m} ik\, A_s A_r \{ \cdot \} e^{-iw\tau} \, d{\bf x}_m\, dw \;.$$ (12)

In this case, the gradient can be derived as

 

$${\bf g}({\bf x};{\bf x}_h) = - {\bf F}'{\bf d}= \int_w \int_{{\bf x}_m} ik\, A_s A_r D e^{-iw\tau} \,d{\bf x}_m\,dw \;.$$ (13)

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

 

$$\H = {\bf F}'{\bf F}+ \epsilon^2 = \int_w \int_{{\bf x}_m} k... ...-iw\tau} \, d{\bf x}_m\, dw + \epsilon^2({\bf x};{\bf x}_h) \;.$$ (14)

This representation of the Hessian requires a complete forward modeling step as well as a complete migration step; i.e., it costs one full forward modeling of all the recorded data, plus one full prestack migration. 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. 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 one single parameter being inverted for here: $\grave{P}\!\acute{P}$, and it is likely to be spatially correlated due to data bandwidth, and $\H$ in (14) is trying to deconvolve that spatial correlation.

In the interest of efficiency, the Hessian can be diagonalized (i.e., ignoring spatial correlations of $\grave{P}\!\acute{P}$), by assuming the ``1'' in the volume integrand of (14) is a spatially uncorrelated delta function $\delta({\bf x}'-{\bf x})$ (e.g., Clayton and Le Bras, 1988). In this case, the diagonalized Hessian reduces to:

 

$$\H = {\bf F}'{\bf F}+ \epsilon^2 \approx \int_w \int_{{\bf x... ... A^2_r \, d{\bf x}_m\, dw + \epsilon^2({\bf x};{\bf x}_h) \;.$$ (15)

Combining (13) and (15), results in an approximate l₂ reflectivity solution under the stationary phase condition:

 

$$\grave{P}\!\acute{P}({\bf x};{\bf x}_h) \approx - \frac{ {\b... ...2_s A^2_r \,d{\bf x}_m\,dw + \epsilon^2({\bf x};{\bf x}_h)} \;.$$ (16)

The approximate diagonal Hessian solution (16) 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 ${\bf g}$ and $\H$ images must be updated separately during the estimation (migration) process. They may only be combined as in (16) after both have been completely imaged. This seems to be the price to pay for accurate $\grave{P}\!\acute{P}(\Theta)$ amplitude recovery.

It should be noted that the form and weight structure of the l₂ inverse solution (16) is somewhat similar, but different in detail, to the trial solution of Bleistein (1987). In particular, (16) involves two separate but simultaneous images to be evaluated, the gradient (numerator) being a weighted Kirchhoff prestack depth migration, and the approximate 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 $\grave{P}\!\acute{P}$ amplitude result in general.