\documentclass {article}
\usepackage[dvips]{graphicx}
\begin{document}
\title{AUTOMATON - A FUZZY AUTOMATIC PICK SYSTEM}
\author{Paul Gettings}
\maketitle
%\tableofcontents
\section*{ABSTRACT}
An automatic first-arrival picking algorithm based on a fuzzy membership
function and empirical wavelet fitting is described.  Synthetic data
can be automatically picked without error by the algorithm.  Field data,
using an elastic wave generator ("thumper"), is more difficult, requiring
multiple manual picks for the shot.  Clean field data can be processed
with at most half, and more likely less than one quarter, the human work
of manual picking.  Computation time is insignificant on a modern desktop,
allowing possible interactive use in the future.

\section*{INTRODUCTION}
\begin{quotation}
Picking of arrivals is a tricky business.  To build a program that can
replace a human is beyond the reach of any mortal programmer.
No number-crunching hardware can approach the pattern-matching
skill of a biological system - to try with modern equipment is to
invite failure.  Therefore, strive not to replace the human picker,
but rather, to augment the human so that picking may be done faster,
better, cheaper.  Strive to build a program that can free the human
from the easy picks; the mundane portion of the data that occupies
so much time and space.  Leave only the hard, complicated, traces
for the human.  There are many fewer of these, and so the human is
happy, and the results are good.

\flushright --Desiderata of an Automatic Picker
\end{quotation}

Travel-time tomography can be a powerful tool for determining
subsurface structure.  However, tomography depends on accurate
determination of the first arrivals at each geophone location for each
shot.  In a large 2d or 3d experiment, there can be tens or hundreds of
thousands of first arrivals that must be picked before tomographic
imaging can be done.  Picking of these arrivals is often more time
consuming than acquiring the data, and hence represents a significant
barrier to wide-spread use of travel-time tomography.  An algorithm
which can automatically determine first arrivals allows for more rapid
data processing, and hence tomography projects of much greater scope
are feasible.

This paper describes a simple automatic picking algorithm based on a
fuzzy membership function and empirical wavelets.  Due to the
assumptions of the algorithm, it has only been tested on controlled
source data, and will likely not be of any benefit in passive-source
situations.

This fuzzy algorithm has been implemented for testing, although the
implementation is certainly not as fully developed as it could be;
a binding for Matlab could be very useful, for example.  Note that
the package was developed on a Unix-like system (GNU/Linux) and hence
subscribes to the Unix philosophy of programs: do one thing, and do
it well.  There is no monolithic interface as part of the this package;
it is meant to be a collection of utilities that are combined in scripts,
custom GUIs, etc.  It is assumed that the code will be running on a
unix-like OS; porting to a Windows-based OS may prove tricky.

\section*{ALGORITHM DESCRIPTION}
The algorithm presented here is quite simple, and is founded on an
examination of how a trained seismologist picks a first arrival: while
looking at a trace, the seismologist finds the first arrival of energy
that is not noise.  The criterion are simple - does the waveform look
"reasonable", is the pick at about the expected time, and does the
waveform have about the "right" amplitude for the geophone's offset.  A
trained seismologist can also interpolate on traces with poor data, and
can use sophisticated knowledge of previous experiments to guide the
choice of waveform, timing, and amplitude.

This algorithm seeks to replicate a seismologist, but does not expect
to replace the seismologist.  Completely automated picking is a very
difficult problem, which is well beyond the scope of this small project.
Instead, the algorithm is designed to pick the simple traces, and leave
the difficult traces for the seismologist.  Hence, part of the algorithm
is a flexible, but reliable, criterion for what is a "good" pick.

Briefly, the algorithm for the picking of the channels of a single shot
is:
\begin{itemize}
\item Compute a normalized empirical wavelet from some known first
arrival, provided by the seismologist
\item Compute membership function, $M$, for each trace:
  \begin{itemize}
  \item[-] Maximize the correlation coefficient from the wavelet to a trace:
    \begin{itemize}
    \item[-] Fit wavelet to each time offset in a data window along trace;
	     keep highest coefficient.
    \end{itemize}
  \item[-] Compare time offset and amplitude of best fit with predicted values
  \end{itemize}
\item If $M$ $>$ threshold, accept as pick;
  else leave the trace for the seismologist
\end{itemize}

Individual steps are described in detail below.

\subsection*{NORMALIZED EMPIRICAL WAVELET}
To compute the wavelet, the seismologist gives the algorithm a trace
with a known first arrival, and the length of the desired wavelet.
Data from the trace is taken and scaled to the range $[-1,1]$.  This
produces a wavelet that can be fit to other traces for varying start
times.

The length of this wavelet is very important.  Longer wavelets are more
unique, and hence the correlation coefficient is a more useful number.
But very long wavelets include too much local information from the
known first arrival, and are very sensitive to changes in waveform.  In
the limited testing done on real field data, a wavelet of at least 20
ms was required for reliably results, and wavelets of over 50 ms caused
many traces to drop below the membership function threshold.

\subsection*{WAVELET FITTING}
Armed with a scaled wavelet, the algorithm can then begin fitting the
wavelet to each trace, for all possible time offsets.  However, to
reduce the number of computations, the data is windowed through a
minimum and maximum expected velocity.  It is assumed that the wavelet
cannot arrive faster than the geophone-source distance divided by the
maximum expected velocity.  Similarly, the wavelet is assumed to arrive
no later than the geophone-source distance divided by the minimum
expected velocity.  For suitably chosen velocities, this procedure can
greatly reduce the number of linear fits required, and hence
dramatically decrease overall computation time.

For each possible time offset in the data window, a linear
least-squares fit is performed for the scaled wavelet to the trace data
using the equation
\begin{displaymath}
T=A \cdot W
\end{displaymath}
where $T$ is the trace data, $A$ is a scalar amplitude, and $W$ is the
scaled wavelet.  This is a simple linear regression problem for a line,
with 0 intercept.  In general, fitting the wavelet to the trace data
should include a phase-shift term.  Since the fitting will be attempted
for each time offset, the phase-shift has been fixed at 0, roughly
halving the number of computations for each fit.

As results of the linear fit, both an amplitude and correlation
coefficient are obtained.  The amplitude will be compared to a
predicted value in the membership function, and the correlation
coefficient is used both in the membership function and in deciding on
the "best" pick for the trace.  The time offset with the largest
coefficient (actually the largest absolute value) is chosen as the
"best" pick for the trace.  This corresponds to choosing the time
offset where the trace data most closely matches the known first arrival
waveform.

Because the wavelet is defined by the known first arrival, this linear
fit will perform poorly on traces with significantly different first
arrival waveforms.  There are two possible ways to handle this problem:
1) break a shot into multiple parts and choose a known first arrival for
each part (this is the current scheme); and 2) use the previous good
pick to compute a new wavelet.  The second method will automatically
handle slowly varying first arrival waveforms, but does require sorting
the data set by distance, and processing in the correct order.  It has
not been implemented or tested, but could greatly improve the overall
efficiency for difficult data.

\subsection*{COMPARISON OF FIT AND PREDICTED TIME AND AMPLITUDE}
After finding the "best" pick for a trace, the time offset and amplitude
of the pick are compared to predicted values.  Amplitude is predicted from
a $r^{-2}$ power law for intensity; amplitude of the known first arrival
is scaled to the shot coordinates to derive a "zero-offset" amplitude
for the source wavelet.  This initial amplitude is then reduced by $r^2$
for each trace, forming a predicted amplitude.  This has the additional
benefit that amplitude units are irrelevant, so long as they are constant.

Time offsets for the first arrival are predicted from a best-guess
velocity for a direct wave, which is supplied by the seismologist.
The start time of the known first arrival is used to compute an initial
time offset at the source position, which is increased by $r/v$, where
$r$ is the source-geophone distance for the trace under analysis, and $v$
is the best-guess velocity.  Only one velocity is used, so large offset
arrivals may result in large time differences as refracted waves
overtake the direct wave.  In this case, either down-weighting of the
time difference, or processing of the shot data in several pieces can
fix the problem.

\subsection*{MEMBERSHIP FUNCTION}
The membership function is a superposition of three functions on the
interval $[0,1]$:
\begin{displaymath}
M = \frac{W_1 |r| + W_2 \exp(-\delta a/C) + W_3 \exp(-\delta t/D)}{\sum
W_i}
\end{displaymath}
where $r$ is the correlation coefficient of the best fit, $\delta a$ is
the difference between the fit and predicted amplitudes, $\delta t$ is the
difference between the fit and predicted time offsets, $W$ is a vector
of weights, $C$ is the amplitude decay constant, and $D$ is the time
decay constant.

From this formula, it is clear that the decay constants scale the delta
values.  For the amplitude function, it is necessary to rescale the
amplitude constant for each trace, so the weighting is done on a relative
measure of misfit, not on the actual value.  Hence, the amplitude constant
is actually:
\begin{displaymath}
  C = A \cdot p
\end{displaymath}
where $p$ is the predicted amplitude, and $A$ is the seismologist-supplied
amplitude decay constant.  By rescaling $C$ for each trace, it is
possible to use the same membership function for all fits, regardless
of the scale of the amplitudes.  Since amplitudes can range over 3-4
orders of magnitude within a shot, this is of critical importance.

The time decay constant, $D$, should be set to approximately twice the
allowed timing error.  This means that a $\delta t$ value within the allowed
timing error will at most evaluate to $\exp(-1/2)$, which will keep
the $\delta t$ contribution significant.

If decay constants are set too small, the membership function values
will be entirely dependent on the correlation coefficient, regardless
of weighting.  If the constants are too large, large deviations in
amplitude and time will be ignored.  In either case, the range of the
membership function will be compressed, which makes distinguishing
"good" and "bad" picks very difficult.

When set correctly, the range of the membership function is maximal,
and hence the spread between "good" and "bad" picks is large, and
easily seen.

\section*{PERFORMANCE TESTING}
\subsection*{ALGORITHM PERFORMANCE}
Performance of the algorithm is measured as a reduction of workload
on the seismologist, compared to hand-picking all channels of all
shots.

A total success is considered a 10-fold reduction in work.  For a shot
with 120 channels, this means the seismologist picks no more than 12
traces.  This is currently thought to be difficult to obtain in practice,
due to several factors: 1) waveforms are clipped near the source, which
makes it basically impossible to automatically pick first arrivals for
$\sim$5 traces next to the source (depending on geophone offsets); 2)
large offset traces often contain significantly different first arrival
waveforms; and 3) large offset traces often have refracted wave first
arrivals, which are significantly earlier than the predicted direct wave.
Between having to subset a shot into several parts for time offset and
waveform concerns, and having to manually pick $\sim$10 traces around
a source point, it is difficult to achieve the 10-fold reduction in the
field data tested so far (120 channels @ 0.5m spacing per shot).

From the limited field data tested so far, it should be trivial to
halve the number of picks, and 4- to 8-fold reduction in picks should
be easy for an experienced seismologist.  Results are sensitive to the
choice of parameters fed to the algorithm, in particular the best-guess
velocity, waveform length, and membership function constants and
weights.  Optimal settings of these parameters varies between
experiments, and possibly between shots within an experiment.  General
guidelines are still being devised using test data sets.

\subsection*{IMPLEMENTATION PERFORMANCE}
The important measure for the implementation is speed of computation.
For the implementation described in this paper, the code is written in
C++ and compiled using GCC 2.95.4 under Debian GNU/Linux.  Compilation
and time tests were done on a dual Athlon MP 1200 MHz with 1 GB DDR RAM
running at 133 MHz.  Due to the large amount of RAM, no memory swapping
was required, and data were read from an Ultra160 SCSI disk.  Data
files varied in size from 14 MB to 27.6 MB.

Running the autopick algorithm on a single shot of 168 channels (14 MB
data file) took a total of 7.48 seconds on the test machine.  Much of
this time was for disk I/O, both the reading of the data file and the
writing of resulting picks and debugging output.  This total time
reduces to 22.46 channels picked per second.  This is almost
unacceptably slow for large problems.  However, running the algorithm
on a data set of 28 shots of 120 channels each (3360 channels total)
took a total of 13.27 seconds including file I/O.  This computes to
253.20 channels per second.

The dramatic rise in channels/second from one test to the other indicates
that much of the total time is spent in disk I/O, not computation.
Hence, it is likely that until data sets reach the tens or hundreds of
thousands of channels, the total run time will be less than 1 minute on
the test machine.  This is fast enough for interactive refinement of
processing parameters.

It should be clear that this computational performance is highly
dependent on the choice of expected velocities, as these set the width
of the data window.  Very broad choices in the velocities will drop the
channels/second rate, as each channel requires more linear wavelet
fits.

Computational performance is also heavily dependent on the compiler
used to generate the machine code.  Test runs were done with an
executable produced from GCC, which is known to produce relatively slow
code to make it more portable.  Anecdotal evidence on the use of the
Intel C++ compiler points to a significant (possibly over twice)
improvement in speed vs. GCC-compiled code.  This has not yet been
tested.

\subsection*{MEMORY USAGE}
This implementation represents both time and amplitude of a seismic
trace explicitly.  SEGY files implicitly express time in the file -
amplitude is stored, and the sample number is multiplied by a time step
to find time.  This means that SEGY files take double the space in RAM
than on disk.  For a modern machine and reasonable SEGY files, this is
not a problem.  A 48 shot survey with 120 channels per shot, 2000
samples per channel, is $\sim$47MB on disk as a SEGY file.  Once loaded, the
autopick routine consumes $\sim$91MB of RAM.

Because of this ballooning of data, it is advisable to run the autopick
algorithm on machines with lots of RAM, and keep the number of
shots and channels small in each file.  Since each shot is processed
individually, there is no need to keep all shots in one SEGY file.

\subsection*{INTEGRATION WITH OTHER PROGRAMS}
For the reference implementation presented here, a functional approach
to the algorithm is used.  Data, picks, and parameters are loaded by the
main function, and then passed to the actual autopick function to do
the work.  This function creates the wavelet, fits the data, and writes
good and bad picks to disk files.

For interfacing with Matlab, etc., revise the autopick function to
take and return data in whatever format is desired.  A wrapper function
around a more general function is probably the best method of devising
the interface.  Wrapping the general autopick function to interface to
Matlab, Java, etc. can then be handled more easily, since the wrapper
function is then only responsible for decoding the incoming data to the
format for the autopick function, and encoding the returning picks for
the calling language.

\section*{ALGORITHM TESTING}

\subsection*{SYNTHETIC DATA}
Initial testing of the algorithm and implementation was done with a
simple synthetic data set.  A simple two-layer model was used, with a
single source at a receiver position.  The computed waveforms were fed
to the autopick algorithm, with a user-supplied known pick based on the
synthetic data.  As expected, the algorithm picks the first arrivals
without error.

This is to be expected for two reasons: 1) waveforms do not alter
across the data as the model has only two homogeneous, isotropic layers;
and 2) trace data is zero for all times except the first arrival, so
there is only one possible wavelet fit with a reasonably large
correlation coefficient.  Neither condition is expected to hold in real
data, but the synthetic data confirm that the algorithm and
implementation are working as designed.

\subsection*{FIELD DATA}
Due to time constraints, only one field data set has been tested so far:
a single shot with 168 channels in a single line at 0.5 m spacing.
Trace data from the shot are shown in figure \ref{trace_fig}.
Data were collected at a 0.5 ms interval on all channels with a
synchronized clock.  Range limitations of the geophones or seismograph
results in clipping of the 5 traces on either side of the source, which
is located at a geophone.
\begin{figure}
  \centerline{
  \includegraphics{traces.eps}
  }
\caption{Plot of trace data for field data test of algorithm. Overlay
numbers denote chunks of data used in testing.}
\label{trace_fig}
\end{figure}

Trace data are reasonably clean, with the exception of the clipping at
near-source geophones.  Because of the clipping, the known first
arrival must be chosen well away from the source; for the algorithm test
trace 64 was used.  Clipping also requires the seismologist to pick the
traces near the source, automatically precluding a 10-fold reduction in
work.  Two examples of the trace data are shown in figure
\ref{trace2_fig}.  The first trace is near the source point, and shows
significant waveform clipping.  The second trace is 30 m from the
source, and hence shows a nice clean waveform.  Waveform clipping
drastically complicates automatic picking, generally to the point where
it is easier to manually pick the first arrivals rather than tweak
the algorithm.
\begin{figure}
  \centerline{
  \includegraphics{trace.eps}
  }
\caption{Plot of trace data for selected traces; (a) is near source
position, note the clipping of the waveform. (b) is $\sim$30 m from
source, note the clean waveform and low noise prior to the first
arrival.}
\label{trace2_fig}
\end{figure}

Even with 0.5 m geophone spacing, offsets become large enough that
the refracted wave overtakes the direct wave inside the survey area.
This means the data should be broken into at least 4 parts for automatic
picking (see figure \ref{trace_fig}): 1) traces from 0 offset to
near-source traces with clipping; 2) traces with clipping; 3) traces
without clipping, with direct wave first arrivals; and 4) traces with
refracted wave first arrivals.  Parts 1, 3, and 4 all require a different
known first arrival for fitting.  Part 2 must be picked entirely manually,
since the clipping destroys too much of the waveform for reliable
automatic picking.  With suitable tweaking of parameters, this allows
for 168 channels to be picked for roughly 20 human picks.

Figure \ref{dt_fig} plots the difference between a manual pick for
each channel and the results of the autopick algorithm.  Note that the
algorithm was run only once, with only a single known pick (at channel 64)
for all channels.  The known first arrival was not optimized by looking at
the data, but chosen nearly at random from among large offset traces; the
time was rounded from the manual pick on the trace.  The heavy horizontal
line marks zero time difference, and the dashed lines delineate $\pm$3
ms, which was taken as the acceptable arrival time error for tomography
processing.  The strong linear trend, and increasing number of rejected
picks, with channel number are to be expected for field data: stations
near the source (at channel 10) are clipped and hence the known arrival
waveform doesn't fit; first arrivals appear to become refracted waves
for channels larger than about 70, and hence the autopick is choosing the
direct wave arrival, not the first arrival; velocity estimates for near-
and far-source channels are inaccurate, so the predicted arrival times
are significantly wrong, which reduces the number of accepted picks.
\begin{figure}
  \centerline{
  \includegraphics{dt.eps}
  }
\caption{Plot of (automatic - manual) first arrival times in ms for
test field data; automatic picks were run as a single data set with the
known arrival at channel 64. Heavy line indicates zero time difference,
and dashed lines mark $\pm$3 ms. Note the strong linear trend
of time difference with channel number (offset), and the increasing
number of rejected picks with larger channel number. Channels below 18
are amplitude clipped, and therefore all picks are rejected.}
\label{dt_fig}
\end{figure}

Bringing the linear trend in figure \ref{dt_fig} flat should be relatively
simple through using multiple manual picks for different subsets of the
data.  Figure \ref{pick_fig} shows the first part of the seismic data
with the automatic and manual picks; there are two sets of automatic
picks: one using channel 64 (32 m offset) as the known pick (green
line), and one using channel 144 (72 m offset) as the known pick.  Note
that the automatic and manual picks are within acceptable error for
several channels around each known.  Figures \ref{pick32m_fig} and
\ref{pick72m_fig} are expanded views around the known pick channels.
For both cases, the automatic and manual pick times diverge as the
waveform changes away from the known pick.  To achieve the 4-fold
reduction in work considered minimally acceptable, it is only necessary
to choose 1 of 4 traces as a manual pick, which is clearly achievable
with the current examples.
\begin{figure}
  \centerline{
  \includegraphics{picks_bw.eps}
  }
\caption{Trace data for test field data; amplitudes have been reversed
to show filling at the first arrivals.  Automatic and manual picks are
shown by the stars and lines on the plot.  Black thin line shows
manual picks; dark grey dashed line shows automatic picks using channel
144; light grey dashed line shows automatic picks using channel 64.}
\label{pick_fig}
\end{figure}

\begin{figure}
  \centerline{
  \includegraphics{picks_32m_bw.eps}
  }
\caption{Trace and pick data for channel 64 as the known arrival (32
m). The same color scheme holds as for figure \ref{pick_fig}.}
\label{pick32m_fig}
\end{figure}

\begin{figure}
  \centerline{
  \includegraphics{picks_72m_bw.eps}
  }
\caption{Trace and pick data for channel 144 as known arrival (72 m
offset). Colors are taken from figure \ref{pick_fig}.}
\label{pick72m_fig}
\end{figure}

\section*{CONCLUSIONS \& FUTURE WORK}
Limited testing of the fuzzy algorithm indicates that suitable parameter
and known first arrival choices can greatly reduce the amount of work
required to pick first arrivals for tomography.  Although the desired
10-fold decrease in work is unlikely with real data, 4-, 6-, or 8-fold
decreases in work can be expected.  Since most tomography experiments
result in thousands to hundreds of thousands of picks, these savings could
be quite significant.  Moreover, the speed of this implementation
allows for an iterative, interactive graphical system where parameters are
changed, shot data broken into subsets, and the picking re-run in near
real-time.  This will allow the seismologist to find optimal parameters,
and check the results immediately.

Several possible improvements, beyond the interactive real-time interface,
are obvious for this implementation.  More good picks might be achieved
by using an adapting wavelet, where the empirical wavelet is recomputed
at each trace with a good pick.  Wavelet construction is trivial, so
the additional computation is insignificant.  Whether the data sorting
and trace ordering problems can be efficiently handled for all cases is
not immediately clear.

Another possible improvement is to try to fit a given trace multiple
times; if the "best" pick is rejected by the membership function, then
re-window the data based on $\delta t$ and try again.  This would greatly
increase the computation time for difficult data sets, as each trace could
be run many times.  However, this should allow the algorithm to eventually
pick the appropriate first arrival for any trace.  This modification
would cause major problems in cases where the best-guess velocity is a
poor approximation of the real effective velocity, however.

The first improvement to the existing implementation should be an
interactive interface, perhaps implemented in Matlab.  Further work
on additional membership function parameters (or different functional
forms) or alternative schemes for determining a "best" pick (e.g. neural
networks, spectral analysis) should wait for more testing and refinement
of this algorithm, to find the practical limitations of the method.
\end{document}
