paper.tex

国外免费地震资料处理软件包

\lefthead{Fomel}
\righthead{Fast marching}
\footer{SEP--95}

\title{A variational formulation \\ of the fast marching eikonal solver}

\email{[e-mail]}

\author{Sergey Fomel}

\maketitle

\begin{abstract}
  I exploit the theoretical link between the eikonal equation and
  Fermat's principle to derive a variational interpretation of the
  recently developed method for fast traveltime computations. This
  method, known as fast marching, possesses remarkable computational
  properties. Based originally on the eikonal equation, it can be
  derived equally well from Fermat's principle.  The new variational
  formulation has two important applications: First, the method can be
  extended naturally for traveltime computation on unstructured
  (triangulated) grids. Second, it can be generalized to handle other
  Hamilton-type equations through their correspondence with
  variational principles.
\end{abstract}

\footnotesize
 \begin{quote}
   Now we are in the rarefied atmosphere of theories of excessive
   beauty and we are nearing a high plateau on which geometry, optics,
   mechanics, and wave mechanics meet on a common ground. Only
   concentrated thinking, and a considerable amount of re-creation,
   will reveal the full beauty of our subject in which the last word
   has not been spoken yet.--Cornelius Lanczos, \emph{The variational
     principles of mechanics}
 \end{quote}
\normalsize

\section{Introduction}

Traveltime computation is one of the most important tasks in seismic
processing (Kirchhoff depth migration and related methods) and modeling.  The
traveltime field of a fixed source in a heterogeneous medium is governed by
the eikonal equation, derived about 150 years ago by Sir William Rowan
Hamilton. A direct numerical solution of the eikonal equation has become a
popular method of computing traveltimes on regular grids, commonly used in
seismic imaging \cite[]{GEO55-05-05210526,GEO56-06-08120821,podvin.gji.91}. A
recent contribution to this field is the \emph{fast marching} level set
method, developed by \cite{paper2} in the general context of level set methods
for propagating interfaces \cite[]{osher,book}. \cite{mihai} report a
successful application of this method in three-dimensional seismic
computations. The fast marching method belongs to the family of upwind
finite-difference schemes aimed at providing the \emph{viscosity} solution
\cite[]{lions}, which corresponds to the first-arrival branch of the
traveltime field. The remarkable stability of the method results from a
specifically chosen order of finite-difference evaluation.  The order
selection scheme resembles the \emph{expanding wavefronts} method of
\cite{GEO57-03-04780487}.  The fast speed of the method is provided by the
heap sorting algorithm, commonly used in Dijkstra's shortest path computation
methods \cite[]{mit}. A similar idea has been used previously in a slightly
different context, in the \emph{wavefront tracking} algorithm of
\cite{GEO59-04-06320643}.
\par
In this paper, I address the question of evaluating the fast marching
method's applicability to more general situations. I describe a simple
interpretation of the algorithm in terms of variational principles
(namely, Fermat's principle in the case of eikonal solvers.) Such an
interpretation immediately yields a useful extension of the method for
unstructured grids: triangulations in two dimensions and tetrahedron
tesselations in three dimensions.  It also provides a constructive way
of applying similar algorithms to solving other eikonal-like
equations: anisotropic eikonal \cite[]{SEG-1991-1530}, ``focusing''
eikonal \cite[]{Biondi.sep.95.biondo1}, kinematic offset continuation
\cite[]{Fomel.sep.84.179}, and kinematic velocity continuation
\cite[]{Fomel.sep.92.159}. Additionally, the variational formulation can
give us hints about higher-order enhancements to the original
first-order scheme.

\section{A brief description of the fast marching method}

For a detailed description of level set methods, the reader is
referred to Sethian's recently published book \cite[]{book}. More
details on the fast marching method appear in articles by
\cite{paper2} and \cite{mihai}. This section serves as a brief
introduction to the main bulk of the algorithm.
\par
The key feature of the algorithm is a carefully selected order of
traveltime evaluation. At each step of the algorithm, every grid point
is marked as either \emph{Alive} (already computed), \emph{NarrowBand}
(at the wavefront, pending evaluation), or \emph{FarAway} (not touched
yet). Initially, the source points are marked as \emph{Alive}, and the
traveltime at these points is set to zero. A continuous band of points
around the source are marked as \emph{NarrowBand}, and their
traveltime values are computed analytically. All other points in the
grid are marked as \emph{FarAway} and have an ``infinitely large''
traveltime value.
\par
An elementary step of the algorithm consists of the following
moves:
 \begin{enumerate}
 \item Among all the \emph{NarrowBand} points, extract the point with
 the minimum traveltime.
 \item Mark this point as \emph{Alive}.
 \item Check all the immediate neighbors of the minimum point and
 update them if necessary.
 \item Repeat.
 \end{enumerate}
\par
An update procedure is based on an upwind first-order approximation to
the eikonal equation. In simple terms, the procedure starts with
selecting one or more (up to three) neighboring points around the
updated point. The traveltime values at the selected neighboring
points need to be smaller than the current value. After the selection,
one solves the quadratic equation
\begin{equation}
\label{eqn:update}
\sum_{j} \left(\frac{t_i-t_j}{\triangle x_{ij}}\right)^2 = s_i^2
\end{equation}
for $t_i$. Here $t_i$ is the updated value, $t_j$ are traveltime
values at the neighboring points, $s_i$ is the slowness at the point
$i$, and $\triangle x_{ij}$ is the grid size in the $ij$ direction.
As the result of the updating, either a \emph{FarAway} point is marked
as \emph{NarrowBand} or a \emph{NarrowBand} point gets assigned a new
value.
\par
Except for the updating scheme (\ref{eqn:update}), the fast marching
algorithm bears a very close resemblance to the famous shortest path
algorithm of \cite{dijkstra}.  It is important to point out that
unlike Moser's method, which uses Dijkstra's algorithm directly
\cite[]{GEO56-01-00590067}, the fast marching approach does not
construct the ray paths from predefined pieces, but dynamically
updates traveltimes according to the first-order difference operator
(\ref{eqn:update}). As a result, the computational error of this
method goes to zero with the decrease in the grid size in a linear
fashion.  The proof of validity of the method (omitted here) is also
analogous to that of Dijkstra's algorithm \cite[]{paper2,book}. As in
most of the shortest-path implementations, the computational cost of
extracting the minimum point at each step of the algorithm is greatly
reduced [from $O (N)$ to $O (\log N)$ operations] by maintaining a
priority-queue structure (heap) for the \emph{NarrowBand} points
\cite[]{mit}.

\inputdir{fmarch}

Figure \ref{fig:salt} shows an example application of the fast marching
eikonal solver on the three-dimensional SEG/EAGE salt model.  The computation
is stable despite the large velocity contrasts in the model. The current
implementation takes about 10 seconds for computing a 100x100x100 grid on one
node of SGI Origin 200. \cite{Alkhalifah.sep.95.tariq5} discuss the
differences between Cartesian and polar coordinate implementations.

\sideplot{salt}{width=3.5in}{Constant-traveltime con\-tours
  of the first-arrival travel\-time, computed in the SEG/EAGE salt
  model. A point source is positioned inside the salt body. The top
  plot is a diagonal slice; the bottom plot, a depth slice.}
\par
The difference equation (\ref{eqn:update}) is a finite-difference approximation
to the continuous eikonal equation
\begin{equation}
\label{eqn:eikonal}
\left(\frac{\partial t}{\partial x}\right)^2 +
\left(\frac{\partial t}{\partial y}\right)^2 +
\left(\frac{\partial t}{\partial z}\right)^2 = s^2 (x,y,z)\;,
\end{equation}
where $x$, $y$, and $z$ represent the spatial Cartesian
coordinates. In the next two sections, I show how the updating
procedure can be derived without referring to the eikonal
equation, but with the direct use of Fermat's principle.

\section{The theoretic grounds of variational principles}

\inputdir{XFig}

This section serves as a brief reminder of the well-known theoretical
connection between Fermat's principle and the eikonal equation.  The
reader, familiar with this theory, can skip safely to the next
section.

\sideplot{fermat}{width=2in}{Illustration of the connection
  between Fermat's principle and the eikonal equation. The shortest
  distance between a wavefront and a neighboring point $M$ is along
  the wavefront normal.}
\par
Both Fermat's principle and the eikonal equation can serve as the
foundation of traveltime calculations. In fact, either one can be
rigorously derived from the other. A simplified derivation of this
fact is illustrated in Figure \ref{fig:fermat}. Following the notation
of this figure, let us consider a point $M$ in the immediate
neighborhood of a wavefront $t (N) = t_N$.  Assuming that the source
is on the other side of the wavefront, we can express the traveltime
at the point $M$ as the sum
\begin{equation}\label{eqn:tm}
  t_M = t_N + l (N,M) s_M\;,
\end{equation}
where $N$ is a point on the front, $l (N,M)$ is the length of the
ray segment between $N$ and $M$, and $s_M$ is the local slowness.
As follows directly from equation (\ref{eqn:tm}),
\begin{equation}\label{eqn:link}
  \left|\nabla t\right| \cos{\theta} = \frac{\partial t}{\partial l}
  = \lim_{M \rightarrow N} \frac{t_M - t_N}{l (N,M)} = s_N\;.
\end{equation}
Here $\theta$ denotes the angle between the traveltime gradient (normal
to the wavefront surface) and the line from $N$ to $M$, and
$\frac{\partial t}{\partial l}$ is the directional traveltime derivative along
that line.
\par
If we accept the local Fermat's principle, which says that the ray
from the source to $M$ corresponds to the minimum-arrival time, then,
as we can see geometrically from Figure \ref{fig:fermat}, the angle
$\theta$ in formula (\ref{eqn:link}) should be set to zero to achieve
the minimum. This conclusion leads directly to the eikonal equation
(\ref{eqn:eikonal}). On the other hand, if we start from the eikonal
equation, then it also follows that $\theta=0$, which corresponds to
the minimum traveltime and constitutes the local Fermat's principle.
The idea of that simplified proof is taken from \cite{lanc},
though it has obviously appeared in many other publications. The
situations in which the wavefront surface has a discontinuous normal
(given raise to multiple-arrival traveltimes) require a more elaborate
argument, but the above proof does work for first-arrival traveltimes
and the corresponding viscosity solutions of the eikonal equation
\cite[]{lions}.
\par
The connection between variational principles and first-order
partial-differential equations has a very general meaning, explained
by the classic Hamilton-Jacobi theory. One generalization of the
eikonal equation is 
\begin{equation}\label{eqn:elips}
  \sum_{i,j} a_{ij} (\mathbf{x})\,
  \frac{\partial \tau}{\partial x_i}\,
  \frac{\partial \tau}{\partial x_j} = 1\;,
\end{equation}