Reflection angle estimation

In principle, an l₂ solution for the reflection angles $\bar{\theta}$ can be estimated directly from the reflection data D by substituting the solution (86) into the normal equation $\partial E / \partial \bar{\theta}= 0$. A difficulty arises from the fact that the kernel functions of (70) are nonlinearly related to $\bar{\theta}$ by the cosine. In this case,

 

$${\bf F}{\bf m}\rightarrow {\bf F}\{\bar{\theta}\} = \int_{... ...{P}\!\acute{P} \cos(\theta_{sr}-2\bar{\theta}) \, d{\cal V}\;,$$ (87)

is a nonlinear forward mapping between the model ${\bf m}=\bar{\theta}$ and the data U.

We solve this nonlinear inverse problem for $\bar{\theta}$ by linearizing (87) with respect to a new model parameter: ${\bf m}=\cos 2\bar{\theta}$.Using the standard cosine expansion,

 

$$\cos\phi_r\equiv \cos(\theta_{sr}-2\bar{\theta}) = \cos\theta_{sr}\cos 2\bar{\theta}+ \sin\theta_{sr}\sin 2\bar{\theta}\;.$$ (88)

and using the small angle approximation:

 

$$\theta_{sr}\,, 2\bar{\theta}\ll \pi/4 \;,$$ (89)

we approximate $\cos\phi_r$ as:

 

$$\cos\phi_r\approx \cos\theta_{sr}\cos 2\bar{\theta}\;.$$ (90)

We note that under this approximation, at the specular reflection (stationary point):

 

$$\cos\phi_r= 1 \approx \cos\theta_{sr}\cos 2\bar{\theta}\;.$$ (91)

We will use this relation (91) later to ``undo'' the small angle approximation after we obtain a linearized solution, as will be discussed later.

Now, under approximation (90), we obtain a linearized forward theory relating the new model parameter ${\bf m}=\cos 2\bar{\theta}$ to the data U:

 

$${\bf F}{\bf m}= {\bf F}\{\cos 2\bar{\theta}\} = \int_{{\cal... ...P}\!\acute{P} \cos\theta_{sr}\cos 2\bar{\theta}\, d{\cal V}\;.$$ (92)

This linearized map is now in standard linear form and can be solved with an analogous Gauss-Newton gradient method as for the case of the $\grave{P}\!\acute{P}$ inverse problem. We derive the new gradient operator ${\bf g}_2$ as:

 

$${\bf g}_2({\bf x};{\bf x}_h) = - {\bf F}'{\bf d}= - \int_{\... ...\grave{P}\!\acute{P}\cos\theta_{sr}D \, d{\bf x}_m\, d\omega\;.$$ (93)

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

 

$$\H_2 = {\bf F}'{\bf F}+ \epsilon^2 = \int_{\omega} \int_{{\b... ...t] \, d{\bf x}_m\, d\omega + \epsilon^2({\bf x};{\bf x}_h) \;,$$ (94)

where

 

$$I_1 = \rho'\omega^2 A'_s A'_r e^{i\omega\tau_{sr}'}\grave{P}\!\acute{P}\cos\theta_{sr}' \{\cdot\} \;.$$ (95)

At this point, $\cos 2\bar{\theta}$ could be estimated as ${\bf m}\approx -\H_2^{-1}{\bf g}_2$ using the gradient and Hessian expressions (93) and (94), and substituting the results of a previous estimation of $\grave{P}\!\acute{P}$ from (86). However, this would be very inefficient for two reasons: (1) the full Hessian of equation (94) requires $N_{{\bf x}}$ forward modeling and migration steps, and (2) the values of $\grave{P}\!\acute{P}$ required by expressions (93) and (94) would need to be calculated by a prior prestack depth migration given by (86). For this reason, we now show how to approximate $\grave{P}\!\acute{P}$ and diagonalize $\H_2$ in (93) and (94), so that a complete estimate of $\grave{P}\!\acute{P}(\bar{\theta})$can be made for the same cost as one single prestack depth migration.

First we diagonalize the Hessian operator (94), by evaluating the ${\cal V}'$ integral in (94) at the stationary point where:

 

$$\int_{V'} \rho'\omega^2 A'_s A'_r e^{i\omega\tau_{sr}'}\grav... ...rho\omega^2 A_s A_r \grave{P}\!\acute{P}\cos 2\bar{\theta} \;,$$ (96)

Now we try to ``undo'' the small angle approximation (90) by noting that at the point of stationary phase, which gives the most contribution to ${\bf g}_2$ and $\H_2$, that

 

$$\cos\theta_{sr}\cos 2\bar{\theta}\rightarrow 1 \;,$$ (97)

as expressed in (91).

With these simplifying assumptions (96) and (97), the diagonal Hessian becomes:

 

$$\H_2 = {\bf F}'{\bf F}+ \epsilon^2 \approx {\bf h}_2 \;,$$ (98)

where the diagonal elements are given as

 

$$\mbox{Diag}\{\H_2\} = \mbox{Diag}\{{\bf h}_2\} = \int_{\omeg... ...sr}} \, d{\bf x}_m\, d\omega+ \epsilon^2({\bf x};{\bf x}_h) \;.$$ (99)

Also, we want to avoid computing a first prestack depth migration for $\grave{P}\!\acute{P}$ to serve as input for a second prestack depth migration estimate of $\cos 2\bar{\theta}$, although one could certainly do this if computer time was not an issue. To economize, we assume a stationary specular reflection coefficient value for $\grave{P}\!\acute{P}$, given from the data as:

 

$$\grave{P}\!\acute{P}\rightarrow \frac{D}{\rho\omega^2 A_s A_r} \;,$$ (100)

since at the stationary point, (87) reduces to:

 

$$D' \rightarrow \rho'\omega^2 A'_s A'_r \grave{P}\!\acute{P}' \;.$$ (101)

Substituting (100) into (93) and (99), we obtain a very efficient estimation solution for the specular reflection angle $\bar{\theta}$ directly from the recorded data D:

 

$$\cos 2\bar{\theta}({\bf x};{\bf x}_h) = -\H_2^{-1}{\bf g}_... ...r}}\, d{\bf x}_m\,d\omega+ \epsilon^2({\bf x};{\bf x}_h) } \;.$$ (102)

Equation (102) is the (implicit) least-squares estimate of $\bar{\theta}$ which we seek.

We now discuss the physical interpretation of (102) to lend insight. Consider a fixed point ${\bf x}$ in the subsurface. As we migrate a constant offset section into ${\bf x}$, the angle $\theta_{sr}$ between the source and receiver rays ranges from $\theta_{sr}\approx 0$ at ${\bf x}_m= {\bf x}_{min}$, to $\theta_{sr} \approx 2\bar{\theta}$ near the specular midpoint, and back to $\theta_{sr}\approx 0$ at ${\bf x}_m={\bf x}_{max}$. Analogously, the diffraction-reflection coefficient $R_{_{P\!P}}=\grave{P}\!\acute{P}\cos\phi_r$ varies from $R_{_{P\!P}}\approx 0$ (pure diffraction) to $R_{_{P\!P}}\approx \grave{P}\!\acute{P}(\bar{\theta})$ (specular reflection), to $R_{_{P\!P}}\approx 0$ (pure diffraction) again, over the same midpoint integration range. Hence, it is apparent that $R_{_{P\!P}}$ will attain a maximal peak amplitude at the specular midpoint, whereupon $R_{_{P\!P}}= \grave{P}\!\acute{P}(\bar{\theta})$ and $\theta_{sr}= 2\bar{\theta}$. Hence, (102) represents a weighted first moment average of the data (squared), which has a central amplitude peak at the midpoint and angle of specular reflection.

It should be noted that the estimate (102) is similar to the result of Bleistein (1987). However, our result contains different WKBJ migration weights, and is a combination of a gradient migration and a normalizing Hessian migration image. Furthermore, our solution is basically two weighted stacks of squared data values, which adds a major robustness advantage, especially at subsurface points ${\bf x}$ where $\vert\grave{P}\!\acute{P}\vert$ tends to be small or zero, where Bleistein's result would tend toward zero division.

As an aside, given that (102) was derived from least-squares optimization theory and involves two sums of squared data values, by analogy we can conjecture that an l₁ estimate of $\bar{\theta}$ would be approximated by:

 

$$l_1: \;\; \cos 2\bar{\theta}({\bf x};{\bf x}_h) \rightarrow ... ...r}}\, d{\bf x}_m\,d\omega+ \epsilon^2({\bf x};{\bf x}_h) } \;.$$ (103)

It is likely that the l₁ angle estimate (103) is even more robust and less sensitive to data outliers than the l₂ result of (102). The relative merits of the l₁ versus the l₂ norm with applications to some seismic problems is discussed by Claerbout and Muir (1973).