paper.tex

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

where $\mathbf{x} = \{x_1, x_2, \ldots\}$ represents the vector of space
coordinates, and the coefficients $a_{ij}$ form a positive-definite
matrix $A$.  Equation (\ref{eqn:elips}) defines the characteristic
surfaces $t = \tau (\mathbf{x})$ for a linear hyperbolic second-order
differential equation of the form
\begin{equation}\label{eqn:wave}
  \sum_{i,j} a_{ij} (\mathbf{x})
  \frac{\partial^2 u}{\partial x_i \partial x_j} +
  F (\mathbf{x}, u, \frac{\partial u}{\partial x_i}) =
  \frac{\partial^2 u}{\partial t^2}\;,
\end{equation}
where F is an arbitrary function.
\par
A known theorem \cite[]{smirnov} states that the propagation rays
[characteristics of equation (\ref{eqn:elips}) and, correspondingly,
bi-characteristics of equation (\ref{eqn:wave})] are geodesic
(extreme-length) curves in the Riemannian metric
\begin{equation}\label{eqn:riman}
  d \tau = \sqrt{\sum_{i,j} b_{ij} (\mathbf{x})\, dx_i\, dx_j}\;,
\end{equation}
where $b_{ij}$ are the components of the matrix $B = A^{-1}$. This
means that a ray path between two points $\mathbf{x}_1$ and $\mathbf{x}_2$
has to correspond to the extreme value of the curvilinear integral
\[
\int_{\mathbf{x}_1}^{\mathbf{x}_2}\,\sqrt{\sum_{i,j} b_{ij} (\mathbf{x})\,
  dx_i\, dx_j}\;.
\]
For the isotropic eikonal equation (\ref{eqn:eikonal}), $a_{ij} =
\delta_{ij}/s^2 (\mathbf{x})$, and metric (\ref{eqn:riman}) reduces to
the familiar traveltime measure
\begin{equation}\label{eqn:simple}
  d \tau = s (\mathbf{x})\, d\sigma\;,
\end{equation}
where $d\sigma = \sqrt{\sum_{i} dx_i^2}$ is the usual Euclidean
distance metric.  In this case, the geodesic curves are exactly
Fermat's extreme-time rays.
\par
From equation (\ref{eqn:riman}), we see that Fermat's principle in the
general variational formulation applies to a much wider class of
situations if we interpret it with the help of non-Euclidean
geometries.

\section{Variational principles on a grid}

In this section, I derive a discrete traveltime computation procedure,
based solely on Fermat's principle, and show that on a Cartesian
rectangular grid it is precisely equivalent to the update formula
(\ref{eqn:update}) of the first-order eikonal solver.

\sideplot{triangle}{width=2in}{A geometrical scheme for the
  traveltime updating procedure in two dimensions.}
\par
For simplicity, let us focus on the two-dimensional case\footnote{A
  very similar analysis applies in three dimensions, but requires a
  slightly more tedious algebra. It is left as an exercise for the
  reader.}.  Consider a line segment with the end points $A$ and $B$,
as shown in Figure \ref{fig:triangle}. Let $t_A$ and $t_B$ denote the
traveltimes from a fixed distant source to points $A$ and $B$,
respectively. Define a parameter $\xi$ such that $\xi=0$ at $A$,
$\xi=1$ at $B$, and $\xi$ changes continuously on the line segment
between $A$ and $B$. Then for each point of the segment, we can
approximate the traveltime by the linear interpolation formula
\begin{equation}\label{eqn:linear}
  t (\xi) = (1-\xi) t_A + \xi t_B\;.
\end{equation}
Now let us consider an arbitrary point $C$ in the vicinity of $AB$. If
we know that the ray from the source to $C$ passes through some point
$\xi$ of the segment $AB$, then the total traveltime at $C$ is
approximately
\begin{equation}\label{eqn:tc}
  t_C = t (\xi) + s_C\,\sqrt{|AB|^2 (\xi-\xi_0)^2 +
  \rho_0^2}\;,
\end{equation}
where $s_C$ is the local slowness, $\xi_0$ corresponds to the
projection of $C$ to the line $AB$ (normalized by the length $|AB|$),
and $\rho_0$ is the length of the normal from $C$ to $\xi_0$.
\par
Fermat's principle states that the actual ray to $C$ corresponds to a
local minimum of the traveltime with respect to raypath perturbations.
According to our parameterization, it is sufficient to find a local
extreme of $t_C$ with respect to the parameter $\xi$.  Equating the
$\xi$ derivative to zero, we arrive at the equation
\begin{equation}\label{eqn:fermat}
 t_B - t_A + \frac{s_C\,|AB|^2\,(\xi-\xi_0)}
 {\sqrt{|AB|^2 (\xi-\xi_0)^2 + \rho_0^2}} = 0\;,
\end{equation}
which has (as a quadratic equation) the explicit solution for $\xi$:
\begin{equation}
  \label{eqn:xi}
  \xi = \xi_0 \pm \frac{\rho_0\,(t_A - t_B)}
  {|AB|\,\sqrt{s_C^2\,|AB|^2 - (t_A - t_B)^2}}\;.
\end{equation}
Finally, substituting the value of $\xi$ from (\ref{eqn:xi}) into
equation (\ref{eqn:tc}) and selecting the appropriate branch of the
square root, we obtain the formula
\begin{eqnarray}\label{eqn:formula}
  c\,t_C & = & \rho_0\,\sqrt{s_C^2 c^2 - (t_A - t_B)^2} +
  c\,t_A\,(1-\xi_0) + c\,t_B\,\xi_0 = \nonumber \\
  & & \rho_0\,\sqrt{s_C^2 c^2 - (t_A - t_B)^2} +
  a\,t_A\,\cos{\beta} + b\,t_B\,\cos{\alpha}\;,
\end{eqnarray}
where $c = |AB|$, $a = |BC|$, $b = |AC|$, angle $\alpha$
corresponds to $\widehat{BAC}$, and angle $\beta$ corresponds to
$\widehat{ABC}$ in the triangle $ABC$ (Figure \ref{fig:triangle}).

\sideplot{square}{width=1.5in}{NR}{A geometrical scheme for
  traveltime updating on a rectangular grid.}
\par
To see the connection of formula (\ref{eqn:formula}) with the eikonal
difference equation (\ref{eqn:update}), we need to consider the case
of a rectangular computation cell with the edge $AB$ being a diagonal
segment, as illustrated in Figure \ref{fig:square}. In this case,
$\cos{\alpha} = \frac{a}{c}$, $\cos{\beta} = \frac{b}{c}$, $\rho_0 =
\frac{ab}{c}$, and formula (\ref{eqn:formula}) reduces to 
\begin{equation}\label{eqn:square}
  t_C = \frac{ab\,\sqrt{s_C^2 (a^2 + b^2) - (t_A - t_B)^2} + a^2\,t_A + b^2\,t_B}
{a^2+b^2}\;.
\end{equation}
We can notice that (\ref{eqn:square}) is precisely equivalent to the
solution of the quadratic equation (\ref{eqn:formula}), which in our
new notation takes the form
\begin{equation}\label{eqn:supdate}
  \left(\frac{t_C - t_A}{b}\right)^2 +
  \left(\frac{t_C - t_B}{a}\right)^2 = s_C^2\;.
\end{equation}
\par
What have we accomplished by this analysis? First, we have derived a
local traveltime computation formula for an arbitrary grid.  The
derivation is based solely on Fermat's principle and a local linear
interpolation, which provides the first-order accuracy.  Combined with
the fast marching evaluation order, which is also based on Fermat's
principle, this procedure defines a complete algorithm of
first-arrival traveltime calculation. On a rectangular grid, this
algorithm is exactly equivalent to the fast marching method of
\cite{paper2} and \cite{mihai}.  Second, the derivation
provides a general principle, which can be applied to derive analogous
algorithms for other eikonal-type (Hamilton-Jacobi) equations and
their corresponding variational principles.

%A very similar analysis applies in 3-D (Appendix A).

\section{Solving the eikonal equation on a triangulated grid}

Unstructured (triangulated) grids have computational advantages over
rectangular ones in three common situations:
 \begin{itemize}
 \item When the number of grid points can be substantially reduced by
   putting them on an irregular grid. This situation corresponds to
   irregular distribution of details in the propagation medium.
 \item When the computational domain has irregular boundaries. One
   possible kind of boundary corresponds to geological interfaces and
   seismic reflector surfaces \cite[]{SEG-1993-0170}. Another type of irregular
   boundary, in application to traveltime computations, is that of
   seismic rays. The method of bounding the numerical eikonal solution
   by ray envelopes has been introduced recently by \cite{SEG-1996-1208}.
 \item When the grid itself needs to be dynamically updated to maintain
   a certain level of accuracy in the computation.
 \end{itemize}
With its computational speed and unconditional stability, the fast
marching method provides considerable savings in comparison with 
alternative, more accurate methods, such as semi-analytical ray
tracing \cite[]{SEG-1991-1497,SEG-1995-1247} or the general Hamilton-Jacobi
solver of \cite{abgrall}.

%\plot{test}{width=6in}{Traveltime contours, computed in the rough
%  Marmousi model (left), the smoothed Marmousi (middle), and the
%  smoothed triangulated Marmousi (right).}
%\par
%Figure \ref{fig:test} shows a comparison between first-arrival
%traveltime computations in regularly gridded and triangulated Marmousi
%models. The two results match each other within the first-order
%accuracy of the fast marching method. However, the cost of the
%triangulated computation has been greatly reduced by constraining the
%number of nodes.
\par
Computational aspects of triangular grid generation are outlined in
Appendix A. A three-dimensional application would follow the same
algorithmic patterns.

\section{Conclusions}

Variational principles have played an exceptionally important role in
the foundations of mathematical physics. Their potential in numerical
algorithms should not be underestimated.
\par
In this paper, I interpret the fast marching eikonal solver with the
help of Fermat's principle. Two important generalizations follow
immediately from that interpretation. First, it allows us to obtain a
fast method of first-arrival traveltime computation on triangulated
grids.  Furthermore, we can obtain a general principle, which extends
the fast marching algorithm to other Hamilton-type equations and their
variational principles. More research is required to confirm these
promises. 
\par
In addition, future research should focus on 3-D implementations and
on increasing the approximation order of the method.

\section{Acknowledgments}
I thank Mihai Popovici and Biondo Biondi for drawing my attention to
the fast marching level set method. Discussions with Bill Symes, Dan
Kosloff, and Tariq Alkhalifah were crucial for developing a general
understanding of the method. Jamie Sethian kindly responded to more
specific questions. The conforming triangulation program was developed
at the suggestion of Leonidas Guibas as a research project for his
Geometrical Algorithms class at Stanford.

\bibliographystyle{segnat}
\bibliography{SEP2,SEG,paper}

%\APPENDIX{A}
%\section{Fermat's principle on a three-dimensional grid}
%\begin{equation}\label{eqn:box}
%  t_D = \frac{\sqrt{3\,\left( s_C^2 (\triangle x)^2 - t_A^2 - t_B^2 - t_C^2\right) +
%      (t_A + t_B + t_C)^2} + (t_A + t_B + t_C)}{3}\;.
%\end{equation}
%\begin{equation}\label{eqn:bupdate}
%  (t_D - t_A)^2 + (t_D - t_B)^2 + (t_D - t_C)^2 = s_D^2 (\triangle x)^2\;.
%\end{equation}
%and so on.

\newpage

%\APPENDIX{A}
\input{delaunay}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: