[cig-commits] r13366 - in seismo/3D/ADJOINT_TOMO/flexwin: . latex
carltape at geodynamics.org
carltape at geodynamics.org
Thu Nov 20 18:28:44 PST 2008
Author: carltape
Date: 2008-11-20 18:28:43 -0800 (Thu, 20 Nov 2008)
New Revision: 13366
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/latex/acknowledgements.tex
seismo/3D/ADJOINT_TOMO/flexwin/latex/figures_paper.tex
seismo/3D/ADJOINT_TOMO/flexwin/latex/flexwin_paper.pdf
seismo/3D/ADJOINT_TOMO/flexwin/latex/method.tex
seismo/3D/ADJOINT_TOMO/flexwin/latex/results.tex
seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90
Log:
Minor paper edits. Also comments for xcorr_calc in seismo_subs.f90.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/acknowledgements.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/latex/acknowledgements.tex 2008-11-21 01:26:10 UTC (rev 13365)
+++ seismo/3D/ADJOINT_TOMO/flexwin/latex/acknowledgements.tex 2008-11-21 02:28:43 UTC (rev 13366)
@@ -8,3 +8,4 @@
We thank the Hi-net 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 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.
+{\bf CHT mod:} Thanks reviewers? Acknowledge SURF funding for Daniel?
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/figures_paper.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/latex/figures_paper.tex 2008-11-21 01:26:10 UTC (rev 13365)
+++ seismo/3D/ADJOINT_TOMO/flexwin/latex/figures_paper.tex 2008-11-21 02:28:43 UTC (rev 13366)
@@ -46,7 +46,7 @@
% CHECK THAT THE MOMENT IS LISTED IN N-M, NOT DYNE-CM
% 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 \\
-200808270646A & -10.49 & 41.4400 & 24.0 & 4.68e+17 & 5.7 & Comoros Region \\
+200808270646A & -10.49 & 41.44 & 24.0 & 4.68e17 & 5.7 & Comoros Region \\
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
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/flexwin_paper.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/method.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/latex/method.tex 2008-11-21 01:26:10 UTC (rev 13365)
+++ seismo/3D/ADJOINT_TOMO/flexwin/latex/method.tex 2008-11-21 02:28:43 UTC (rev 13366)
@@ -6,7 +6,7 @@
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. All the synthetic seismograms presented in
-this paper have been generated using the SPECFEM3D package \citep{KomatitschEtal2002}.
+this paper have been generated using the SPECFEM3D package \citep{KomatitschEtal2002,KomatitschEtal2004}.
The window selection process has five stages, each of which is discussed individually
below: {\em \stga:} pre-processing; {\em \stgb:} definition of preliminary
@@ -251,6 +251,28 @@
synthetic seismograms within these windows, and to reject
those that fail basic fit-based criteria.
+{\bf CHT mod:}
+In this section we consider records that have been windowed as
+%
+\begin{eqnarray}
+d_k(t) &=& W_k(t)\,d(t)
+\\
+s_k(t) &=& W_k(t)\,s(t)
+\end{eqnarray}
+%
+where $k$ denotes the index of the measurement window in the full dataset, and $W_k(t)$ is the corresponding boxcar window between times $t_{1k}$ and $t_{2k}$:
+%
+\begin{align}
+W_k(t) &=
+ \begin{cases}
+ 0 & \text{$t < t_{1k}$}, \\
+ 1 & \text{$t_{1k} \le t \leq t_{2k}$}, \\
+ 0 & \text{$t > t_{2k}$}.
+ \end{cases}
+\end{align}
+%
+To avoid clutter, we will omit the $k$-subscript, but the reader is reminded that the data and synthetics are zeroed outside the window $[t_1,\,t_2]$.
+
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 $\mathrm{CC}$,
@@ -260,6 +282,7 @@
$A_{\rm noise}$ are the maximum absolute values of the observed seismogram $|d(t)|$ in the window
and in the noise time-span respectively (the noise time-span is the same as that for equation~\ref{eq:noise}). The cross-correlation value $\mathrm{CC}$ is defined as the maximum value of the
cross-correlation function $\mathrm{CC}={\rm max}[\Gamma(t)]$, where
+{\bf CHT: I DO NOT THINK THAT THIS REFLECTS THE IMPLEMENTATION.}
\begin{equation}
\Gamma(t) = \frac {\int s(t'-t)d(t')\,\rmd t'}{ \left[\int s^2(t')\,\rmd t' \int d^2(t')\,\rmd t' \right]^{1/2}}
\end{equation}
@@ -275,7 +298,7 @@
%\end{equation}
\begin{equation}
\Delta\ln{A} = \ln(A_{\rm obs}/A_{\rm syn})
-= 0.5 \ln \left[ \frac{\int d(t)^2 \,\rmd t}{\int s(t)^2 \,\rmd t} \right] \label{eq:dlnA_def}
+= 0.5 \ln \left[ \frac{\int d^2(t) \,\rmd t}{\int s^2(t) \,\rmd t} \right] \label{eq:dlnA_def}
\end{equation}
%
(note that \citet[][Eq.~3]{DahlenBaig2002} is the first-order approximation of equation~\ref{eq:dlnA_def}).
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/results.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/latex/results.tex 2008-11-21 01:26:10 UTC (rev 13365)
+++ seismo/3D/ADJOINT_TOMO/flexwin/latex/results.tex 2008-11-21 02:28:43 UTC (rev 13366)
@@ -272,7 +272,7 @@
We test the windowing code using three period ranges: \trange{6}{30}, \trange{3}{30}, and \trange{2}{30}. The parameters we use for the windowing code are listed in Table~\ref{tb:example_params}. Figures~\ref{fg:socal_CLC} and~\ref{fg:socal_FMP} show examples of the output from the windowing algorithm for event 9818433 listed in Table~\ref{tb:events} recorded at two different stations, while Figure~\ref{fg:socal_rs_T06} shows a summary plot for event 9983429 in the \trange{6}{30} period range.
The windowing algorithm tends to identify five windows on each set of three-component longer-period seismograms (Figures~\ref{fg:socal_CLC} and~\ref{fg:socal_rs_T06}): on the vertical and radial components the first window corresponds to the body-wave arrival and the second to the Rayleigh wave, while windows on the transverse component capture the Love wave.
-The shorter-period synthetic seismograms do not agree well with the observed seismograms, especially in the later part of the signal, leading to fewer picked windows. In Figure~\ref{fg:socal_CLC}c, only three windows are selected by the algorithm: the P arrival recorded on the radial component, the S arrival on the transverse component, and the Love-wave arrival on the transverse component. The P-wave arrival on the vertical component is rejected because the cross-correlation value within the time window did not exceed the specified minimum value of 0.85 (Table~\ref{tb:example_params}). The small peak at 38~s in the STA:LTA curve for the transverse component identifies P-wave energy that does not have a large enough signal-to-noise ratio to be picked. However, it highlights the possibility of measuring subtle phases that may be present in 3D synthetics.
+The shorter-period synthetic seismograms do not agree well with the observed seismograms, especially in the later part of the signal, leading to fewer picked windows. In Figure~\ref{fg:socal_CLC}c, only three windows are selected by the algorithm: the P arrival recorded on the radial component, the S arrival on the transverse component, and the Love-wave arrival on the transverse component. The P arrival (PmP or Pn) in fact appears on all three components on both data and synthetics. On the vertical component it is rejected because the cross-correlation value within the time window did not exceed the specified minimum value of 0.85 (Table~\ref{tb:example_params}). On the transverse component it does not have a large enough signal-to-noise ratio to be picked, but it is evident as a small peak at 36~s in the STA:LTA curve, and more conspicuous when zooming into the synthetics and data. The presence of the P arrival on the transverse component highlights the possibility of measuring subtle phases that may be present in 3D synthetics.
Figure~\ref{fg:socal_FMP} shows results for the same event as Figure~\ref{fg:socal_CLC}, but for a different station, FMP, situated 52~km from the event and within the Los Angeles basin. Comparison of the two figures highlights the characteristic resonance caused by the thick sediments within the basin. This resonance is beautifully captured by the transverse component synthetics (Figure~\ref{fg:socal_FMP}b, record T), thanks to the inclusion of the basin in the model \citep{KomatitschEtal2004}. In order to pick such long time windows with substantial frequency-dependent measurement differences, we are forced to lower the minimum cross-correlation value $\mathrm{CC}_0$ for the entire dataset (0.71 in Table~\ref{tb:example_params}) and increase $c_{4b}$ to capture the slow decay in the STA:LTA curves (Figure~\ref{fg:socal_FMP}b, record T). It is striking that although these arrivals look nothing like the energy packets typical for the global case, the windowing algorithm is still able to determine the proper start and end times for the windows. In Figure~\ref{fg:socal_FMP}c the windowing algorithm selects three short-period body-wave time windows with superb agreement between data and synthetics.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90 2008-11-21 01:26:10 UTC (rev 13365)
+++ seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90 2008-11-21 02:28:43 UTC (rev 13366)
@@ -675,7 +675,7 @@
! local variables
integer :: nlen
integer :: i_left, i_right, i, j, id_left, id_right
- double precision :: cc, norm
+ double precision :: cc, norm, norm_s
! initialise shift and cross correlation to zero
ishift = 0
@@ -689,28 +689,49 @@
! length of window (number of points, including ends)
nlen = i2 - i1 + 1
+ ! power of synthetic signal in window
+ norm_s = sqrt(sum(s(i1:i2)*s(i1:i2)))
+
! left and right limits of index (time) shift search
! NOTE: This looks OUTSIDE the time window of interest to compute TSHIFT and CC.
- ! THIS IS A VERY IMPORTANT DECISION!
+ ! How far to look outside, in theory, should be another parameter.
+ ! However, it does not matter as much if the data and synthetics are
+ ! zeroed outside the windows, as currently done in calc_criteria.
i_left = -1*int(nlen/2.0)
i_right = int(nlen/2.0)
- ! i -> shift (to be applied to DATA (d) in cc search)
+ ! i is the index to shift to be applied to DATA (d)
do i = i_left, i_right
+
+ ! normalization factor varies as you take different windows of d
+ id_left = max(1,i1+i) ! left index for data window
+ id_right = min(npts,i2+i) ! right index for data window
+ norm = norm_s * sqrt(sum(d(id_left:id_right)*(d(id_left:id_right))))
+
+ ! cc as a function of i
cc = 0.
- id_left = max(1,i1+i) ! left-most point on the data that will be treated
- id_right = min(npts,i2+i) ! right-most point on the data that will be treated
- norm = sqrt(sum(s(i1:i2)*s(i1:i2)) * sum(d(id_left:id_right)*(d(id_left:id_right))))
- do j = i1, i2
- if((j+i).ge.1 .and. (j+i).le.npts) cc = cc + s(j)*d(j+i)
+ do j = i1, i2 ! loop over full window length
+ if((j+i).ge.1 .and. (j+i).le.npts) cc = cc + s(j)*d(j+i) ! d is shifted by i
enddo
cc = cc/norm
+
+ ! keeping cc-max only
if (cc .gt. cc_max) then
cc_max = cc
ishift = i
endif
enddo
+ ! EXAMPLE: consider the following indexing:
+ ! Two records are from 1 to 100, window is i1=20 to i2=41.
+ ! --> nlen = 22, i_left = -11, i_right = 11
+ ! i i1+i i2+i id_left id_right
+ ! -11 9 30 9 30
+ ! -5 15 36 15 36
+ ! 0 20 41 20 41
+ ! 5 25 46 25 46
+ ! 10 31 52 31 52
+
end subroutine xcorr_calc
!------------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list