Reflectivity estimation

Under the stationary phase approximation, the $\grave{P}\!\acute{P}$ equation is in the standard linear form:

 

$${\bf d}= {\bf F}{\bf m}\;,$$ (76)

where ${\bf m}=\grave{P}\!\acute{P}$, ${\bf d}=D$, and F is the forward theoretical mapping of $\grave{P}\!\acute{P}$ to U, given (66). This linear form 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}\;,$$ (77)

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 inner product:

 

$$< {\bf d}, {\bf F}{\bf m}\gt \; = \; < {\bf F}'{\bf d}, {\bf m}\gt \;.$$ (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:

 

$${\bf F}'\{\cdot \} = \int_{\omega} \int_{{\bf x}_m} \rho\ome... ... \{ \cdot \} e^{-i\omega\tau_{sr}} \, d{\bf x}_m\, d\omega\;.$$ (79)

In this case, the gradient can be derived as

 

$${\bf g}({\bf x};{\bf x}_h) = - {\bf F}'{\bf d}= - \int_{\o... ...mega^2 A_s A_r D e^{-i\omega\tau_{sr}} \,d{\bf x}_m\,d\omega\;.$$ (80)

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

 

$$\H({\bf x},{\bf x}^{^{\prime}};{\bf x}_h) = {\bf F}'{\bf F... ... dV' \right] e^{-i\omega\tau_{sr}} \, d{\bf x}_m\, d\omega\;.$$ (81)

Evaluation of one ``column'' of the Hessian operator requires a complete forward modeling response to an impulse located at a single subsurface location ${\bf x}^{^{\prime}}$, followed by a complete prestack depth migration of that impulse response back to all subsurface positions ${\bf x}$. Therefore, the complete evaluation of the full Hessian operator requires $N_{{\bf x}}$ forward modeling runs plus $N_{{\bf x}}$ prestack depth migrations, where $N_{{\bf x}}$ 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 ${\bf x}$ due to ${\bf x}^{^{\prime}}$. 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 ($\grave{P}\!\acute{P}$), and it is likely to be spatially correlated due to finite data bandwidth. The Hessian $\H$ 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 $\grave{P}\!\acute{P}$), by assuming the ``'' in the volume integrand of (81) is a spatially uncorrelated delta function $\delta({\bf x}^{^{\prime}}-{\bf x})$ (e.g., Le Bras and Clayton, 1988). In this case, we find a diagonal operator ${\bf h}({\bf x},{\bf x}^{^{\prime}};{\bf x}_h)$ which is the diagonalized Hessian such that:

 

$${\bf h}({\bf x},{\bf x}^{^{\prime}};{\bf x}_h) \approx \H({\bf x},{\bf x}^{^{\prime}};{\bf x}_h) \;,$$ (82)

where the elements along the diagonal are given by:

 

$${\bf h}({\bf x}={\bf x}^{^{\prime}};{\bf x}_h) = \H({\bf x}=... ..._r \, d{\bf x}_m\, d\omega+ \epsilon^2({\bf x};{\bf x}_h) \;,$$ (83)

where the $\epsilon^2$ term has been added for stability. The inverse diagonal Hessian is then simply the inverse of each diagonal element of ${\bf h}({\bf x},{\bf x}^{^{\prime}})$:

 

$$\H^{-1}({\bf x},{\bf x}^{^{\prime}};{\bf x}_h) \approx {\bf h}^{-1}({\bf x},{\bf x}^{^{\prime}};{\bf x}_h) \;,$$ (84)

where ${\bf h}^{-1}$ is another diagonal operator with elements:

 

$$\mbox{Diag}\{{\bf h}^{-1}\} = \left( \int_{\omega} \int_{{... ... x}_m\, d\omega+ \epsilon^2({\bf x};{\bf x}_h) \right)^{-1} \;.$$ (85)

Applying the approximate inverse diagonal Hessian operator (85) to the gradient (80) results in an approximate l₂ solution for the $\grave{P}\!\acute{P}$ geometric reflection coefficient:

 

$$\grave{P}\!\acute{P}({\bf x};{\bf x}_h) = -\H^{-1}{\bf g}\ap... ...A^2_r \,d{\bf x}_m\,d\omega+ \epsilon^2({\bf x};{\bf x}_h)} \;.$$ (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 ${\bf g}$ and ${\bf h}$ 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₂ $\grave{P}\!\acute{P}(\bar{\theta})$ amplitude estimate.

It should be noted that the form and weight structure of the l₂ 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 $\grave{P}\!\acute{P}(\bar{\theta})$ amplitude result in general.