[cig-commits] r12011 - seismo/3D/automeasure/latex

[alessia at geodynamics.org]

Author: alessia
Date: 2008-05-23 08:54:13 -0700 (Fri, 23 May 2008)
New Revision: 12011

Added:
   seismo/3D/automeasure/latex/manual_introduction.tex
   seismo/3D/automeasure/latex/manual_method.tex
   seismo/3D/automeasure/latex/manual_technical.tex
   seismo/3D/automeasure/latex/manual_tuning.tex
Removed:
   seismo/3D/automeasure/latex/installation.tex
Modified:
   seismo/3D/automeasure/latex/Makefile
   seismo/3D/automeasure/latex/def_base.tex
   seismo/3D/automeasure/latex/flexwin_manual.tex
   seismo/3D/automeasure/latex/flexwin_paper.pdf
   seismo/3D/automeasure/latex/flexwin_paper.tex
Log:
Started writing flexwin manual

Modified: seismo/3D/automeasure/latex/Makefile
===================================================================
--- seismo/3D/automeasure/latex/Makefile	2008-05-23 07:44:17 UTC (rev 12010)
+++ seismo/3D/automeasure/latex/Makefile	2008-05-23 15:54:13 UTC (rev 12011)
@@ -4,16 +4,19 @@
 SRC_FILES= AM-allcitations.bib \
 abstract.tex \
 acknowledgements.tex \
+appendix.tex \
 conclusion.tex \
 def_base.tex \
 discussion.tex \
+figures_manual.tex \
 figures_paper.tex \
 flexwin_manual.tex \
 flexwin_paper.tex \
-installation.tex \
 introduction.tex \
+manual_introduction.tex \
+manual_method.tex \
+manual_technical.tex \
 method.tex \
-appendix.tex \
 results.tex
 
 MANUAL = flexwin_manual.pdf

Modified: seismo/3D/automeasure/latex/def_base.tex
===================================================================
--- seismo/3D/automeasure/latex/def_base.tex	2008-05-23 07:44:17 UTC (rev 12010)
+++ seismo/3D/automeasure/latex/def_base.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -26,4 +26,4 @@
 \def\degW{$\rm^{\circ}W$}
 \def\degcpkm{${\rm ^{\circ}C\; km}^{-1}$}
 
-
+\newcommand{\trange}[2]{\mbox{{#1}--{#2}~s}}

Modified: seismo/3D/automeasure/latex/flexwin_manual.tex
===================================================================
--- seismo/3D/automeasure/latex/flexwin_manual.tex	2008-05-23 07:44:17 UTC (rev 12010)
+++ seismo/3D/automeasure/latex/flexwin_manual.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -14,15 +14,17 @@
 
 \input{def_base}
 \begin{document}
-\title{FLEXWIN: A flexible automated data-window selection algorithm}
+\title{FLEXWIN: An automated time-window selection algorithm for seismic tomography}
 \author{Alessia Maggi}
 \date{}
 \maketitle
 \tableofcontents
 
-\input{introduction}
-\input{installation}
-\input{method}
+\input{manual_introduction}
+\input{manual_technical}
+\input{manual_method}
+\input{manual_tuning}
+\input{figures_manual}
 
 \bibliographystyle{plainnat}
 \bibliography{AM-allcitations}

Modified: seismo/3D/automeasure/latex/flexwin_paper.pdf
===================================================================
(Binary files differ)

Modified: seismo/3D/automeasure/latex/flexwin_paper.tex
===================================================================
--- seismo/3D/automeasure/latex/flexwin_paper.tex	2008-05-23 07:44:17 UTC (rev 12010)
+++ seismo/3D/automeasure/latex/flexwin_paper.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -18,7 +18,7 @@
 
 \input{def_base}
 
-\newcommand{\trange}[2]{\mbox{{#1}--{#2}~s}}
+%\newcommand{\trange}[2]{\mbox{{#1}--{#2}~s}}
 
 %=============================================
 

Deleted: seismo/3D/automeasure/latex/installation.tex
===================================================================
--- seismo/3D/automeasure/latex/installation.tex	2008-05-23 07:44:17 UTC (rev 12010)
+++ seismo/3D/automeasure/latex/installation.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -1,68 +0,0 @@
-\section{FLEXWIN installation\label{sec:install}}
-
-\subsection{System requirements}
-In order to run, FLEXWIN requires:
-\begin{itemize}
-\item UNIX operating system (Linux, Solaris, MacOS \ldots)
-\item [TODO] Memory requirements (measure what is needed!)
-\item [TODO] Other code (sac / saclib / etc?)
-\end{itemize}
-
-Note: The barebones of this code is not dependent on any input format for data.  In order to use a different data format, just write your own i/o routines.
-
-Note: The filtering and enveloping functionality used in this code comes from
-SacLib - this should be distributed with the code or made available separately.
-
-
-\subsection{Obtaining the code}
-
-[TODO] Write this better once structure of code (and packages that will be
-delivered) is finalised.
-
-The code is available as a gzipped tarball on {\tt ftp:\\ftp-eost.u-strasbg.fr/pub/alessia/code/flexwin.tgz} [NOTE: change this to what everyone agrees with...]
-The tarball is unpacked using the follwing command 
-\begin{verbatim}
-tar xvzf flexwin.tgz
-\end{verbatim}
-to obtain a {\tt FLEXWIN} directory containing the Fortran90 source code ({\tt
-.f90} files), a pdf copy of this manual ({\tt flexwin\_manual.pdf}), {\tt
-INSTALL} and {\tt README} text files with a bare bones instruction set to
-install and run the code, and a {\tt test\_data} directory containing a test
-data set to run the code on.
-
-\subsection{Installation}
-
-[TODO] - Write this when you have cleaned up the distribution
-
-FLEXWIN is made up of a library of the subroutines which actually do the work,
-and an example front-end program which strings the relevant subroutines
-together into a full program.  This front-end program is only an example, showing the correct order of operations.  Feel free to make your own front-end.  
-
-The makefiles contain the basic compilation instruction for the basic code
-(using sac i/o).  It is straightforward to modify them.
-
-make -f makefile clean
-make -f makefile executablename
-
-\subsection{Running the Test case}
-
-[TODO]
-
-\subsection{Input / Output formats}
-
-\subsubsection{Input formats}
-Input formats for data.  In standard version this is SAC evenly spaced files.
-Data and synthetics must be synchronized (i.e. have the same reference time)
-and co-sampled (there is a 1/1000 [CHECK VALUE] built in).  If you wish to use a different input format, here is the list of subroutines you must rewrite/substitute:
-\begin{verbatim}
-[TODO] list of subroutines
-\end{verbatim}
-
-\subsubsection{Output formats}
-The code outputs a set of text files per pair of observed and synthetic
-seismograms analysed.  All text files have a number of header lines, prefixed by {\tt \#} with a set of name-value pairs.
-
-[TODO] list files here - and their formats
-
-[TODO] output in sac format?
-

Added: seismo/3D/automeasure/latex/manual_introduction.tex
===================================================================
--- seismo/3D/automeasure/latex/manual_introduction.tex	                        (rev 0)
+++ seismo/3D/automeasure/latex/manual_introduction.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -0,0 +1,57 @@
+\chapter{Introduction}
+
+Seismic tomography - the process of imaging the 3D structure of the Earth using
+seismic recordings - has been transformed by recent advances in methodology.
+Finite-frequency approaches are been used instead of ray-based approaches, and 3D reference models instead of 1D reference models.
+
+It is common practice in tomography to work only with certain subsets of the
+available seismic data. The choices made in selecting these subsets are
+inextricably linked to the assumptions made in the tomographic method.  For example, ray-based
+travel-time tomography deals
+only with high frequency body wave arrivals, while great-circle
+surface wave tomography must satisfy the path-integral approximation,
+and only considers surface waves that present no evidence of multipathing.  
+In both these examples, a large proportion of the information contained within the seismograms goes to waste.
+The emerging 3D-3D tomographic methods take advantage of
+full wavefield simulations and numeric finite-frequency kernels,
+thereby reducing
+the data restrictions required when using approximate forward modelling and simplified descriptions
+of sensitivity.  
+In order to exploit the
+full power of 3D-3D tomographic methods, we require a new data selection
+strategy that does not exclude such complex arrivals.
+This data selection strategy should define measurement time windows that cover as much of a given seismogram as possible, whilst avoiding portions of
+the waveform that are dominated by noise.
+
+Our answer to the time-window selection problem is FLEXWIN (FLEXible WINdowing).
+
+From a signal processing point of view, the simplest way to avoid serious
+contamination by noise is to select and measure strong signals, which in
+seismology correspond to seismic arrivals.  FLEXWIN selects time windows on the synthetic seismogram within which the waveform
+contains a distinct energy arrival, then requires an adequate correspondence
+between observed and synthetic waveforms within these windows.  This selection
+paradigm is general, and can be applied to synthetic seismograms regardless of
+how they have been obtained.  It is clear, however, that a synthetic seismogram
+obtained by 3D propagation through a good 3D Earth model will provide a better
+fit to the observed seismogram over a greater proportion of its length than
+will be the case for a more approximate synthetic seismogram.
+
+In order to the isolate changes in amplitude or frequency content susceptible of being
+associated with distinct energy arrivals, we need to analyse the character of the synthetic waveform itself.  This analysis is similar to that used on observed waveforms
+in automated phase detection algorithms for the routine location of
+earthquakes.  FLEXWIN takes a tool used in this detection process ---
+the long-term / short-term ratio --- and applies it to the definition of
+time windows around distinct seismic phases.  
+
+The choices made in time-window selection for tomography are
+interconnected with all aspects of the tomographic inversion process,
+from the waveform simulation method (direct problem), through the choice of
+measurement method, to the method
+used to obtain sensitivity kernels.  One of the major difficulties in defining
+a general data selection strategy is the great range of possible choices open
+to the tomographer.  FLEXWIN is a configurable data selection process
+that can be adapted to different tomographic scenarios by tuning a handful of
+parameters (see Table~\ref{tb:params}).  Although we have designed our
+algorithm for use in adjoint tomography, its inherent flexibility should make
+it useful in many data-selection applications.
+

Added: seismo/3D/automeasure/latex/manual_method.tex
===================================================================
--- seismo/3D/automeasure/latex/manual_method.tex	                        (rev 0)
+++ seismo/3D/automeasure/latex/manual_method.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -0,0 +1,522 @@
+\chapter{FLEXWIN, the algorithm\label{sec:algorithm}} 
+
+FLEXWIN 
+operates on pairs of
+observed and synthetic single component seismograms.  There is no restriction
+on the type of simulation used to generate the synthetics, though realistic
+Earth models and more complete propagation theories yield waveforms that are more similar to the observed
+seismograms, and thereby allow the definition of measurement windows
+covering more of the available data.  The input seismograms can be measures of
+displacement, velocity or acceleration, indifferently.  There is no requirement
+for horizontal signals to be rotated into radial and transverse directions.
+
+The window selection process has five phases, each of which is discussed individually
+below: {\em phase 0:} pre-processing; {\em phase A:} definition of preliminary
+measurement windows; {\em phase B:} rejection of preliminary windows based on
+the content of the synthetic seismogram alone; {\em phase C:} rejection of
+preliminary windows based on the differences between observed and synthetic
+seismograms; {\em phase D:} resolution of preliminary window overlaps.  The parameters that permit tuning of the
+window selection towards a specific tomographic scenario are all contained in a
+simple parameter file (see Table~\ref{tb:params}).  More complexity and finer
+tuning can be obtained by rendering some of these parameters time dependent, via user defined functions that can depend on the source parameters (e.g. event location or depth).
+
+\begin{table}
+\begin{tabular}{lp{0.8\linewidth}}
+\hline
+\multicolumn{2}{l}{Standard tuning parameters:} \\[5pt]
+$T_{0,1}$     & band-pass filter corner periods \\
+$r_{P,A}$     & signal to noise ratios for whole waveform \\
+$r_0(t)$      & signal to noise ratios single windows \\
+$w_E(t)$      & water level on short-term:long-term ratio \\
+$CC_0(t)$          & acceptance level for normalized cross-correlation\\
+$\Delta\tau_0(t)$  & acceptance level for time lag \\
+$\Delta\ln{A}_0(t)$   & acceptance level for amplitude ratio \\ 
+\hline
+\multicolumn{2}{l}{Fine tuning parameters:} \\ [5pt]
+$c_0$ & for rejection of internal minima \\
+$c_1$ & for rejection of short windows \\
+$c_2$ & for rejection of un-prominent windows \\
+$c_{3a,b}$  & for rejection of multiple distinct arrivals \\
+$c_{4a,b}$  & for curtailing of windows with emergent starts and/or codas \\
+$w_{CC}\quad w_{\rm len}\quad w_{\rm nwin}$ & for selection of best non-overlapping window combination \\
+\hline
+\end{tabular}
+\caption{\label{tb:params}
+Overview of standard tuning parameters, and of fine
+tuning parameters.  Values are defined in a parameter file, and the
+time dependence of those that depend on time is described by user-defined functions.
+} 
+\end{table}
+
+
+
+%----------------------
+
+\pagebreak
+\section{Phase 0 -- Pre-processing \label{sec:phase0}}
+%{\em Parameters used: $T_{0,1}$.}
+The purpose of this phase is to pre-process input seismograms, to reject
+noisy records, and to set up a secondary waveform (the short-term / long-term average ratio) derived from the envelope of the synthetic seismogram.  This STA:LTA waveform will be used later to define preliminary
+measurement windows.
+
+%----------------------
+
+%\subsubsection{Pre-processing}
+
+We apply minimal and identical pre-processing to both observed and synthetic
+seismograms: band-pass filtering with a non-causal Butterworth
+filter, whose
+short and long period corners we denote as $T_0$ and $T_1$ respectively. 
+Values of these corner periods should reflect the information content of the data,
+the quality of the Earth model, and the accuracy of the simulation used to generate the synthetic seismogram.
+All further references to ``seismograms'' in this paper will refer to these filtered waveforms.
+
+%----------------------
+
+%\subsubsection{Seismogram rejection on the basis of noise in observed seismogram}
+
+Our next step is to reject seismograms that are dominated by noise.  This rejection is based on two signal-to-noise criteria that compare the power and amplitude of the signal to those of the background noise (given by the observed waveform before the first $P$-wave arrival).  The power signal-to-noise ratio is defined as
+${\rm SNR}_P = P_{\rm signal}/P_{\rm noise},$ 
+where the time-normalized power in the signal and noise portions of the data are defined respectively by
+\begin{align}
+P_{\rm signal} & = \frac{1}{t_E-t_A} \int_{t_A}^{t_E}d^2(t)dt, \\ 
+P_{\rm noise}  & = \frac{1}{t_A-t_0} \int_{t_0}^{t_A}d^2(t)dt, \label{eq:noise}
+\end{align}
+where $d(t)$ denotes the observed seismogram, $t_0$ is its start time, $t_A$ is
+set to be slightly before the time of the first arrival, and $t_E$ is the end
+of the main signal (a good choice for $t_E$ is the end of the dispersed surface
+wave).  The amplitude signal-to-noise ratio is defined analogously as
+${\rm SNR}_A = A_{\rm signal}/A_{\rm noise}$, 
+where $A_{\rm signal}$ and $A_{\rm noise}$ are the maximum values of $|d(t)|$
+in the signal and noise time-spans respectively.  The limits for these two
+signal-to-noise ratios are given by the parameters $r_P$ and $r_A$ in Table~\ref{tb:params}.  We reject any record for which
+${\rm SNR}_P < r_P$ or ${\rm SNR}_A < r_A$.
+
+%----------------------
+
+%\subsubsection{Construction of STA:LTA timeseries}
+
+Detection of seismic phase arrivals is routinely performed by automated
+earthquake location algorithms.  We have taken a tool used in this
+standard detection process --- the short-term long-term average ratio (STA:LTA)
+--- and adapted it to the task of defining time windows around seismic phases.  Given a synthetic seismogram $s(t)$, we derive its
+STA:LTA timeseries using an iterative algorithm.
+If we denote the Hilbert transform of the synthetic seismogram by
+$\mathcal{H}[s(t)]$, its envelope $e(t)$ is given by:
+\begin{equation}
+e(t) = | s(t) + i \mathcal{H}[s(t)] |.
+\end{equation}
+In order to create the STA:LTA waveform $E(t)$, we discretize the envelope time
+series with timestep $\delta t$, calculate its short term average
+$S(t_i)$ and its long term average $L(t_i)$ as follows, 
+\begin{align}
+S(t_i) & = C_S \; S(t_{i-1}) + e(t_i) \\
+L(t_i) & = C_L \; L(t_{i-1}) + e(t_i) , 
+\end{align}
+and obtain their ratio:
+$E(t_i) = S(t_i)/L(t_i)$. 
+The constants $C_S$ and $C_L$ determine the decay of the relative
+weighting of earlier parts of the signal in the calculation of the current
+average.  This decay is necessarily longer for the long term average than
+for the short term average, implying that $C_S < C_L < 1$.  The choice of these
+constants determines the sensitivity of the STA:LTA timeseries.  
+\citet{BaiKennett2001} used a similar timeseries to
+analyse the character of broad-band waveforms, and allowed the constants
+$C_S$ and $C_L$ to depend on the dominant period of the waveform under
+analysis.  We have followed their lead in setting
+\begin{equation}
+C_S = 10^{- \delta t / T_0} \qquad {\rm and} \qquad C_L = 10^{-\delta t / 12 T_0}, 
+\end{equation}
+where the use of $T_0$, the low-pass corner period of our band-pass filter,
+substitutes that of the dominant period.  
+
+An example of a synthetic seismogram and its corresponding envelope and STA:LTA timeseries $E(t)$ is
+shown in Figure~\ref{fg:stalta}.  Before the first arrivals on the synthetic
+seismogram, the $E(t)$ timeseries warms up and rises to a plateau.  At each
+successive seismic arrival on the synthetic, $E(t)$ rises to a
+local maximum.  We can see from Figure~\ref{fg:stalta} that these local maxima
+correspond both in position and in width to the seismic phases in the
+synthetic, and that the local minima in $E(t)$ correspond to the
+transitions between one phase and the next.  In the following sections we shall
+explain how we use these correspondences to define time windows.
+
+\begin{figure}
+\center \includegraphics[width=6in]{figures/050295B.050-150/ABKT_II_LHZ_seis_nowin.pdf}
+\caption{\label{fg:stalta}
+Synthetic seismogram and its corresponding envelope and STA:LTA timeseries.
+The seismogram was calculated using SPECFEM3D and the
+Earth model S20RTS \citep{RitsemaEtal2004} for the CMT catalog event
+050295B, whose details can be found in Table~\ref{tb:events}.  The
+station, ABKT, is at an epicentral distance of 14100~km and at an azimuth of 44
+degrees from the event.  The top panel shows the vertical component synthetic
+seismogram, filtered between periods of 50 and 150 seconds. The center panel shows its envelope, and the bottom panel
+shows the corresponding STA:LTA waveform.  The dashed line overlaid on
+the STA:LTA waveform is the water level $w_E(t)$.
+} 
+
+\end{figure}
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\clearpage
+\pagebreak
+\section{Phase A -- Preliminary measurement windows \label{sec:phaseA}}
+%{\em Parameters used: $w_E(t)$}.
+
+The correspondence between local maxima in the STA:LTA waveform $E(t)$ and the
+position of the seismic phases in the synthetic seismogram suggests that we
+should center time windows around these local maxima.  The
+correspondence between the local minima in $E(t)$ and the transition between
+successive phases suggests the time windows should start and end at these local
+minima.  In the case of complex phases, there may be several local maxima and
+minima within a short time-span.  In order to correctly window these complex
+phases, we must determine rules for deciding when adjacent local maxima
+should be part of a single window.  From an algorithmic point
+of view, it is simpler to create all possible combinations of adjacent windows
+and subsequently reject the unacceptable ones, than to consider expanding
+small, single-maximum windows into larger ones. 
+
+We start by defining a water level on $E(t)$ via the time dependent parameter
+$w_E(t)$ in Table~\ref{tb:params}.  The water level shown in
+Figure~\ref{fg:stalta} corresponds to $w_E=0.08$ for the duration of the main
+seismic signal.  Typical values for $w_E$ vary between $0.05$ and $0.25$ depending on the seismological scenario and
+the desired sensitivity.  Once set for typical seismograms for a given
+seismological scenario, it is not necessary to change $w_E$ for each
+seismogram.  This is also true of all the other parameters in
+Table~\ref{tb:params}: once the system has been tuned,
+these parameters remain unchanged and are used for all seismic events in the same scenario. The functional forms of the time-dependent parameters are defined by the user, can depend on
+information about the earthquake source and the receiver, and also
+remain unchanged once the system has been tuned (see Appendix~\ref{ap:user_fn}). 
+For the example in Figure~\ref{fg:stalta}, we have required the water level
+$w_E(t)$ to double after the end of the surface wave arrivals (as defined by
+the epicentral distance and a group velocity of $3.2$~\kmps) so as to avoid
+creating time windows after $R1$.  All local maxima that lie above $w_E(t)$
+are used for the creation of candidate time windows.
+
+We take each acceptable local maximum in turn as a seed maximum, and create all
+possible candidate windows that contain it, as illustrated by
+Figure~\ref{fg:win_composite}a.  Each candidate window is defined by three times: its
+start time $t_S$, its end time $t_E$ and the time of its seed maximum $t_M$.
+The start and end times correspond to local minima in $E(t)$.  It is important
+to note that in many of the window rejection algorithms, $t_M$ will be significant.  For $N$ local maxima that lie above $w_E(t)$, the number of preliminary candidate windows defined in this manner is
+\begin{equation}
+N_{\rm win} = \sum_{n=1}^N \left[nN - (n-1)^2\right] \sim O(N^3).
+\end{equation}
+
+\begin{figure}
+\center \includegraphics[width=6in]{figures/fig/window_composite.pdf}
+\caption{\label{fg:win_composite}
+(a)~Window creation process.  The thick black line represents the STA:LTA
+waveform $E(t)$, and the thick horizontal dashed line its water level $w_E(t)$.
+Local maxima are indicated by alternating red and blue dots, windows are
+indicated by two-headed horizontal arrows.  The time of the local maximum used
+as the window seed $t_M$ is denoted by the position of the dot. Only windows for the fourth local maximum are shown.  (b)~Rejection of candidate windows based on the amplitude of the local minima.  The two deep
+local minima indicated by the grey arrows form virtual barriers. All candidate
+windows that cross these barriers are rejected.
+(c)~Rejection of candidate
+windows based on the prominence of the seed maximum.  The local maxima
+indicated by the grey arrows are too low compared to the local minima
+adjacent to them.  All windows that have these local maxima as their seed are
+rejected (black crosses over the window segments below the time series).
+(d)~Shortening of long coda windows.  The grey bar indicates the maximum coda
+duration $c_{4b} T_0$.  Note that after the rejection based on prominence represented in (c) and before shortening of long coda windows represented in (d), the algorithm rejects candidate windows based on the separation of distinct phases, a process that is illustrated in Figure~\ref{fg:separation}.
+}
+\end{figure}
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\clearpage
+\pagebreak
+\section{Phase B -- Rejection based on the synthetic \label{sec:phaseB}}
+%{\em Parameters used: $T_0$, $w_E(t)$, $c_{0-4}$.}
+
+After having created a complete suite of candidate time windows in the manner
+described above, we start the rejection process.  We reject windows based on
+two sets of criteria concerning respectively the shape of the STA:LTA waveform $E(t)$,
+and the similarity of the observed and synthetic waveforms
+$d(t)$ and $s(t)$ within each window.   Here we describe the first set of
+criteria; the second set is described in the following section.
+
+The aim of shape-based window rejection is to retain the set of candidate
+time windows within which the synthetic waveform $s(t)$ contains well-developed seismic phases or groups of phases. The
+four rejection criteria described here are parameterized by the constants
+$c_{0-3}$ in Table~\ref{tb:params}, and are scaled in time by $T_0$ and in
+amplitude by $w_E(t)$.  We apply these criteria sequentially.
+
+Firstly, we reject all windows that contain internal local minima of $E(t)$
+whose amplitude is less than $c_0 w_E(t)$.  We have seen above that local
+minima of $E(t)$ tend to lie on the transitions between seismic phases.  By
+rejecting windows that span deep local minima, we are in fact forcing partition
+of unequivocally distinct seismic phases into separate time windows (see Figure~\ref{fg:win_composite}b).
+Secondly, we reject windows whose length is less than $c_1 T_0$.  By
+rejecting short windows, we are requiring that time windows be long enough to
+contain useful information.
+Thirdly, we reject windows whose seed maximum $E(t_M)$ rises by less than
+$c_2 w_E(t)$ above either of its adjacent minima.  Subdued local maxima of
+this kind represent minor changes in waveform character, and should not be used
+to anchor time windows.  They may, however, be considered as part of a time window with a more prominent maximum (see Figure~\ref{fg:win_composite}c).
+Lastly, we reject windows that contain at least
+one strong phase arrival that is well separated in time from $t_M$.  The
+rejection is performed using the following criterion:
+\begin{equation}
+%h/h_M > f(\frac{\Delta T}{T_0}; c_{3a},c_{3b}),
+h/h_M > f(\Delta T/T_0; c_{3a},c_{3b}),
+\end{equation}
+where $h_M$ is the height of the seed maximum $E(t_M)$ above the deepest
+minimum between itself and another maximum, $h$ is the height of this other
+maximum above the same minimum, and $f$ is a function of the time
+separation $\Delta T$ between the two maxima (see Figure~\ref{fg:separation}).
+The function $f(\Delta T)$ has the following form:
+\begin{equation}
+f(\Delta T) = 
+\begin{cases}
+c_{3a} & \text{ $\Delta T/T_0 \leq c_{3b}$} \\
+c_{3a}\exp{[-(\Delta T/T_0-c_{3b})^2/c_{3b}^2]} & \text{ $\Delta T/T_0 > c_{3b}$.} 
+\end{cases}
+\label{eq:sep}
+\end{equation}
+If we take
+as an example $c_{3a}=1$, this criterion leads to the automatic rejection of
+windows containing a local maximum that is higher than the seed maximum; it also leads to the rejection of windows containing a local maximum that is 
+lower than the seed maximum if it is also sufficiently distant in time from
+$t_M$.  This criterion allows us to distinguish unseparable phase groups from
+distinct seismic phases that arrive close in time. 
+
+The candidate windows that remain after application of these four rejection
+criteria are almost ready to be passed on to the next stage, in which they will
+be evaluated in terms of the similarity between observed and synthetic
+waveforms within the window limits.  Special precautions may have to be taken,
+however, in the case of windows that contain long coda waves: the
+details of codas are often poorly matched by synthetic seismogram calculations,
+as they are essentially caused by multiple scattering processes.   In order to
+avoid rejecting a nicely fitting phase because of a poorly fitting coda or a
+poorly fitting emergent start, we introduce the $c_4$ tuning parameters, which
+permit shortening of windows starting with monotonically increasing $E(t)$
+or ending with monotonically decreasing $E(t)$.  
+These windows are shortened on the left if they start earlier than $c_{4a} T_0$
+before their first local maximum, and on the right if they end later than
+$c_{4b} T_0$ after their last local maximum (see Figure~\ref{fg:win_composite}d).  
+
+Figures~\ref{fg:win_composite} and~\ref{fg:separation} illustrate the shape based
+rejection procedure (Phase B) on a schematic $E(t)$ time series.  Each
+successive criterion reduces the number of acceptable candidate windows.  A
+similar reduction occurs when this procedure is applied to real $E(t)$ time series, as shown
+by the upper portion of Figure~\ref{fg:win_rej_data}.
+
+\begin{figure}
+\center \includegraphics[width=5in]{figures/fig/window_rejection_separation.pdf}
+\caption{\label{fg:separation}
+Rejection of candidate windows based on the separation of distinct phases.
+(a)~Heights of pairs of local maxima above their intervening minimum.
+(b)~The black line represents $f(\Delta T/T_0)$ from
+equation~(\ref{eq:sep}) with $c_{3a}=c_{3b}=1$.  Vertical bars represent
+$h/h_M$ for each pair of maxima.  Their position along the horizontal axis is
+given by the time separation $\Delta T$ between the maxima of each pair.  The
+color of the bar is given by the color of the seed maximum corresponding to $h_M$.  Bars whose height
+exceeds the $f(\Delta T/T_0)$ line represent windows to be rejected.
+(c)~The windows that have been rejected by this criterion are indicated by black
+crosses.
+}
+\end{figure}
+
+\begin{figure}
+\center \includegraphics[width=6in]{figures/fig/window_rejection_global_data.pdf}
+\caption{\label{fg:win_rej_data} 
+Window rejection applied to real data.
+Top panel: observed (black) and synthetic (red) seismograms for the 050295B event
+recorded at ABKT (see Figure~\ref{fg:stalta}).
+Subsequent panels: candidate windows at different stages, separated into Phase B (shape based rejection) and
+Phase C (fit based rejection).  Each candidate window is indicated by a black
+segment.  The number of windows at each stage is shown to the left of the
+panel.
+}
+\end{figure}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\clearpage
+\pagebreak
+\section{Phase C -- Rejection based on seismogram differences \label{sec:phaseC}}
+%{\bf User parameters: $CC_0(t)$, $\Delta\tau_0(t)$, $\Delta\ln{A}_0(t)$}
+
+After having greatly reduced the number of candidate windows by rejection based
+on the shape of the STA:LTA time series $E(t)$, we are now left with a set of
+windows that contain well-developed seismic phases or
+groups of phases on the synthetic seismogram.  
+The next stage is to evaluate the degree of similarity between the observed and
+synthetic seismograms within these windows, and to reject
+those that fail basic fit-based criteria.  A similar kind of rejection is
+performed by most windowing schemes.  
+
+The quantities we use to define well-behavedness of data within a window are
+signal
+to noise ratio ${\rm SNR}_W$, normalised cross-correlation value between
+observed and synthetic seismograms $CC$,
+cross-correlation time lag $\Delta \tau$, and amplitude ratio $\Delta \ln
+A$.  The signal to noise ratio for single windows is defined as an amplitude
+ratio, ${\rm SNR}_W=A_{\rm window}/A_{\rm noise}$, where $A_{\rm window}$ and
+$A_{\rm noise}$ are the maximum absolute values of the observed seismogram $|d(t)|$ in the window
+and in the noise time-span respectively (see equation~\ref{eq:noise}).  The cross-correlation value $CC$ is defined as the maximum value of the
+cross-correlation function ${\rm CC}={\rm max}[\Gamma(t^\prime)]$, where
+\begin{equation}
+\Gamma(t^\prime) = \int s(t-t^\prime)d(t)dt, 
+\end{equation}
+and
+quantifies the similarity in shape between the $s(t)$ and $d(t)$
+waveforms.  The time lag $\Delta \tau$ is defined as the value of $t^\prime$
+at which $\Gamma$ is maximal, and quantifies the delay in time between a
+synthetic and observed phase arrival. The amplitude ratio $\Delta \ln A$ is
+defined as the amplitude ratio between observed and synthetic
+seismograms \citep{DahlenBaig2002}
+\begin{equation}
+\Delta\ln{A} = \left[ \frac{\int d(t)^2 dt}{\int s(t)^2 dt}   \right]^{1/2} - 1. \label{eq:dlnA_def}
+\end{equation}
+The limits that trigger rejection of windows based on the values of these four
+quantities are the time dependent parameters $r_0(t)$, $CC_0(t)$, $\Delta
+\tau_0(t)$ and $\Delta \ln A_0(t)$ in Table~\ref{tb:params}.
+As for the STA:LTA water level $w_E(t)$ used in above, the functional form of
+these parameters is defined by the user, and can depend on source and receiver
+parameters such as epicentral distance and earthquake depth.   
+Figure~\ref{fg:criteria} shows the time
+dependence of $CC_0$ , $\Delta \tau_0$ and $\Delta \ln A_0$ for an example seismogram.  
+
+We only accept candidate windows that satisfy all of the following:
+\begin{align}
+{\rm SNR}_W & \geq r_0(t_M), \label{eq:snr_win} \\
+{\rm CC} & \geq {\rm CC}_0(t_M), \label{eq:cc} \\
+|\Delta\tau| & \leq \Delta\tau_0(t_M), \label{eq:tau} \\
+|\Delta\ln{A}| & \leq \Delta\ln{A}_0(t_M), \label{eq:dlnA}
+\end{align}
+where $t_M$ is the time of the window's seed maximum.  In words, we only accept
+windows in which the observed signal is above the noise level, the observed and
+synthetic signals are reasonably similar in shape, their arrival times
+differences are small, and their amplitudes are broadly compatible.  When the synthetic and observed
+seismograms are similar, the fit-based criteria of
+equations~(\ref{eq:cc})-(\ref{eq:dlnA}) reject only a few of the candidate data
+windows (see lower portion of Figure~\ref{fg:win_rej_data}).  They are
+essential, however, in eliminating problems due secondary events (natural or
+man-made), diffuse noise sources, or instrumental glitches. 
+
+
+\begin{figure}
+\center \includegraphics[width=6in]{figures/050295B.050-150/ABKT_II_LHZ_criteria.pdf}
+\caption{\label{fg:criteria}
+Time dependent fit based criteria 
+for the 050295B event recorded at ABKT. The time-dependence of these criteria
+is given by the formulae in Appendix~\ref{ap:user_global}. The lower limit on
+acceptable cross-correlation value, $CC_0$ (solid line), is
+0.85 for most of the duration of the seismogram; it is lowered to 0.75 during
+the approximate surface wave window  defined by the group velocities 4.2\kmps\
+and 3.2\kmps, and is raised to 0.95 thereafter.  The upper limit on time lag,
+$\tau_0$ (dotted line), is 21~s for the whole seismogram.  The upper limit on amplitude
+ratio, $\Delta \ln A_0$ (dashed line), is 1.0 for most of the seismogram; it is reduced to
+1/3 of this value after the end of the surface waves.  
+}
+\end{figure}
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\clearpage
+\pagebreak
+\section{Phase D -- Overlap resolution \label{sec:phaseD}}
+%{\em User parameters: $w_{CC}$, $w_{\rm len}$.}
+
+After having rejected candidate data windows that fail any of the shape or
+similarity based criteria described above, we are left with a small number of
+windows, each of which taken singly would be an acceptable time window for
+measurement.  As can be seen from Figure~\ref{fg:win_composite}d and the last
+panel of Figure~\ref{fg:win_rej_data}, the remaining windows may
+overlap partially or totally with their neighbours.  Such overlaps are
+problematic for automated measurement schemes, as they lead to multiple
+measurements of those features in the seismogram that lie within the overlapping
+portions.  Resolving this overlap problem is the last step in the
+windowing process.
+
+Overlap resolution can be seen as a set of choices leading to
+the determination of an optimal set of time windows.  What do we mean by
+optimal?  For our purposes, an optimal set of time windows contains only windows that
+have passed all previous tests, that do not overlap with other windows in the set,
+and that cover as much of the seismogram as possible.  When choosing between
+candidate windows, we favour those within which the
+observed and synthetic seismograms are most similar (high values of $CC$).
+Furthermore, should we have the choice between two short windows and a longer,
+equally well-fitting one covering the same time-span, we may wish to favour
+the longer window as this poses a stronger constraint on the tomographic inversion. 
+
+The condition that optimal windows should have passed all previous tests
+removes the straightforward solution of merging overlapping windows.  Indeed, given any two
+overlapping windows, we know that the window defined by their merger 
+existed in the complete list of candidate windows obtained at the end of
+Phase~A, and that its absence from the current list means it was rejected
+either because of the shape of its $E(t)$ time-series (Phase~B), or because of
+an inadequate similarity between observed and synthetic waveforms (Phase~C).
+It would therefore be meaningless to re-instate such a window at this stage.
+Any modification of current candidate windows would be disallowed by similar
+considerations.  We must therefore choose between overlapping
+candidates.  
+
+We make this choice by constructing all possible non-overlapping subsets of
+candidate windows, and scoring each subset on three criteria: length of
+seismogram covered by the windows, average cross-correlation value for the windows,
+and total number of windows.  These criteria often work against each other. For
+example, a long window may have a lower $CC$ than two shorter ones, if the two
+short ones have different time lags $\Delta\tau$.  An optimal weighting of the
+three scores is necessary, and is controlled by the three weighting parameters
+$w_{CC}$, $w_{\rm len}$ and $w_{\rm nwin}$ in Table~\ref{tb:params}. 
+
+As can be seen in Figure~\ref{fg:phaseD}, the generation of subsets is
+facilitated by first grouping candidate windows such that no group overlaps
+with any other group.  The selection of the optimal subsets can then be
+performed independently within each group.  We score each non-overlapping
+subset of windows within a group using the following three metrics:
+\begin{align}
+S_{CC} &= \sum_i^{N_{\rm set}} CC_i / N_{\rm set},\\
+S_{\rm len} &= [\sum_i^{N_{\rm set}} t^e_i - t^s_i]/[t^e_g - t^s_g], \\
+S_{\rm nwin} & = 1 - N_{\rm set}/N_{\rm group},
+\end{align}
+where $CC_i$ is the cross-correlation value of the $i$th window in
+the subset, $N_{\rm set}$ is the number of windows in the subset, $N_{\rm
+group}$ is the number of windows in the group, and $t^s_i$, $t^e_i$, $t^s_g$
+and $t^e_g$ are respectively the start and end times of the $i$th candidate
+window in the set, and of the group itself.  The three scores
+are combined into one using the weighting parameters:
+\begin{equation}
+S = \frac{w_{CC}S_{CC}+w_{\rm len}S_{\rm len}+w_{\rm nwin}S_{\rm nwin}}{w_{CC}+w_{\rm len}+w_{\rm nwin}}.
+\label{eq:score}
+\end{equation}
+The best subset of candidate windows within each group is the one with the
+highest combined score $S$.  The final, optimal set of windows is
+given by concatenating the best subsets of candidate windows for each group.
+Figure~\ref{fg:res_abkt} shows an example of optimal windows selected on real
+data.
+
+\begin{figure}
+\center \includegraphics[width=5in]{figures/fig/window_overlap.pdf}
+\caption{\label{fg:phaseD}
+The selection of the best non-overlapping window
+combinations.  Each grey box represents a distinct group of windows.
+Non-overlapping subsets of windows are shown on separate lines.  Only one
+line from within each group will be chosen, the one corresponding to the
+highest score obtained in equation~(\ref{eq:score}).  The resulting optimal set
+of data windows is shown by thick arrows.}
+\end{figure}
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{figure}
+\center \includegraphics[width=6in]{figures/fig/window_results.pdf}
+\caption{\label{fg:res_abkt}
+Window selection results for event 050295B
+from Table~\ref{tb:events} recorded at ABKT ($37.93$\degN,
+$58.11$\degE, $\Delta=127$\deg, vertical
+component).
+(a)~Top: observed and synthetic seismograms (black and red
+traces); bottom: STA:LTA timeseries $E(t)$.  Windows chosen by the algorithm
+are shown using light blue shading.  The phases contained these windows are:
+(1) $PP$, (2) $PS+SP$, (3) $SS$, (4) $SSS$, (5) $S5$, (6) $S6$, (7) fundamental
+mode Rayleigh wave.
+(b)~Ray paths corresponding to the body wave phases present in the data windows.
+}
+\end{figure}
+%

Added: seismo/3D/automeasure/latex/manual_technical.tex
===================================================================
--- seismo/3D/automeasure/latex/manual_technical.tex	                        (rev 0)
+++ seismo/3D/automeasure/latex/manual_technical.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -0,0 +1,71 @@
+\chapter{Obtaining and installing FLEXWIN\label{sec:install}}
+Here is where you find basic information for obtaining and installing the code.
+For details of the algorithm, see chapter~\ref{sec:algorithm}.  For details of
+how to tune the algorithm to your system, see chapter~\ref{sec:tuning}.
+
+\section{System requirements}
+In order to run, FLEXWIN requires:
+\begin{itemize}
+\item UNIX operating system (Linux, Solaris, MacOS \ldots)
+\item [TODO] Memory requirements (measure what is needed!)
+\item [TODO] Other code (sac / saclib / etc?)
+\end{itemize}
+
+Note: The barebones of this code is not dependent on any input format for data.  In order to use a different data format, just write your own i/o routines.
+
+Note: The filtering and enveloping functionality used in this code comes from
+SacLib - this should be distributed with the code or made available separately.
+
+
+\section{Obtaining the code}
+
+[TODO] Write this better once structure of code (and packages that will be
+delivered) is finalised.
+
+The code is available as a gzipped tarball on {\tt ftp:\\ftp-eost.u-strasbg.fr/pub/alessia/code/flexwin.tgz} [NOTE: change this to what everyone agrees with...]
+The tarball is unpacked using the follwing command 
+\begin{verbatim}
+tar xvzf flexwin.tgz
+\end{verbatim}
+to obtain a {\tt FLEXWIN} directory containing the Fortran90 source code ({\tt
+.f90} files), a pdf copy of this manual ({\tt flexwin\_manual.pdf}), {\tt
+INSTALL} and {\tt README} text files with a bare bones instruction set to
+install and run the code, and a {\tt test\_data} directory containing a test
+data set to run the code on.
+
+\section{Installation}
+
+[TODO] - Write this when you have cleaned up the distribution
+
+FLEXWIN is made up of a library of the subroutines which actually do the work,
+and an example front-end program which strings the relevant subroutines
+together into a full program.  This front-end program is only an example, showing the correct order of operations.  Feel free to make your own front-end.  
+
+The makefiles contain the basic compilation instruction for the basic code
+(using sac i/o).  It is straightforward to modify them.
+
+make -f makefile clean
+make -f makefile executablename
+
+\section{Running the Test case}
+
+[TODO]
+
+\section{Input / Output formats}
+
+\subsection{Input formats}
+Input formats for data.  In standard version this is SAC evenly spaced files.
+Data and synthetics must be synchronized (i.e. have the same reference time)
+and co-sampled (there is a 1/1000 [CHECK VALUE] built in).  If you wish to use a different input format, here is the list of subroutines you must rewrite/substitute:
+\begin{verbatim}
+[TODO] list of subroutines
+\end{verbatim}
+
+\subsubsection{Output formats}
+The code outputs a set of text files per pair of observed and synthetic
+seismograms analysed.  All text files have a number of header lines, prefixed by {\tt \#} with a set of name-value pairs.
+
+[TODO] list files here - and their formats
+
+[TODO] output in sac format?
+

Added: seismo/3D/automeasure/latex/manual_tuning.tex
===================================================================
--- seismo/3D/automeasure/latex/manual_tuning.tex	                        (rev 0)
+++ seismo/3D/automeasure/latex/manual_tuning.tex	2008-05-23 15:54:13 UTC (rev 12011)
@@ -0,0 +1 @@
+\chapter{Tuning FLEXWIN for your problem\label{sec:tuning}}