Application to velocity estimation

The velocity domain representation of seismic data is an alternative to the standard CMP presentation. Transformation of CMP data into the velocity domain (producing a velocity model or panel of the data) exhibits clearly the moveout inherent in the data and therefore, forms a convenient basis for velocity analysis as a linear inverse problem. The velocity transform A from the model space (velocity domain) into the data space (CMP gathers) stretches the velocities back in the offset plane (superposition of hyperbolas) whereas the adjoint operation (A') squeezes the data (summation over hyperbolas):

$$\bold{A}=\bold{HS},$$

with

$$\bold{S}m(t,x) = \displaystyle \left. \sum_{s}\frac{t}{\tau}w(s,x,\tau)m(\tau,s) \right\vert _{\tau=\sqrt{t^2-s^2x^2}}\$$

$$\bold{A'}=\bold{S'H'},$$

with

$$\bold{S'}d(\tau,s) = \displaystyle \left. \sum_{x}w(s,x,\tau)d(t,x) \right\vert _{t=\sqrt{\tau^2+s^2x^2}}\$$

where $w(s,x,\tau)$ is a weighting function, $\bold{H}$ is a filter that we define later. $\bold{A}$ is related to the velocity stack as defined by Taner and Koehler (1969).

The problem is: given a CMP gather can we find a velocity panel which synthesizes it $\it{via}$$\bold{A}$? In equations, given data $\bold{d}$, we want to solve for model $\bold{m}$:

$$\bold{Am} \approx \bold{d}.$$

A simple way to solve this problem is to find a model $\bold{m}$ that minimizes the mean square misfit

$$(\bold{Am}-\bold{d})'(\bold{Am}-\bold{d}).$$

This optimization problem is equivalent to the linear system (``normal equations'')

$$\bold{A'Am} = \bold{A'd}.$$

This system is easy to solve if $\bold{A'A} \approx \bold{I}$, i.e if $\bold{A}$ is close to unitary: then $\bold{m}=\bold{A'd}$.In general, $\bold{A}$ is far from an unitary operator for many reasons. However, the choice of a weighting function compensates to some extent for geometrical spreading and other effects Claerbout and Black (1997):

$$w(s,x,\tau)=\frac{1}{\sqrt{(\tau^2+s^2x^2)^{1/2}}}\frac{\tau}{\sqrt{\tau^2+s^2x^2}}\sqrt{xs}.$$

The summation in the velocity space boosts low frequencies. Claerbout and Black (1997) suggest that a good choice of filter $\bold{H}$ is a half derivative operator ($\sqrt{i\omega}$). These choices for $\bold{H}$ and $w(s,x,\tau)$ bring $\bold{A}$ closer to being an unitary operator.

Since the data is noisy, the modeling operator is not unitary and the numbers of equations and unknowns may be large, an iterative data-fitting approach seems reasonable:

$$\rm{min}_\bold{m}(\it{E}(\bold{Am - d}))$$

where $\bold{m}$ is the model, $\bold{d}$, the data we want to fit, $\bold{A}$ the modeling operator, and E a misfit measurement function we have to choose. We have already presented one possibility, namely that E is the l² norm (least squares inversion). A convenient iterative method in this case is to solve the normal equation using conjugate gradient iteration. We refer to this approach as ``CG'' or ``l²''. An alternative approach is to take for E the Huber function introduced in the first section. With this Huber misfit measure, the velocity transform inverse problem is no longer equivalent to a linear system. We choose to solve it using a general-purpose nonlinear optimizer, as mentioned above, rather than one of several special-purpose methods invented for this type of problem. We refer to this approach as''Huber'' or the ``Huber solver''.

The next two parts of this paper compare the performance of the CG algorithm to Huber in the velocity analysis application.