\lefthead{Velocity}
\righthead{Assignment 4}
\noemailaddress
\nokeywords
\title{Lab 4 - Interval Velocity estimation}
\email{lomask@sep.stanford.edu}
\author{TA: Jesse Lomask}

\begin{center}
{\bf Due Date: Friday, February 27, 2004 - 10:00 a.m. } \\
{\bf Your Name: John George Diefenbaker}
\end{center}

\footer{IEE}

\ABS{
This lab will guide you through most of the steps - and pitfalls - of
interval velocity estimation, perhaps the most important problem in
exploration geophysics.  You will apply the least squares velocity scan
of Lab 2, and then pick RMS velocities using an autopicker and by hand.
You will apply the multiple realizations concept to bound uncertainty
on your picked RMS velocity.  You will solve the Dix equation as a least
squares problem to estimate interval velocity in 1-D and make some 
preliminary steps in the direction of 2-D.
}

\section{INTRODUCTION}

Exploration seismologists propagate waves through the earth's surface,
and record the traveltime.  In oversimpfied terms, the traveltime is the path
integral of the interval slowness (1/velocity) over the raypath between 
source and receiver.  ``Inversion'' of the traveltime for the underlying 
velocity field is perhaps the most important problem in our branch
of geophysics.

Under suitable assumptions, we can directly estimate root-mean-square 
(RMS) velocity from multioffset seismic data.  You have seen the theory
in {\em Basic Earth Imaging} and experimented with it in Lab 2 of this
course.

However, in order to justify our salaries, we geophysicists must go
beyond RMS velocity,  and estimate a field of interval velocities,
which is needed for seismic depth conversion.

For a stratified ($v(z)$) earth, we can easily write the relationship
between RMS velocity ($V$) and interval velocity ($v$):
\begin{equation}
	V^2(\tau) = \frac{1}{\tau}\  \sum_i v^2_i \Delta\tau_i \label{eq:vrms}
\end{equation}
I say ``easily'' because the above relationship is implied by the words
{\em root-mean-square}.  Equation (\ref{eq:vrms}) has an analytic inverse, 
known as the Dix equation, after the geophysicist who developed it
\cite{dix}.
Let us write the Dix equation directly:
\begin{equation}
	v_k^2 = \frac{ \tau_k V_k^2 - \tau_{k-1} V_{k-1}^2 }{ \tau_k} \label{eq:dix}
\end{equation}
For a detailed derivation, please see {\em Basic Earth Imaging}.  In 
practice, the Dix equation often fails to produce a pleasing result, 
particularly in regions of complex stratigraphy.  

In this lab, you will use your expertise in optimization to cast
the Dix equation as a least squares problem, in order to more robustly 
estimate an interval velocity function from a marine CMP 
gather.

Equation (\ref{eq:vrms}) states the nonlinear relationship between RMS
velocity and interval velocity.  \longcite{harlan} linearizes
the relationship and solves the problem using a Gauss-Newton nonlinear
iteration.  However, notice that equation (\ref{eq:vrms}) is a 
linear mapping between {\em squared} RMS velocity and interval velocity.
This ``linear'' problem is solved in the section of the textbook labeled 
{\em Null Space and Interval Velocity}.

\subsection{Estimating RMS velocity}

Equation (\ref{eq:vrms}) treats (squared) RMS velocity as a known 
quantity.  In fact, the source of this ``data'' is the output of another
important processing step: picking maximum energy from a velocity scan
panel.  You worked in Lab 2 to make a better velocity scan; now you will
use it for something.  A ``perfect'' velocity scan is contingent on the fulfillment 
of a long list of requirements: perfect hyperbolic moveout, infinitely long
recorded offsets, horizontal reflectors, vertical rays, and no multiple
reflections, among others.

Even if we were able to get a perfect velocity scan, the RMS velocity
function which we derive from the velocity scan is only perfectly well 
defined if we have a reflection for each $\tau$!  Oftentimes, picking
RMS velocities is like playing a game of ``connect-the-dots''.  The
picked RMS velocity is the known data in equation (\ref{eq:vrms}), so
in effect, we treat the {\em interpolated} data as data.  Obviously,
RMS picks obtained by interpolation or from weak reflectors should
not be given the same weight as RMS picks from a strong reflection. 

\longcite{harlan} does not use a variable residual weighting, although
he discusses the problem in detail. \longcite{Clapp.sep.97.bob1} 
discuss the use of stack power and semblance as residual weights.
In this lab you will first explore the issue of residual weighting a bit
more and then design a weighting scheme of your own.

\subsection{Interval velocity}

In the second portion of this lab, you will apply least squares 
optimization to ``invert'' equation (\ref{eq:vrms}), and thus compute
the interval velocity ($v(t)$) which best fits the RMS velocity.

First review the {\em Null Space and Interval Velocity} section of the
textbook to familiarize yourself with the notation.  Let us rewrite
equation (\ref{eq:vrms}) as a least squares fitting goal:
\begin{equation}
	\bf W ( Cu - d ) \approx 0  \label{eq:vrms-ls1}
\end{equation}
$\bf u$ is the unknown model, a vector of squared interval velocities.
$\bf d$ is the known data, a vector of squared RMS velocities.
$\bf C$ is the simple causal integration operator.
$\bf W$ is a residual weighting function, assumed to be the inverse
variance of the known data.  Equation (\ref{eq:vrms-ls1}) is notoriously
unstable to high frequency variations in RMS velocity, and moreover, it is
underdetermined in the sense that only strong reflections really qualify
as ``data''.  Therefore, we supplement
the system with a regularization term (first derivative) which penalizes 
``wiggliness'':
\begin{eqnarray}
	\bf W ( Cu - d ) \approx 0 \nonumber \\
	\bf \epsilon  Du \approx 0  \label{eq:vrms-ls2}
\end{eqnarray}
\longcite{GEO64-03-06590661} implement a similar approach for the 
inversion of VSP data for interval velocity.  The textbook solves a 
similar preconditioned system.

\subsection{ ``Blocky'' Models -- Adaptive Regularization }

Without seismic reflections, we would be unable to pick RMS 
velocities.  Ironically, the lithologic discontinuities implied by
reflections are the main cause of instability in Dix inversion.
Ideally, we would like to somehow ``turn off'' the regularization
at the location of a reliable reflector for two reasons.  First,
note that the regularization smoothes the interval velocities,
which compromises the quality of the data fit in the process.  Strong 
reflections imply high data quality, so we want to fit these picks
accurately.  Second, if we turn off the regularization penalty at the
reflector, we allow the velocity to develop discontinuities there,
and thus honor the lithologic discontinuities which lead to 
reflections in the first place.

\longcite{GEO64-03-06590661} mention this strategy.  To quote,

{\em ``Steps in velocity can be accounted for by modifying the
matrix $\bf D$ and thus removing the contribution to the penalty
function across a given depth interval.'' }

We supplement the regularization of equation (\ref{eq:vrms-ls2})
with a ``bad reflection'' selector matrix, $\bf{B}$, which is
zero wherever we have a strong reflection, one elsewhere.
\begin{eqnarray}
	\bf W ( Cu - d ) \approx 0 \nonumber \\
	\bf \epsilon BDu \approx 0  \label{eq:vrms-ls3}
\end{eqnarray}

\section{ YOUR ASSIGNMENT}

\begin{enumerate}

	\item Type {\tt make velscan0.view}.  You will see a CMP gather,
		a velocity scan computed by the code used in Lab 2, and the
		envelope of the velocity scan, computed by SEPLib's {\tt Envelope}.

	\plot{velscan0}{width=6.0in,height=4.0in}{ Left to Right: CMP gather,
		velocity scan (iterative least squares), envelope of velocity scan. }

	\item The velocity scan may not look as you expect.  Adjust the parameters
		in {\tt velscan0.P} to improve the results.  Look in program {\tt Veltran.f90}
		for explanations of the parameters.

	\item What causes the high-velocity (small-slowness) features in the 
		velocity scan?

	\boxit{ YOUR ANSWER: 
	}

	\item Why is it common practice to scan slownesses, rather than velocities?
		{\bf Hint:} Read Harlan's \shortcite{harlan} paper.

	\boxit{ YOUR ANSWER:
	}

	\item Type {\tt make vrms0.view}.  In the upper-left-hand panel, we have the
		``semblance'' derived from the velocity scan of Figure \ref{fig:velscan0},
		overlain by the RMS velocity function picked automatically by the program
		{\tt Pick_rms.rs90}.  Also notice the picking ``fairway'', denoted as the
		region between the blue and yellow lines.  {\tt Pick_rms.rs90} has plenty
		of options--review the code to familiarize yourself with them.  
	
		The upper-right-hand panel contains the ``Stack Power'', or the value of 
		the semblance taken at the picked velocity.  The lower panels show the
		input CMP gather before and after NMO with the picked RMS velocity.

	\plot{vrms0}{width=5.0in,height=5.0in}{ Autopicked RMS velocity, stack power,
 		and NMO result.  }

	\item Has the autopicker done its job correctly?  Does the RMS velocity look 
		as you would hope?  Either way, explain your answer clearly.  If you think
		that the picking fairway should be changed, modify the parameters in file
		{\tt pick0.P}.

	\boxit{ YOUR ANSWER:
	}

	\item You may look at the auto-picked RMS velocity function and think to 
		youself, ``Gee, I could have done a lot better.''  Now you have your
		chance.

		Use {\tt Envelope} to estimate a semblance-type
		plot from {\tt velscan0.H}.  Load the result into the {\tt Ricksep} program
		and pick the semblance peaks.  Save the picks into a file named 
		{\tt picked_reflectors0.txt}.  Type {\tt make picked_reflectors0.H} to place
		your picks into a SEP history file.  {\bf Hint:} {\tt Cat} your semblance 
		panel to itself, then load into {\tt Ricksep}.  For picking hints, look in 
		the menus or ask another student.   If you are using GNOME, then 1)log out  2)log in session failsafe mode  3)type {\tt fvwm2} \verb|&|   4)use Ricksep.

	\item Now you must take these picks and interpolate smoothly between them to
		get an RMS velocity for all times.  You are given a dummy main program called 
		{\tt InterpRMS.f90}.  The program does the I/O for you.  Add the following 
		functionality to the program:
		\begin{itemize}
			\item Enable the program to solve the following least squares missing data 
				problem:
				\begin{eqnarray}
					\bf Bv_{rms} &\approx& \bf d \nonumber \\
					\bf \epsilon \nabla v_{rms} &\approx& \bf 0  \label{eq:vrms-interp}
				\end{eqnarray}
				$\bf v_{rms}$ is the output RMS velocity function.  $\bf B$ is a nearest-neighbor
				binning operator.  $\bf d$ is an irregularly-sampled array of RMS picks 
				({\tt picked_reflectors0.H}).  {\bf Hint:} You are provided with a module, 
				{\tt bin1.r90}, to do 1-D nearest-neighbor binning.  Use it.  
			\item Emulate the code in the previous assignment ({\tt Lapfill2.f90, lapfill2.r90}) 
				to enable your program to do 10 multiple realizations.
				Recall the modified fitting 
				goals:
				\begin{eqnarray}
					\bf Bv_{rms} &\approx& \bf d \nonumber \\
					\bf \epsilon \nabla v_{rms} &\approx& \bf \sigma v  \label{eq:vrms-interp-mr}
				\end{eqnarray}
				$\bf v$ is a vector of random numbers, while $\sigma$ is a scaling factor.
			\item Compute the ``minimum variance'' solution, i.e., the solution to 
				equation (\ref{eq:vrms-interp}).  Use this solution to compute the 
				variance of the suite of models:
				\begin{equation}
					\mbox{var}[\bf v_{rms}] = \frac{1}{N} \sum_{i=1}^N \left( v_{rms}^i - v_{rms}^0 \right)^2
				\end{equation}
				The summation loops over realizations.  $\bf v_{rms}^i$ is the $i^{th}$
				realization; $\bf v_{rms}^0$ is the minimum variance solution.
			\item Use the model variance to design a residual weight, which you will use
				to solve equation (\ref{eq:vrms-ls3}).  Write this weight to disk, using
				the file tag ``weight''.  
		\end{itemize}

	\item Type {\tt make multreal0.view} to view the multiple realizations.  Adjust
		the {\tt randscale} parameter in {\tt mr.P} until each of the realizations looks like
		a plausible RMS velocity, i.e., the output of an autopicker.

	\item Type {\tt make vrms-interp0.view} to view the 
		fruits of your labor.  Hopefully, as evidenced by improved flattening after
		NMO, your interpolated velocity function has improved the situation.

	\item In the figure, you see two residual weights ($\bf W$) for equation
		(\ref{eq:vrms-ls3}).  One is the ``Stack Power'', while the other is the weight 
		you computed with {\tt InterpRMS.f90}.	

		What have we hoped to obtain by using multiple realizations?  In theoretical 
		terms, which of these functions seems like the better choice for ($\bf W$)?  
		Enumerate all the parameters that affect the character of the weight that you 
		computed with  {\tt InterpRMS.f90}.  Explain all your answers in detail.

	\boxit{ YOUR ANSWER:
	}

	\plot{vrms-interp0}{width=6.5in,height=7.0in}{ Autopicked RMS slowness versus
		Handpicked+Interpolated RMS slowness.  Weighting functions derived from 
		stack power of autopicked slowness and from variance of multiple realizations
		of handpicked+interpolated slowness. }
	
	\sideplot{multreal0}{width=2.5in,height=4.0in}{ ``Minimum variance'' RMS velocity
		(eq. \ref{eq:vrms-interp}) in solid line.  Two realizations (eq. 
		\ref{eq:vrms-interp-mr}) in dashed lines. RMS picks superimposed as horizontal
		lines. }

	\item Now we move on to the real goal, which is the estimation of interval
		velocities using equation (\ref{eq:vrms-ls3}).  Type {\tt make vint0.view}.
		You see the result of solving equation (\ref{eq:vrms-ls3}) using autopicked
		and handpicked RMS velocities.  Since you have presumably picked reflections,
		your picks are used to set up the $\bf B$ operator of equation 
		(\ref{eq:vrms-ls3}).

	\plot{vint0}{width=4.5in,height=6.0in}{ Interval velocity functions estimated
		using autopicked and handpicked RMS velocities.  ``Blocky'' result, with
		discontinuities allowed at RMS pick locations. }

	\item Why not just use Dix equation (\ref{eq:dix}) directly?

	\boxit{ YOUR ANSWER:
	}

	\item Keep your copy of {\tt Ricksep} open so that you can easily add/remove
		picks easily.  Experiment with different picking strategies and comment on
		what you learn. The placement and number of picks strongly affects the inversion
		results.

	\boxit{ YOUR ANSWER:
	}

	\item Can you think of certain lithologies that are well suited to use of the 
		blocky model?  What about the smooth model?  Is the blocky model appropriate
		in the case of this data?
	\boxit{ YOUR ANSWER:
	}

	\item You've probably noticed that preconditioning has gone unmentioned. 
		Equation (\ref{eq:vrms-ls-prec}) is currently identical to equation
		(\ref{eq:vrms-ls3}).  Modify the LaTeX fitting goals below to reflect 
		the ``usual'' preconditioned fitting goals.  Fully explain any changes
		you make and any new variables that you add.

	\boxit{
		YOUR ANSWER:

		\begin{eqnarray}
			\bf W ( Cu - d ) \approx 0 \nonumber \\
			\bf \epsilon BDu \approx 0  \label{eq:vrms-ls-prec}
		\end{eqnarray}
	}

	\item Read \longcite{harlan} and explain how his strategy is similar to the
		preconditioned approach in the textbook.

	\boxit{ YOUR ANSWER:
	}

	\item Now we'll go beyond 1-D.  Type {\tt make cmp2d.H; Grey < cmp2d.H | Tube}.
		You'll see 40 CMP's from the 2-D prestack dataset used frequently in {\em Basic
		Earth Imaging}.  Now type {\tt make vrms2d.view}.  This will take a while;
		first, we must perform velocity analysis on each CMP, and then use our autopicker 
		(yes, it has a loop over {\tt n2}) to pick the best RMS velocity at each CMP 
		location.

		Once the panel comes up, you'll see the autopicked RMS velocity and a stack
		of the CMP's.

	\plot{vrms2d}{width=4.5in,height=6.0in}{ Autopicked RMS velocity of one CMP 
		from a 2-D prestack dataset, raw RMS velocity, the CMP in question, and the
		stacked data using the raw RMS velocity.  }

	\item  In general, to compute interval velocities, the 1-D RMS velocity function 
		at each CMP location is inverted independently.  From a programming standpoint,
		we need only add a loop over {\tt n2} in the inversion program.  To see
		what happens when we take this simple tack, type {\tt make vcomp2d.view}.

	\sideplot{vcomp2d}{width=4.0in,height=2.5in}{ Interval velocity computed by
		1-D CMP-by-CMP inversion of the RMS velocity in Figure \ref{fig:vrms2d}.  }

	\item The first application of your new interval velocity might be simple 
		vertical depth conversion of the stacked section in Figure \ref{fig:vrms2d}.  
		Is the estimated interval velocity suitable for this purpose?  Explain why
		or why not.

	\boxit{ YOUR ANSWER:
	}

	\item Look at the RMS velocity.  Describe in detail ways in which the data 
		are correlated, i.e., the dependence of a given 1-D velocity function on its
		neighbors and itself.

	\boxit{ YOUR ANSWER:
	}

	\item Now imagine an ``ideal'' interval velocity, the one you {\em wish} your
		inversion code could produce.  Likewise, describe in detail how you think 
		this ideal interval velocity should be correlated.

	\boxit{ YOUR ANSWER:
	}

	\item Taking your answers to the last two questions into account, clearly 
		propose, in words and equations, a method for using least squares to 
		intelligently extend the interval velocity estimation problem to 2-D.
		You are given the same equations as equation (\ref{eq:vrms-ls3}).
		Fully explain any changes you make and any new variables/operators that 
		you add.

	\boxit{
		YOUR ANSWER:
	}

	\item You may have preconditioned your new equation (\ref{eq:vrms-ls-2d}).
		If you did not, but you think it is possible to make the ``usual'' change
		of variables, write the change of variables now.  If your new equations
		make it impossible to apply the usual preconditioning, explain other ways
		in which you might speed convergence.

	\boxit{ YOUR ANSWER:
	}
\end{enumerate}

\section{DONE}
When you are all done, type \texttt{make paper.print} to 
	print your paper on the SEP's 4th floor printer.
Clean up your directory by typing \texttt{make clean}. 

\bibliographystyle{sep}
\bibliography{SEP,paper}

