In principle, an l2 solution for the reflection angles can be estimated directly from the reflection data D by substituting the solution (86) into the normal equation . A difficulty arises from the fact that the kernel functions of (70) are nonlinearly related to by the cosine. In this case,
is a nonlinear forward mapping between the model and the data U.
We solve this nonlinear inverse problem for by linearizing (87) with respect to a new model parameter: .Using the standard cosine expansion,
and using the small angle approximation:
we approximate as:
We note that under this approximation, at the specular reflection (stationary point):
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 to the data U:
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 inverse problem. We derive the new gradient operator as:
Neglecting second-derivative terms, the Hessian can be evaluated in the classic Gauss-Newton sense as:
At this point, could be estimated as using the gradient and Hessian expressions (93) and (94), and substituting the results of a previous estimation of from (86). However, this would be very inefficient for two reasons: (1) the full Hessian of equation (94) requires forward modeling and migration steps, and (2) the values of 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 and diagonalize in (93) and (94), so that a complete estimate of can be made for the same cost as one single prestack depth migration.
First we diagonalize the Hessian operator (94), by evaluating the integral in (94) at the stationary point where:
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 and , that
as expressed in (91).
With these simplifying assumptions (96) and (97), the diagonal Hessian becomes:
where the diagonal elements are given as
Also, we want to avoid computing a first prestack depth migration for to serve as input for a second prestack depth migration estimate of , although one could certainly do this if computer time was not an issue. To economize, we assume a stationary specular reflection coefficient value for , given from the data as:
since at the stationary point, (87) reduces to:
Substituting (100) into (93) and (99), we obtain a very efficient estimation solution for the specular reflection angle directly from the recorded data D:
Equation (102) is the (implicit) least-squares estimate of which we seek.
We now discuss the physical interpretation of (102) to lend insight. Consider a fixed point in the subsurface. As we migrate a constant offset section into , the angle between the source and receiver rays ranges from at , to near the specular midpoint, and back to at . Analogously, the diffraction-reflection coefficient varies from (pure diffraction) to (specular reflection), to (pure diffraction) again, over the same midpoint integration range. Hence, it is apparent that will attain a maximal peak amplitude at the specular midpoint, whereupon and . 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 where 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 l1 estimate of would be approximated by:
It is likely that the l1 angle estimate (103) is even more robust and less sensitive to data outliers than the l2 result of (102). The relative merits of the l1 versus the l2 norm with applications to some seismic problems is discussed by Claerbout and Muir (1973).