[cigcommits] r14389  seismo/3D/ADJOINT_TOMO/flexwin_paper/latex
alessia at geodynamics.org
alessia at geodynamics.org
Thu Mar 19 03:11:33 PDT 2009
Author: alessia
Date: 20090319 03:11:30 0700 (Thu, 19 Mar 2009)
New Revision: 14389
Removed:
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/figures_manual.tex
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_method.tex
Modified:
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/AMallcitations.bib
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/Makefile
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/acknowledgements.tex
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_manual.pdf
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_manual.tex
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_paper.pdf
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_introduction.tex
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_other.tex
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_technical.tex
seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_tuning.tex
Log:
Fixed manual
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/AMallcitations.bib
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/AMallcitations.bib 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/AMallcitations.bib 20090319 10:11:30 UTC (rev 14389)
@@ 4182,12 +4182,14 @@
Volume = 250,
Year = 2006}
 at article{MaggiEtal2008,
+ at article{MaggiEtal2009,
Author = {Maggi, A. and Tape, C. and Chen, M. and Chao, D. and Tromp, J.},
Journal = gji,
Title = {An automated timewindow selection algorithm for seismic tomography},
Volume = XX,
 Year = 2008}
+ Year = 2009,
+ note = {(in press)}
+ }
@article{MJJ99,
Author = {Mahatsente, R. and Jentzsch, G. and Jahr, T.},
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/Makefile
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/Makefile 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/Makefile 20090319 10:11:30 UTC (rev 14389)
@@ 8,13 +8,11 @@
conclusion.tex \
def_base.tex \
discussion.tex \
figures_manual.tex \
figures_paper.tex \
flexwin_manual.tex \
flexwin_paper.tex \
introduction.tex \
manual_introduction.tex \
manual_method.tex \
manual_tuning.tex \
manual_technical.tex \
manual_other.tex \
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/acknowledgements.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/acknowledgements.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/acknowledgements.tex 20090319 10:11:30 UTC (rev 14389)
@@ 8,6 +8,6 @@
Additional global scale data were provided by the GEOSCOPE network.
We thank the Hinet Data Center (NIED), especially Takuto Maeda and Kazushige Obara, for their help in providing the seismograms used in the Japan examples.
For the southern California examples, we used seismograms from the Southern California Seismic Network, operated by California Institute of Technology and the U.S.G.S.
The FLEXWIN code makes use of filtering and enveloping algorithms that are part of SAC (Seismic Analysis Code, Lawerence Livermore National Laboratory) provided for free to IRIS members. We thank Brian Savage for adding interfaces to these algorithms in recent SAC distributions.
+The FLEXWIN code makes use of filtering and enveloping algorithms that are part of SAC (Seismic Analysis Code, Lawerence Livermore National Laboratory) provided for free to IRIS members. We thank Brian Savage for adding interfaces to these algorithms in recent SAC distributions. We thank Vala Hjorleifsdottir for her constructive suggestions during the development of the code.
We thank Jeroen Ritsema and an anonymous reviewer for insightful comments that
helped improve the manuscript.
Deleted: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/figures_manual.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/figures_manual.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/figures_manual.tex 20090319 10:11:30 UTC (rev 14389)
@@ 1,300 +0,0 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Table captions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\begin{table}
\begin{tabular}{lrrrrrl}
\hline
Identifier & Latitude & Longitude & Depth, km & Moment, N m & $M_w$ & Location \\
\hline
\multicolumn{7}{c}{Global} \\ \hline
% CHECK THAT THE MOMENT IS LISTED IN NM, NOT DYNECM
% CARL HAS FORMULAS TO CONVERT FROM A MOMENT TENSOR TO M0 TO MW
101895B & 28.06 & 130.18 & 18.5 & 5.68e19 & 7.1 & Ryukyu Islands \\
050295B & 3.77 & 77.07 & 112.8 & 1.27e19 & 6.7 & Northern Peru \\
060994A & 13.82 & 67.25 & 647.1 & 2.63e21 & 8.2 & Northern Bolivia \\
\hline
\multicolumn{7}{c}{Japan} \\ \hline
051502B & 24.66 & 121.66 & 22.4 & 1.91e18 & 6.1 & Taiwan \\
200511211536A & 30.97 & 130.31 & 155.0 & 2.13e18 & 6.2 & Kyuhu, Japan \\
091502B & 44.77 & 130.04 & 589.4 & 4.24e18 & 6.4 & Northeastern China \\
\hline
\multicolumn{7}{c}{Southern California} \\ \hline
9983429 & 35.01 & 119.14 & 13.5 & 9.19e15 & 4.6 & Wheeler Ridge, California \\
9818433 & 33.91 & 117.78 & 9.4 & 3.89e15 & 4.3 & Yorba Linda, California \\
\hline
\end{tabular}
\caption{\label{tb:events}
Example events used in this study. The identifier refers to the CMT catalog for global events and Japan events, and refers to the Southern California Earthquake Data Center catalog for southern California events.
}
\end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\begin{table}
\begin{center}
\begin{tabular}{lccccc}
\hline
 & Global & \multicolumn{2}{c}{Japan} & \multicolumn{2}{c}{S. California} \\
\hline
$T_{0,1}$ & 50, 150 & 24, 120 & 6, 30 & 6, 40 & 2, 40 \\
$r_{P,A}$ & 3.5, 3.0 & 3.5, 3.0 & 3.5, 3.0 & 3.5, 3.0 & 3.5, 2.5 \\
$r_0$ & 2.5 & 1.5 & 3.0 & 2.5 & 4.0 \\
$w_E$ & 0.08 & 0.11 & 0.12 & 0.22 & 0.07 \\
$CC_0$ & 0.85 & 0.70 & 0.70 & 0.74 & 0.85 \\
$\Delta\tau_0$ & 15 & 12.0 & 3.0 & 3.0 & 2.0 \\
$\Delta\ln{A}_0$& 1.0 & 1.0 & 1.0 & 1.5 & 1.0 \\
\hline
$c_0$ & 0.7 & 0.7 & 0.7 & 0.7 & 1.0 \\
$c_1$ & 4.0 & 3.0 & 3.0 & 2.0 & 4.0 \\
$c_2$ & 0.3 & 0.0 & 1.0 & 0.0 & 0.0 \\
$c_{3a,b}$ & 1.0, 2.0 & 1.0, 2.0 & 1.0, 2.0 & 3.0, 2.0 & 4.0, 2.5 \\
$c_{4a,b}$ & 3.0, 10.0 & 3.0, 25.0 & 3.0, 12.0 & 2.5, 12.0 & 2.0, 6.0 \\
$w_{CC}, w_{\rm len}, w_{\rm nwin}$
 & 1, 1, 1 & 1, 1, 1 & 1, 1, 1 & 1, 0, 0 & 1, 0, 0.5 \\
\hline
\end{tabular}
\caption{\label{tb:example_params}
Values of standard and finetuning parameters for the three seismological
scenarios discussed in this study.
}
\end{center}
\end{table}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure captions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\begin{figure}
\center \includegraphics[width=6in]{figures/fig/examples_global.pdf}
\caption{\label{fg:examples}
(a)~Window selection results for event 101895B from Table~\ref{tb:events} recorded
at LBTB ($25.01$\degS, $25.60$\degE, $\Delta=113$\deg, radial component).
Phases contained within selected windows:
(1)~$SKS$, (2)~$PS+SP$, (3)~$SS$, (4)~fundamental mode Rayleigh wave (5) unidentified late phase.
(b)~Body wave ray paths corresponding to data windows in (a).
(c)~Window selection results for event 060994A from Table~\ref{tb:events} recorded at WUS ($41.20$\degN, $79.22$\degE, $\Delta=140$\deg, transverse component).
Phases contained within selected windows: (1)~$S_{\rm diff}$, (2)~$sS_{\rm diff}$, (3)~$SS$, (4)~$sSS$ followed by $SSS$, (5)~$sS5+S6$, (6)~$sS6+S7$ followed by $sS7$, (7)~major arc $sS4$, (8)~major arc $sS6$.
(d)~Body wave ray paths corresponding to data windows in (c).
}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\begin{figure}
\center \includegraphics[width=6in]{figures/fig/composites_global.pdf}
\caption{\label{fg:composites}
(a)(c)~Summary plots of windowing results for event 101895B in Table~\ref{tb:events}.
(a)~Global map showing greatcircle paths to stations.
(b)~Histograms of number of windows as a function of normalised crosscorrelation $CC$, timelag $\tau$ and amplitude ratio $\Delta \ln A$; these give information about systematic trends in time shift and amplitude scaling.
(c)~Record sections of selected windows for the vertical, radial and transverse components. The filled portions of the each record in the section indicate where windows have been selected by the algorithm.
(d)(f)~Summary plots of windowing results for event 060994A in Table~\ref{tb:events}.
}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JAPAN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\clearpage
%\begin{figure}
%%\center
%\includegraphics[width=5.7in]{figures/japan/ERM_II_051502B}
%\caption{\label{fg:ERM_II_051502B}
%Window selection results for event 051502B from Table~\ref{tb:events} recorded at station ERM ($42.01$\degN, $143.16$\degE, $\Delta=24.83$\deg).
%(a)~Event and station map: event is 051502B indicated by the beach ball with the
%CMT focal mechanism, station ERM is marked as red triangles and all the other stations
%which recorded this event are marked by grey triangles.
%(b)~Results for station ERM for the period range \trange{24}{120}.
%Vertical (Z), radial (R), and transverse (T) records of data (black, left column) and synthetics (red, left column), as well as the STA/LTA records (right column) used to produce the window picks.
%(c)~Results for station ERM for the period range \trange{6}{30}.
%}
%\end{figure}
%\clearpage


\begin{figure}
%\center
\includegraphics[width=5.7in]{figures/japan/KIS_BO_091502B}
\caption{\label{fg:KIS_BO_091502B}
Window selection results for event 091502B from Table~\ref{tb:events} recorded at station KIS ($33.87$\degN, $135.89$\degE, $\Delta=11.79$\deg).
(a)~Event and station map: event is 091502B indicated by the beach ball with the
CMT focal mechanism, station KIS is marked as red triangles and all the other stations
which recorded this event are marked by grey triangles.
(b)~Results for station KIS for the period range \trange{24}{120}.
Vertical (Z), radial (R), and transverse (T) records of data (black, left column) and synthetics (red, left column), as well as the STA/LTA records (right column) used to produce the window picks.
(c)~Results for station KIS for the period range \trange{6}{30}.
}
\end{figure}

\clearpage
\begin{figure}
%\center
\includegraphics[width=5.7in]{figures/japan/SHR_BO_200511211536A}
\caption{\label{fg:SHR_BO_200511211536A}
Window selection results for event 20051121536A from Table~\ref{tb:events} recorded at station SHR ($44.06$\degN, $144.99$\degE, $\Delta=17.47$\deg).
(a)~Event and station map: event 20051121536A is indicated by the beach ball with the
CMT focal mechanism, station SHR is marked as red triangles and all the other stations
which recorded this event are marked by grey triangles.
(b)~Results for station SHR for the period range \trange{24}{120}.
Vertical (Z), radial (R), and transverse (T) records of data (black, left column) and synthetics (red, left column), as well as the STA/LTA records (right column) used to produce the window picks.
(c)~Results for station SHR for the period range \trange{6}{30}.
Note that corresponding lowfrequency bandpassed filtered version (b) has longer record length (800~s).
}
\end{figure}

\clearpage
\begin{figure}
%\center
\includegraphics[width=6in]{figures/japan/200511211536A_T06_rs}
\caption{\label{fg:200511211536A_T06_rs}
Summary plots of windowing results for event 200511211536A in Table~\ref{tb:events}, for the period range \trange{6}{30}.
(a)~Map showing paths to each station with at least one measurement window.
(b)(d)~Histograms of number of windows as a function of normalised crosscorrelation $CC$, timelag $\tau$ and amplitude ratio $\Delta \ln A$.
(e)(g)~Record sections of selected windows for the vertical, radial and transverse components.
}
\end{figure}

\clearpage
\begin{figure}
%\center
\includegraphics[width=6in]{figures/japan/200511211536A_T24_rs}
\caption{\label{fg:200511211536A_T24_rs}
Summary plots of windowing results for event 200511211536A in Table~\ref{tb:events}, for the period range \trange{24}{120}.
}
\end{figure}

\clearpage
\begin{figure}
%\center
\includegraphics[width=6in]{figures/japan/stats_T06}
\caption{\label{fg:T06_rs}
Summary statistics of windowing results for events 051502B, 200511211536A and 091502B in Table~\ref{tb:events}, for the period range \trange{6}{30}.
}
\end{figure}


%
%\clearpage
%\begin{figure}
%%\center
%\includegraphics[width=6in]{figures/japan/051502B_T06_rs}
%\caption{\label{fg:051502B_T06_rs}
%Summary plots of windowing results for event 051502B in Table~\ref{tb:events}, for the period range \trange{6}{30}.
%Same as Figure~\ref{fg:200511211536A_T06_rs}.
%}
%\end{figure}
%
%\clearpage
%\begin{figure}
%%\center
%\includegraphics[width=6in]{figures/japan/091502B_T06_rs}
%\caption{\label{fg:091502B_T06_rs}
%Summary plots of windowing results for event 091502B in Table~\ref{tb:events}, for the period range \trange{6}{30}.
%Same as Figure~\ref{fg:200511211536A_T06_rs}.
%}
%\end{figure}

%\clearpage
%\begin{figure}
%%\center
%\includegraphics[width=6in]{figures/japan/091502B_T06_rs}
%\caption{\label{fg:091502B_T06_rs}
%Summary plots of windowing results for event 051502B in Table~\ref{tb:events},
%for the period range \trange{6}{30}. Same as Figure~\ref{fg:200511211536A_T06_rs).
%}
%\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SOUTHERN CALIFORNIA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\begin{figure}
%\center
\includegraphics[width=6in]{figures/socal/9818433_CLC_window.pdf}
\caption{\label{fg:socal_CLC}
Window selection results for event 9818433 from Table~\ref{tb:events} recorded at station CLC.
(a)~Source and station information for event 9818433 and station CLC.
(b)~Paths to each station with at least one measurement window for the period range \trange{6}{40}.
There are a total of 341 windows picked within 310 records.
Triangle denotes station CLC.
(c)~Paths to each station with at least one measurement window for the period range \trange{2}{40}.
There are a total of 190 windows picked within 193 records.
Triangle denotes station CLC.
(d)~Results for station CLC for the period range \trange{6}{40}.
Vertical (Z), radial (R), and transverse (T) records of data (black, left column) and synthetics (red, left column), as well as the STA:LTA records (right column) used to produce the window picks.
(e)~Results for station CLC for the period range \trange{2}{40}.
Note that corresponding lowerpassed filtered versions are shown in (d).
}
\end{figure}

\clearpage
\begin{figure}
%\center
\includegraphics[width=6in]{figures/socal/9818433_FMP_window.pdf}
\caption{\label{fg:socal_FMP}
Window selection results for event 9818433 from Table~\ref{tb:events} recorded at station FMP.
Same caption as Figure~\ref{fg:socal_CLC}, only for a different station.
}
\end{figure}

\clearpage
\begin{figure}
%\center
\includegraphics[width=6in]{figures/socal/9983429_T06_rs.pdf}
\caption{\label{fg:socal_rs_T06}
Summary plots of windowing results for event 9983429 in Table~\ref{tb:events}, for the period range \trange{6}{40}.
(a)~Map showing paths to each station with at least one measurement window.
(b)(d)~Histograms of number of windows as a function of normalised crosscorrelation $CC$, timelag $\tau$ and amplitude ratio $\Delta \ln A$.
(e)(g)~Record sections of selected windows for the vertical, radial and transverse components.
The two branches observed on the vertical and radial components correspond to the bodywave arrivals and the Rayleigh wave arrivals.
}
\end{figure}

%\clearpage
%\begin{figure}
%%\center
%\includegraphics[width=7in]{figures/socal/9983429_T02_rs.pdf}
%\caption{\label{fg:socal_rs_T02}
%(THIS FIGURE COULD IN THEORY BE CUT OUT, IF SPACE IS SHORT.)
%Summary plots of windowing results for event 9983429 in Table~\ref{tb:events}, for the period range \trange{2}{40}.
%Same as Figure~\ref{fg:socal_rs_T06}, only the windowing code has been run using a different set of parameters (Table~\ref{tb:example_params}), so that primarily only the bodywave arrivals are selected.
%}
%\end{figure}


%\clearpage
%\begin{figure}
%%\center
%\includegraphics[width=7in]{figures/socal/9818433_T06_CLC_adj.pdf}
%\caption{\label{fg:socal_adj}
%Adjoint sources constructed based on the windows picked in Figure~\ref{fg:socal_CLC}d, with the specification of a crosscorrelation traveltime measurement. The adjoint sources for this measurement are simply a weighted version of the synthetic velocity traces. The number to the left of each subplot is the $\pm$ height of the $y$axis. The crosscorrelation measurements for traveltime ($\Delta T$) and amplitude ($\Delta \ln A$) are listed above each time window.
%}
%\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_manual.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_manual.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_manual.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_manual.tex 20090319 10:11:30 UTC (rev 14389)
@@ 6,6 +6,7 @@
\usepackage{setspace}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
+\usepackage{url}
\usepackage{natbib}
\pagestyle{headings}
\usepackage[noend]{algorithmic}
@@ 17,7 +18,7 @@
\input{def_base}
\begin{document}
\title{FLEXWIN: An automated timewindow selection algorithm for seismic tomography}
+\title{FLEXWIN User Manual}
\author{Alessia Maggi}
\date{}
\maketitle
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/flexwin_paper.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_introduction.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_introduction.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_introduction.tex 20090319 10:11:30 UTC (rev 14389)
@@ 24,6 +24,6 @@
the algorithm was designed for use in 3D3D adjoint tomography, its inherent
flexibility should make it useful in many dataselection applications.
For a detailed introduction to FLEXWIN as applied to seismic tomography, please consult \cite{MaggiEtal2008} ({\tt flexwin/latex/flexwin\_paper.pdf}). If you use FLEXWIN for your own research, please cite \cite{MaggiEtal2008}.
+For a detailed introduction to FLEXWIN as applied to seismic tomography, please consult \cite{MaggiEtal2009}. If you use FLEXWIN for your own research, please cite \cite{MaggiEtal2009}.
Deleted: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_method.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_method.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_method.tex 20090319 10:11:30 UTC (rev 14389)
@@ 1,521 +0,0 @@
\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:} preprocessing; {\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}$ & bandpass 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 shortterm:longterm ratio \\
$CC_0(t)$ & acceptance level for normalized crosscorrelation\\
$\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 unprominent 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 nonoverlapping 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 userdefined functions.
}
\end{table}


%

\pagebreak
\section{Phase 0  Preprocessing \label{sec:phase0}}
%{\em Parameters used: $T_{0,1}$.}
The purpose of this phase is to preprocess input seismograms, to reject
noisy records, and to set up a secondary waveform (the shortterm / longterm average ratio) derived from the envelope of the synthetic seismogram. This STA:LTA waveform will be used later to define preliminary
measurement windows.

%

%\subsubsection{Preprocessing}

We apply minimal and identical preprocessing to both observed and synthetic
seismograms: bandpass filtering with a noncausal 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 signaltonoise 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 signaltonoise ratio is defined as
${\rm SNR}_P = P_{\rm signal}/P_{\rm noise},$
where the timenormalized power in the signal and noise portions of the data are defined respectively by
\begin{align}
P_{\rm signal} & = \frac{1}{t_Et_A} \int_{t_A}^{t_E}d^2(t)dt, \\
P_{\rm noise} & = \frac{1}{t_At_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 signaltonoise 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 timespans respectively. The limits for these two
signaltonoise 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 shortterm longterm 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_{i1}) + e(t_i) \\
L(t_i) & = C_L \; L(t_{i1}) + 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 broadband 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 lowpass corner period of our bandpass 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.050150/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 timespan. 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, singlemaximum 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 timedependent 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  (n1)^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 twoheaded 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_{04}$.}

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 shapebased window rejection is to retain the set of candidate
time windows within which the synthetic waveform $s(t)$ contains welldeveloped seismic phases or groups of phases. The
four rejection criteria described here are parameterized by the constants
$c_{03}$ 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_0c_{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 welldeveloped 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 fitbased criteria. A similar kind of rejection is
performed by most windowing schemes.

The quantities we use to define wellbehavedness of data within a window are
signal
to noise ratio ${\rm SNR}_W$, normalised crosscorrelation value between
observed and synthetic seismograms $CC$,
crosscorrelation 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 timespan respectively (see equation~\ref{eq:noise}). The crosscorrelation value $CC$ is defined as the maximum value of the
crosscorrelation function ${\rm CC}={\rm max}[\Gamma(t^\prime)]$, where
\begin{equation}
\Gamma(t^\prime) = \int s(tt^\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 fitbased 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
manmade), diffuse noise sources, or instrumental glitches.


\begin{figure}
\center \includegraphics[width=6in]{figures/050295B.050150/ABKT_II_LHZ_criteria.pdf}
\caption{\label{fg:criteria}
Time dependent fit based criteria
for the 050295B event recorded at ABKT. The timedependence of these criteria
is given by the formulae in Appendix~\ref{ap:user_global}. The lower limit on
acceptable crosscorrelation 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 wellfitting one covering the same timespan, 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)$ timeseries (Phase~B), or because of
an inadequate similarity between observed and synthetic waveforms (Phase~C).
It would therefore be meaningless to reinstate 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 nonoverlapping subsets of
candidate windows, and scoring each subset on three criteria: length of
seismogram covered by the windows, average crosscorrelation 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 nonoverlapping
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 crosscorrelation 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 nonoverlapping window
combinations. Each grey box represents a distinct group of windows.
Nonoverlapping 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}
%
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_other.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_other.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_other.tex 20090319 10:11:30 UTC (rev 14389)
@@ 3,12 +3,16 @@
To report bugs or suggest improvements to the code, please send an email to the CIG Computational Seismology Mailing List (cigseismo at geodynamics.org) or Alessia Maggi (alessia at sismo.ustrasbg.fr), and/or use our online bug tracking system Roundup (www.geodynamics.org/roundup).
\section{Notes and Acknowledgments}
[FIXME] The filtering routines used in {\tt seismo\_subs.f90} are based on the SacLib libraries constructed by Brian Savage from the original source code of SAC (developed at Lawrence Livermore). What about SAC licences??
+The main developers of the FLEXWIN source code are Alessia Maggi and Carl Tape. The following individuals (listed in alphabetical order) have also contributed to the development of the source code: Daniel Chao, Min Chen, Vala Hjorleifsdottir, Qinya Liu, Jeroen Tromp. The following individuals (listed in alphabetical order) contributed to this manual: Sue Kientz, Alessia Maggi, Carl Tape.
The main developers of the FLEXWIN source code are Alessia Maggi and Carl Tape. The following individuals (listed in alphabetical order) have also contributed to the development of the source code: Daniel Chao, Min Chen, Jeroen Tromp. The following individuals (listed in alphabetical order) contributed to this manual: Sue Kientz, Alessia Maggi, Carl Tape \ldots
+The FLEXWIN code makes use of filtering and enveloping algorithms that are part of SAC (Seismic Analysis Code, Lawerence Livermore National Laboratory) provided for free to IRIS members. We thank Brian Savage for adding interfaces to these algorithms in recent SAC distributions.
+We acknowledge support by the National Science Foundation under grant EAR0711177.
+Daniel Chao received additional support from a California Institute of Technology Summer Undergraduate Reseach Fellowship.
+
+
\section{Copyright}
Copyright XXX 2008, by the California Institute of Technology (U.S.) and University of Strasbourg (France). ALL RIGHTS RESERVED. U.S. Government Sponsorship Acknowledged.
+Copyright 2009, by the California Institute of Technology (U.S.) and University of Strasbourg (France). ALL RIGHTS RESERVED. U.S. Government Sponsorship Acknowledged.
Any commercial use must be negotiated with the Office of Technology Transfer at
the California Institute of Technology. This software may be subject to U.S.
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_technical.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_technical.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_technical.tex 20090319 10:11:30 UTC (rev 14389)
@@ 12,75 +12,29 @@
\end{itemize}
FLEXWIN requires the following libraries external to the package in order to
compile and run: libsacio.a and libSacLib.a. Eventually, both libraries
will be distributed by IRIS as part of the SAC package (at the moment only
libsacio is distributed this way). For the time being, you should compile
libSacLib.a using the source code in the SacLib directory that accompanies
flexwin.
+compile and run: {\tt libsacio.a} and {\tt libsac.a}. Both libraries
+are distributed by IRIS as part of the SAC package (version 101.2 and above).
+The IRIS download site (as of 30Sept2008) is here:
+\url{http://www.iris.edu/software/sac/sac.request.htm}.
+(To check your version, type sac.)
{\bf Important note}: The SacLib directory is a temporary fix. The SAC source code
from which the SacLib library is compiled is proprietary and should not be
distributed by anyone other than IRIS. Brian Savage  the author of SacLib 
is currently working on a new version of the library that will be distributed
with future versions of SAC. The official release of flexwin will require
linkage to this new library.
\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 from CIG (Computational Infrastructure for Geodynamics, \url{http://www.geodynamics.org}). The tarball is unpacked by typing {\tt tar xvzf flexwin.tgz}.
The code is available as a gzipped tarball from CIG (Computational Infrastructure for Geodynamics, {\tt http://www.geodynamics.org}). The tarball is unpacked by typing {\tt tar xvzf flexwin.tgz}.

The package contains the flexwin code and documentation, as well as a set of
test data, examples of user files for different scenarios, and a set of utility
scripts that may be useful for running flexwin on large datasets.
The contents of the flexwin directory are as follows:
{\small
\begin{verbatim}
flexwin
 Makefile.in
 PAR_FILE
 README
 TODOs
 configure
 configure.ac
 distaz.f
 flexwin.f90
 io_subs.f90
 latex
 make_gfortran
 make_intel
 make_intel_caltech
 maxima.f90
 measure_windows_xcorr.f90
 measurement_module.f90
 scripts
 seismo_subs.f90
 select_windows_stalta2.f90
 test_data
 travel_times.f90
 ttimes_mod
 user_files
 user_functions.f90
 user_parameters.f90
` xcorrmeasure.f90
\end{verbatim}
}

\section{Compilation}
[TODO]  Rewrite this for the official release.

{\bf Note}: Do NOT use the configure script for beta test compilation. It will not
work.

If your compiler of choice is gfortran, then you should be able to use the
{\tt make\_gfortran} makefiles with only minor modifications (notably you may need to
change the search path for the {\tt libsacio.a} library). If you prefer another
compiler, you should modify the OPT and FC lines in the makefiles accordingly.
+compiler, you should modify the OPT and FC lines in the makefiles accordingly. We tested the code using gfortran version 4.1.2
+(To check your version, type{\tt gfortran version}.)
{\bf Important note}: All the code is compiled with the m32 option, which makes
32bit binaries. This option is currently required to enable compatibility with
@@ 89,19 +43,11 @@
Steps to compile the flexwin package:
\begin{enumerate}
\item Compile {\tt libSacLib.a}. In the {\tt SacLib} directory (which is outside the {\tt flexwin}
directory) type: {\tt make f make\_gfortran}.
\item Compile {\tt libtau.a} and create {\tt iasp91.hed} and {\tt iasp91.tbl}. In the
{\tt flexwin/ttimes\_mod} directory type: {\tt make f make\_gfortran}. This will compile
{\tt libtau.a}, and two programs, {\tt remodl} and {\tt setbrn}. The makefile will also run
{\tt remodl} and {\tt setbrn} to create the {\tt iasp91.hed} and {\tt iasp91.tbl} files. You should
then type {\tt make f make\_gfortran install} to install the iasp91 files.
\item Compile {\tt flexwin}. Edit the {\tt make\_gfortran} file in the root directory to ensure the {\tt SACLIBDIR} variable points to the location of your SAC libraries (by default {\tt /opt/sac/lib}). Then type {\tt make f make\_gfortran}.
+\item Compile {\tt libtau.a} and create {\tt iasp91.hed} and {\tt iasp91.tbl}. In the {\tt flexwin/ttimes\_mod} directory type: {\tt make f make\_gfortran}. This will compile {\tt libtau.a}, and two programs, {\tt remodl} and {\tt setbrn}. The makefile will also run {\tt remodl} and {\tt setbrn} to create the {\tt iasp91.hed}and {\tt iasp91.tbl} files. You should then type {\tt make f make\_gfortran install} to install the iasp91 files.
+\item Compile flexwin. Edit the {\tt make\_gfortran} file in the flexwin root directory to ensure the {\tt SACLIBDIR} environment variable points to the location of your SAC libraries (by default {\tt \$SACHOME/lib}). Then type {\tt make f make\_gfortran}.
\end{enumerate}
You should end up with the {\tt flexwin} executable. The program requires the iasp91
files (or links to them) to be present in the directory from which the code is
launched.
+You should end up with the {\tt flexwin} executable. The program requires the {\tt iasp91.hed} and {\tt iasp91.tbl} files (or symbolic links to them) to be present in the directory from which the code is launched.
\section{Running the Test case}
@@ 118,7 +64,7 @@
file. Your result should be identical to that shown in Figure~\ref{fg:test_data}.
\begin{figure}
\center \includegraphics[width=4in]{../test_data/MEASURE.orig/ABKT_II_LHZ_seis.pdf}
+\center \includegraphics[width=4in]{manual_figures/ABKT_II_LHZ_seis.pdf}
\caption{\label{fg:test_data}
Windowing results for the test data set, plotted using the {\tt ./plot\_seismos\_gmt.sh} script.
}
@@ 171,10 +117,8 @@
subroutines in {\tt io\_subs.f90}.
\section{Scripts}
Several plotting routines ({\tt plot\_*.sh}) are provided as examples for
plotting seismograms, measurements and adjoint sources. All plotting is
done in gmt. These scripts will need to be modified to suit your
particular plotting needs.
+Several plotting routines ({\tt plot\_*.sh}) are provided in the {\tt scripts} subdirectory as examples for plotting seismograms, measurements and adjoint sources. All plotting is
+done using GMT (Generic Mapping Tools). These scripts will need to be modified to suit your particular plotting needs.
The script {\tt extract\_event\_windowing\_stats.sh} extracts statistical
information on the window selection process, on the measurements. Again,
Modified: seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_tuning.tex
===================================================================
 seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_tuning.tex 20090319 04:46:00 UTC (rev 14388)
+++ seismo/3D/ADJOINT_TOMO/flexwin_paper/latex/manual_tuning.tex 20090319 10:11:30 UTC (rev 14389)
@@ 2,19 +2,21 @@
FLEXWIN is adapted to your specific problem by modifying the values of the parameters in Table~\ref{tb:params}, and the functional form of those parameters that are timedependent. We consider the algorithm to be correctly adapted when false positives (windows around undesirable features of the seismogram) are minimized, and true positives (window around desirable features) are maximized. The choice of what makes an adequate set of windows remains subjective, as it depends strongly on the quality of the input model, the quality of the data, and the region of the Earth the tomographic inversion aims to constrain.
The base values of the various parameters are set in the {\tt PAR\_FILE}, which is read at run time. The functional forms of the time dependent parameters may be adjusted by modifying {\tt user\_parameters.f90}, and recompiling the code.
+The base values of the various parameters are set in the {\tt PAR\_FILE}, which is read at run time. Examples of base parameter values for the three tomographic scenarios discussed by \cite{MaggiEtal2009} can be found in Table~\ref{tb:example_params}. The functional forms of the time dependent parameters may be adjusted by modifying {\tt user\_parameters.f90} (see next section), and recompiling the code.
\begin{table}
\begin{tabular}{lp{0.8\linewidth}}
\hline
\multicolumn{2}{l}{Standard tuning parameters:} \\[5pt]
$T_{0,1}$ & bandpass filter corner periods \\
+$T_{0,1}$ & bandpass 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 shortterm:longterm ratio \\
$CC_0(t)$ & acceptance level for normalized crosscorrelation\\
+$\mathrm{CC}_0(t)$ & acceptance level for normalized crosscorrelation\\
$\Delta\tau_0(t)$ & acceptance level for time lag \\
$\Delta\ln{A}_0(t)$ & acceptance level for amplitude ratio \\
+$\Delta\tau_{\rm ref}$ & reference time lag \\
+$\Delta\ln{A}_{\rm ref}$ & reference amplitude ratio \\
\hline
\multicolumn{2}{l}{Fine tuning parameters:} \\ [5pt]
$c_0$ & for rejection of internal minima \\
@@ 22,7 +24,7 @@
$c_2$ & for rejection of unprominent 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 nonoverlapping window combination \\
+$w_{\mathrm{CC}}\quad w_{\rm len}\quad w_{\rm nwin}$ & for selection of best nonoverlapping window combination \\
\hline
\end{tabular}
\caption{\label{tb:params}
@@ 32,18 +34,20 @@
}
\end{table}
\section{User modifiable parameters}
The main usermodifiable parameters in the {\tt PAR\_FILE} are:
+\section{User parameters}
+The main user parameters in the {\tt PAR\_FILE} are:
\begin{description}
\item[{\tt WIN\_MIN\_PERIOD}]Corresponds to $T_0$ in Table~\ref{tb:params}, and is the short wavelength cutoff for the bandpass filter applied to the raw synthetic and observed seismograms.
\item[{\tt WIN\_MAX\_PERIOD}]Corresponds to $T_1$ in Table~\ref{tb:params}, and is the long wavelength cutoff for the bandpass filter applied to the raw synthetic and observed seismograms.
\item[{\tt SNR\_INTEGRATE\_BASE}]Corresponds to $r_P$ in Table~\ref{tb:params}, and is the minimum signal to noise ratio on the power of the observed seismogram for windowing to continue.
\item[{\tt SNR\_MAX\_BASE}]Corresponds to $r_A$ in Table~\ref{tb:params}, and is the minimum signal to noise ratio on the modulus of the observed seismogram for windowing to continue.
\item[{\tt WINDOW\_AMP\_BASE}]Corresponds to $r_0$ in Table~\ref{tb:params}, and is the minimum signal to noise ratio for a window on the observed seismogram to be acceptable.
+\item[{\tt WINDOW\_S2N\_BASE}]Corresponds to $r_0$ in Table~\ref{tb:params}, and is the minimum signal to noise ratio for a window on the observed seismogram to be acceptable.
\item[{\tt STALTA\_BASE}]Corresponds to $w_E$ in Table~\ref{tb:params}, and is the water level to be applied to the synthetic shortterm/longterm average waveform in order to generate candidate time windows. See Figure~\ref{fg:win_composite}a.
\item[{\tt CC\_BASE}]Corresponds to $CC_0$ in Table~\ref{tb:params}, and is the minimum normalized crosscorrelation value between synthetic and observed seismogram for a window to be acceptable.
\item[{\tt TSHIFT\_BASE}]Corresponds to $\Delta\tau_0$ in Table~\ref{tb:params}, and is the maximum crosscorrelation lag (in seconds) between synthetic and observed seismogram for a window to be acceptable.
\item[{\tt DLNA\_BASE}]Corresponds to $\Delta\ln{A}_0$ in Table~\ref{tb:params}, and is the maximum amplitude ratio ($\Delta\ln{A}$ or $\Delta A/A$) between synthetic and observed seismogram for a window to be acceptable.
+\item[{\tt TSHIFT\_REFERENCE}]Corresponds to $\Delta\tau_{\rm ref}$ in Table~\ref{tb:params}, and allows for a systematic traveltime bias in the synthetics.
+\item[{\tt TSHIFT\_REFERENCE}]Corresponds to $\Delta\ln{A}_{\rm ref}$ in Table~\ref{tb:params}, and allows for a systematic amplitude bias in the synthetics.
\item[{\tt C\_0}]Corresponds to $C_0$ in Table~\ref{tb:params}, and is expressed as a multiple of $w_E$. No window may contain a local minimum in its STA:LTA waveform that falls below the local value of $C_0 w_E$. See Figure~\ref{fg:win_composite}b.
\item[{\tt C\_1}]Corresponds to $C_1$ in Table~\ref{tb:params}, and is expressed as a multiple of $T_0$. No window may be shorter than $C_1 T_0$.
\item[{\tt C\_2}]Corresponds to $C_2$ in Table~\ref{tb:params}, and is expressed as a multiple of $w_E$. A window whose seed maximum on the STA:LTA waveform rises less than $C_2 w_E$ above either of its adjacent minima is rejected. See Figure~\ref{fg:win_composite}c.
@@ 56,6 +60,38 @@
\item[{\tt WEIGHT\_N\_WINDOWS}]Corresponds to $w_{\rm nwin}$ in Table~\ref{tb:params}, and is the weight given to the total number of windows in the process of resolving window overlaps.
\end{description}
+\begin{table}
+\begin{center}
+\begin{tabular}{lcccccc}
+\hline
+ & Global & \multicolumn{2}{c}{Japan} & \multicolumn{3}{c}{S. California} \\
+\hline
+$T_{0,1}$ & 50, 150 & 24, 120 & 6, 30 & 6, 30 & 3, 30 & 2, 30 \\
+$r_{P,A}$ & 3.5, 3.0 & 3.5, 3.0 & 3.5, 3.0 & 3.0, 2.5 & 2.5, 3.5 & 2.5, 3.5 \\
+$r_0$ & 2.5 & 1.5 & 3.0 & 3.0 & 4.0 & 4.0 \\
+$w_E$ & 0.08 & 0.11 & 0.12 & 0.18 & 0.11 & 0.07 \\
+$\mathrm{CC}_0$ & 0.85 & 0.70 & 0.73 & 0.71 & 0.80 & 0.85 \\
+$\Delta\tau_0$ & 15 & 12.0 & 3.0 & 8.0 & 4.0 & 3.0 \\
+$\Delta\ln{A}_0$ & 1.0 & 1.0 & 1.5 & 1.5 & 1.0 & 1.0 \\
+$\Delta\tau_{\rm ref}$ & 0.0 & 0.0 & 0.0 & 4.0 & 2.0 & 1.0 \\
+$\Delta\ln{A}_{\rm ref}$& 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
+\hline
+$c_0$ & 0.7 & 0.7 & 0.7 & 0.7 & 1.3 & 1.0 \\
+$c_1$ & 4.0 & 3.0 & 3.0 & 2.0 & 4.0 & 5.0 \\
+$c_2$ & 0.3 & 0.0 & 0.6 & 0.0 & 0.0 & 0.0 \\
+$c_{3a,b}$ & 1.0, 2.0 & 1.0, 2.0 & 1.0, 2.0 & 3.0, 2.0 & 4.0, 2.5 & 4.0, 2.5 \\
+$c_{4a,b}$ & 3.0, 10.0 & 3.0, 25.0 & 3.0, 12.0 & 2.5, 12.0 & 2.0, 6.0 & 2.0, 6.0 \\
+$w_{\mathrm{CC}}, w_{\rm len}, w_{\rm nwin}$
+ & 1, 1, 1 & 1, 1, 1 & 1, 1, 1 & 0.5,1.0,0.7 & 0.70,0.25,0.05 & 1,1,1 \\
+\hline
+\end{tabular}
+\caption{\label{tb:example_params}
+Values of standard and finetuning parameters for the three seismological
+scenarios discussed \cite{MaggiEtal2009}. This table is identical to Table~3 of that study.
+}
+\end{center}
+\end{table}
+
\begin{figure}
\center \includegraphics[width=6in]{figures/fig/window_composite.pdf}
\caption{\label{fg:win_composite}
@@ 105,7 +141,7 @@
\section{Time dependence of user parameters}
A subset of the FLEXWIN parameters from Table~\ref{tb:params} are timedependent (where time is measured along the seismogram). This feature enables the user to exercise fine control of the windowing algorithm. The user can modulate the timedependence of these parameters by editing the {\tt set\_up\_criteria\_arrays} subroutine in the {\tt user\_functions.f90} file. This subroutine is called after the seismograms have been read in, and the following variables have been set:
\begin{description}
\item[{\tt npts, dt, b, npts}] Number of points, time step, time of first point with respect to the reference time of both seismograms. The observed and synthetic seismograms should have identical values of these three quantities.
+\item[{\tt npts, dt, b}] Number of points, time step, time of first point with respect to the reference time of both seismograms. The observed and synthetic seismograms should have identical values of these three quantities.
\item[{\tt evla, evlo, evdp, stla, stlo}] Event latitude, event longitude, event depth (km), station latitude, station longitude, read from the observed seismogram.
\item[{\tt azimuth, backazimuth, dist\_deg, dist\_km}] Calculated from the event and station locations above.
\item[{\tt kstnm, knetwk, kcmpnm}] Station name, network name, component name, read from the observed seismogram.
@@ 181,3 +217,204 @@
The above examples illustrate the power of the {\tt user\_functions.f90} file. The user can choose to include/exclude any portion of the seismogram, and to make the rejection criteria for windows more or less stringent on any other portion of the seismogram. All the seismogramdependent variables whose values are known when the {\tt set\_up\_criteria\_arrays} subroutine is executed may be used to inform these choices, leading to an infinite number of windowing possibilities. The careful user will use knowledge of the properties of the observed data set, the limitations of the synthetic waveforms, and the final use to which the selected windows will be put in order to tailor the subroutine to the needs of each study.
For a given set of data and synthetics, the {\tt PAR\_FILE} and {\tt user\_functions.f90} files uniquely determine the windowing results.
+
+\subsection{Examples of user functions\label{ap:user_fn}}
+
+As concrete examples of how the time dependence of the tuning parameters can be exploited, we present here the functional forms of the time dependencies used for the three example tomographic scenarios (global, Japan and southern California) described in \cite{MaggiEtal2009}. In each example we use predicted arrival times derived from 1D Earth models to help modulate certain parameters. Note, however, that the actual selection of individual windows is based on the details of the waveforms, and not on information from 1D Earth models.
+
+\subsubsection{Global scenario\label{ap:user_global}}
+
+In the following, $h$ indicates earthquake depth, $t_Q$ indicates the approximate start of the Love wave predicted by a group wave speed of 4.2~\kmps, and $t_R$ indicates the approximate end of the Rayleigh wave predicted by a group wave speed of 3.2~\kmps. In order to reduce the number of windows picked beyond R1, and to ensure that those selected beyond R1 are a very good match to the synthetic waveform, we raise the water level on the STA:LTA waveform and impose stricter criteria on the signaltonoise ratio and the waveform similarity after the approximate end of the surfacewave arrivals. We allow greater flexibility in crosscorrelation time lag $\Delta\tau$ for intermediate depth and deep earthquakes. We lower the crosscorrelation value criterion for surfacewaves in order to retain windows with a slight mismatch in dispersion characteristics.
+
+We therefore use the following time modulations:
+\begin{align}
+w_E(t) & =
+ \begin{cases}
+ w_E \text{$t \leq t_R$} ,\\
+ 2 w_E \text{$t > t_R$},
+ \end{cases}
+\\
+r_0(t) & =
+ \begin{cases}
+ r_0 & \text{$t \leq t_R$}, \\
+ 10r_0 & \text{$t > t_R$} ,
+ \end{cases}
+\\
+\mathrm{CC}_0(t) & =
+ \begin{cases}
+ \mathrm{CC}_0 & \text{$t \leq t_R$}, \\
+ 0.9 \mathrm{CC}_0 & \text{$t_Q < t \leq t_R$}, \\
+ 0.95 & \text{$t > t_R$} ,
+ \end{cases}
+\\
+\Delta\tau_0(t) & =
+ \begin{cases}
+ \begin{cases}
+ \tau_0 & \text{$t \leq t_R$}, \\
+ \tau_0/3 & \text{$t > t_R$} ,
+ \end{cases}
+ & \text{$h \leq$ 70~km} \\
+ 1.4\tau_0 & \text{70~km $< h <$ 300~km}, \\
+ 1.7\tau_0 & \text{$h \geq$ 300~km},
+ \end{cases}
+ \\
+\Delta \ln A_0(t) & =
+ \begin{cases}
+ \Delta \ln A_0 & \text{$t \leq t_R$}, \\
+ \Delta \ln A_0/3 & \text{$t > t_R$} .
+ \end{cases}
+\end{align}
+
+%
+
+\subsubsection{Japan scenario\label{ap:user_japan}}
+In the following, $t_P$ and $t_S$ denote the start of the time windows for $P$ and $S$ waves, as predicted by the 1D IASPEI91 model \citep{KennettEngdahl1991}, and $t_{R1}$ indicates the end of the surfacewave time window. For the \trange{24}{120} data, we consider the waveform between the start of the $P$ wave to the end of the surfacewave. We therefore modulate $w_E(t)$ as follows:
+
+%
+\begin{align}
+w_E(t) & =
+ \begin{cases}
+ 10 w_E & \text{$t < t_P$}, \\
+ w_E & \text{$t_P \le t \leq t_{R1}$}, \\
+ 10 w_E & \text{$t > t_{R1}$}.
+ \end{cases}
+\end{align}
+
+For the \trange{6}{30} data, the fit between the synthetic and observed surfacewaves is expected to be poor, as the 3D model used to calculate the synthetics cannot produce the required complexity. We therefore want to concentrate on bodywave arrivals only, and avoid surfacewave windows altogether by modulating $w_E(t)$ as follows:
+%
+\begin{align}
+w_E(t) & =
+ \begin{cases}
+ 10 w_E & \text{$t < t_P$}, \\
+ w_E & \text{$t_P \le t \leq t_S$}, \\
+ 10 w_E & \text{$t > t_S$}.
+ \end{cases}
+\end{align}
+
+We use constant values of $r_0(t)=r_0$, $\mathrm{CC}_0(t)=\mathrm{CC}_0$ and $\Delta \ln A_0(t)=\Delta \ln A_0$ for both period ranges. In order to allow greater flexibility in crosscorrelation time lag $\Delta\tau$ for intermediate depth and deep earthquakes we use:
+
+\begin{align}
+\Delta\tau_0(t) & =
+ \begin{cases}
+ 0.08 \text{$t_P$} & \text{$h \leq$ 70~km}, \\
+ \max(0.05 \text{$t_P$}, 1.4\tau_0) & \text{70~km $< h <$ 300~km}, \\
+ \max(0.05 \text{$t_P$}, 1.7\tau_0) & \text{$h \geq$ 300~km}.
+ \end{cases}
+\end{align}
+%
+
+\subsubsection{Southern California scenario\label{ap:user_socal}}
+
+In the following, $t_P$ and $t_S$ denote the start of the time windows for the crustal P wave and the crustal S wave, computed from a 1D layered model appropriate to Southern California \citep{Wald95}. The start and end times for the surfacewave time window, $t_{R0}$ and $t_{R1}$, as well as the criteria for the time shifts $\Delta\tau_0(t)$, are derived from formulas in \cite{KomatitschEtal2004}. The sourcereceiver distance (in km) is denoted by $\Delta$.
+
+%CHT modified
+
+For the \trange{6}{40} and \trange{3}{40} data, we use constant values of $r_0(t)=r_0$, $\mathrm{CC}_0(t)=\mathrm{CC}_0$, $\Delta\tau_0(t)=\Delta\tau_0$, and $\Delta \ln A_0(t)=\Delta \ln A_0$. We exclude any arrivals before the $P$ wave and after the Rayleigh wave. This is achieved by the boxcar function for $w_E(t)$:
+%
+\begin{align}
+w_E(t) & =
+ \begin{cases}
+ 10 w_E & \text{$t < t_P$}, \\
+ w_E & \text{$t_P \le t \leq t_{R1}$}, \\
+ 10 w_E & \text{$t > t_{R1}$},
+ \end{cases}
+\end{align}
+%For the \trange{6}{40} data, we exclude any arrivals before the $P$ wave and reduce the number of windows picked beyond R1 by modulating $w_E(t)$. We use constant values of $r_0(t)=r_0$, $\mathrm{CC}_0(t)=\mathrm{CC}_0$ and $\Delta \ln A_0(t)=\Delta \ln A_0$, but modulate the crosscorrelation time lag criterion so that it is less strict at larger epicentral distances and for surfacewaves. We therefore use:
+%
+%\begin{align}
+%w_E(t) & =
+% \begin{cases}
+% 10 w_E & \text{$t < t_P$}, \\
+% w_E & \text{$t_P \le t \leq t_{R1}$}, \\
+% 2 w_E & \text{$t > t_{R1}$},
+% \end{cases}
+%\\
+%\Delta\tau_0(t) & =
+% \begin{cases}
+% 3.0 + \Delta/80.0 & \text{$t \le t_{R0}$}, \\
+% 3.0 + \Delta/50.0 & \text{$t > t_{R0}$},
+% \end{cases}
+%\end{align}
+
+For the \trange{2}{40} data, we avoid selecting surfacewave arrivals as the 3D model used to calculate the synthetics cannot produce the required complexity. The waterlevel criteria then becomes:
+
+\begin{align}
+w_E(t) & =
+ \begin{cases}
+ 10 w_E & \text{$t < t_P$}, \\
+ w_E & \text{$t_P \le t \leq t_S$}, \\
+ 10 w_E & \text{$t > t_S$}.
+ \end{cases}
+%\\
+%\Delta\tau_0(t) & = \Delta\tau_0.
+\end{align}
+
+
+%
+
+
+\section{Tuning considerations}
+FLEXWIN is not a blackbox application, and as such cannot be applied blindly
+to any given dataset or tomographic scenario. The data windowing required by
+any given problem will differ depending on the inversion method, the scale of
+the problem (local, regional, global), the quality of the data set and that of
+the model and method used to calculate the synthetic seismograms. The user
+must configure and tune the algorithm for the given problem. Here we
+shall discuss general considerations the user should bear in mind during
+the tuning process.
+
+We suggest the following as a practical starting sequence for tuning the algorithm
+(the process may need to be repeated and refined several times before
+converging on the optimal set of parameters for a given problem and dataset).
+
+$T_{0,1}$ : In setting the corner periods of the bandpass filter, the
+user is deciding on the frequency content of the information to be used in the
+tomographic problem. 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. The
+frequency content in the data depends on the spectral characteristics of the
+source, on the instrument responses, and on the attenuation
+characteristics of the medium. As $T_{0,1}$ depend on the source and station
+characteristics, which may be heterogeneous in any given dataset, these filter
+periods can be modified dynamically by constructing an appropriate user
+function (e.g. {\em if station is in list of stations with instrument X then
+reset T0 and T1 to new values}).
+
+$r_{P,A}$ : In setting the signaltonoise ratios for the entire seismogram the
+user is applying a simple quality control on the data. Note that these criteria
+are applied after filtering. No windows will be defined on data that fail this
+quality control.
+
+$w_E(t)$ : The shortterm average longterm average ratio $E(t)$ of a constant signal
+converges to a constant value when
+the length of the timeseries is greater than the effective averaging length of
+the longterm average. This value is 0.08 for the shortterm average longterm average ratio used in FLEXWIN (it has a small dependence on $T_0$, which can be ignored in most applications). We suggest the user start with a constant
+level for $w_E(t)$ equal to this convergence value. The time dependence of
+$w_E(t)$ should then be adjusted to exclude those portions of the waveform the
+user is not interested in, by raising $w_E(t)$ (e.g. to exclude the fundamental
+mode surfacewave: {\em if t $>$ fundamental mode surfacewave arrival time then set $w_E(t)=1$}).
+We suggest finer adjustments to $w_E(t)$ be made after $r0(t)$,
+$CC_0(t)$, $\Delta T_0(t)$ and $\Delta \ln A_0(t)$ have been configured.
+
+$r_0(t)$, $\mathrm{CC}_0(t)$, $\Delta \tau_{\rm ref}$, $\Delta
+\tau_0(t)$, $\Delta \ln A_{\rm ref}$ and $\Delta \ln A_0(t)$ : These parameters 
+window signaltonoise ratio, normalized crosscorrelation value between
+observed and synthetic seismograms, crosscorrelation time lag, and amplitude
+ratio  control the degree of wellbehavedness of the data within accepted
+windows. The user first sets constant values for these four parameters, then
+adds a time dependence if required. Considerations that should be taken into
+account include the quality of the Earth model used to calculate the synthetic
+seismograms, the frequency range, the dispersed nature of certain arrivals (e.g.
+{\em for t corresponding to the group velocities of surfacewaves, reduce
+$CC_0(t)$}), and {\em a priori} preferences for picking certain smallamplitude seismic phases
+(e.g. {\em for t close to the expected arrival for $P_{\rm diff}$, reduce $r_0(t)$}).
+$\Delta \tau_{\rm ref}$ and $\Delta \ln A_{\rm ref}$ should be set to zero at first, and only
+reset if the synthetics contain a systematic bias in traveltimes or amplitudes.
+
+
+$c_{04}$ : These parameters control the process by which the suite of all possible data windows is pared down using criteria on the shape of the STA:LTA $E(t)$ waveform alone. We suggest the user start by setting these values to those used in our global example (see Table~\ref{tb:example_params}). Subsequent minimal tuning should be performed by running the algorithm on a subset of the data and closely examining the lists of windows rejected at each stage to make sure the user agrees with the choices made by the algorithm.
+
+$w_{\mathrm{CC}}$, $w_{\rm len}$ and $w_{\rm nwin}$ : These parameters control the overlap resolution stage of the algorithm. Values of $w_{\mathrm{CC}}= w_{\rm len} = w_{\rm nwin} = 1$ should be reasonable for most applications.
+
+The objective of the tuning process summarily described here should be to maximize the selection of windows around desirable features in the seismogram, while minimizing the selection of undesirable features, bearing in mind that the desirability or undesirability of a given feature is subjective, and depends on how the user subsequently intends to use the information contained within he data windows.
+
More information about the CIGCOMMITS
mailing list