[cig-commits] r13392 - seismo/3D/ADJOINT_TOMO/flexwin/latex
carltape at geodynamics.org
carltape at geodynamics.org
Tue Nov 25 10:06:27 PST 2008
Author: carltape
Date: 2008-11-25 10:06:27 -0800 (Tue, 25 Nov 2008)
New Revision: 13392
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/latex/cc_notes.pdf
seismo/3D/ADJOINT_TOMO/flexwin/latex/cc_notes.tex
Log:
Updated CC notes.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/cc_notes.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/flexwin/latex/cc_notes.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/latex/cc_notes.tex 2008-11-25 14:40:32 UTC (rev 13391)
+++ seismo/3D/ADJOINT_TOMO/flexwin/latex/cc_notes.tex 2008-11-25 18:06:27 UTC (rev 13392)
@@ -1,4 +1,4 @@
-\documentclass[11pt]{article}
+\documentclass[10pt]{article}
%\documentclass[extra]{gji}
\usepackage{graphicx}
\usepackage{rotating}
@@ -22,28 +22,66 @@
\section*{Notes on normalized cross-correlation}
-Please see both pages.
-
\subsection*{Cross-correlation formula}
-First, we list the following formula:
+In the paper, we have the following.
+
+The cross-correlation value $\mathrm{CC}$ is defined as the maximum value of the normalised
+cross-correlation function, $\mathrm{CC}={\rm max}[\Gamma(t)]$, where
%
\begin{equation}
\Gamma(t) = \frac {\int \widetilde{s}(t') \widetilde{d}(t'-t)\,\rmd t'}
{ \left[\int \widetilde{s}^2(t')\,\rmd t' \int \widetilde{d}^2(t'-t)\,\rmd t' \right]^{1/2}}
- \label{eq:paper}
+\label{eq:paper}
\end{equation}
%
-My point about this has nothing to do with what we implement. It's that, in the absence of integration limits for convolution-like expressions, one would assume $[-\infty,\infty]$ limits. In that case, it doesn't matter how you shift $d^2(t)$, because the result is the same:
+quantifies the similarity in shape between the $\widetilde{s}(t)$ and $\widetilde{d}(t)$ waveforms, and the integration limits are the start and end times of the window. The time lag $\Delta \tau$ is defined as the value of $t$ at which $\Gamma$ is maximal, and quantifies the delay in time between a synthetic and observed phase arrival.
+
+%My point about this has nothing to do with what we implement. It's that, in the absence of integration limits for convolution-like expressions, one would assume $[-\infty,\infty]$ limits. In that case, it doesn't matter how you shift $d^2(t)$, because the result is the same:
+%%
+%\begin{equation}
+%\int_{\infty}^{\infty} \widetilde{d}^2(t'-t)\,\rmd t'
+%= \int_{\infty}^{\infty} \widetilde{d}^2(t'+t)\,\rmd t'
+%= \int_{\infty}^{\infty} \widetilde{d}^2(t')\,\rmd t'
+%\end{equation}
%
+%That is why the limits are important, and the limits are essentially included in Equation~(2) in the PDF that you sent.
+
+{\em AM:
+The effective integration limits in (\ref{eq:paper}) are the start and end times of the window $t_S$ and $t_E$ (which are times that have been chosen on the synthetic seismogram). In the numerator, as $\widetilde{s}(t)$ is zero outside these times, then integrating over $-\infty +\infty$ is the same as integrating over $t_S$ to $t_E$. The same is true for the $\widetilde{s}(t)$ integral in the denominator. For the $\widetilde{d}(t)$ integral in the denominator, using the synthetic $t_S$ and $t_E$ as integration limits essentially extracts that portion of the windowed data that overlaps in time with the synthetic. I have added a line in the text of the paper to say what the integration limits are.
+}
+
+Next, let's add some details to (\ref{eq:paper}), just for the notes. I put the $[t_S,t_E]$ limits on the integrations with $\widetilde{s}(t')$, even though they would be the same as using $[-\infty,\infty]$, since $\widetilde{s}(t')$ is zero outside $[t_S,t_E]$.
+%
\begin{equation}
-\int_{\infty}^{\infty} \widetilde{d}^2(t'-t)\,\rmd t'
-= \int_{\infty}^{\infty} \widetilde{d}^2(t'+t)\,\rmd t'
-= \int_{\infty}^{\infty} \widetilde{d}^2(t')\,\rmd t'
+\Gamma(t) = \frac {\int_{t_S}^{t_E} \widetilde{s}(t') \widetilde{d}(t'-t)\,\rmd t'}
+ { \left[\int_{t_S}^{t_E} \widetilde{s}^2(t')\,\rmd t' \int_{t_S}^{t_E} \widetilde{d}^2(t'-t)\,\rmd t' \right]^{1/2}}
+\label{eq:paper2}
\end{equation}
%
-That is why the limits are important, and the limits are essentially included in Equation~(2) in the PDF that you sent. What would (2) look like for us? In our case, I suppose we have
+I think I get it, finally. In the algorithm we slide the integration bounds (boxcar window) across $\widetilde{d}(t)$, but in the formula we slide $\widetilde{d}(t)$ passed the boxcar. Now in the implementation, we do not allow the entire window to slide past the other. So one way to write this would be:
+%
+\begin{align}
+\Gamma(t) &=
+ \begin{cases}
+ \frac {\int_{t_S}^{t_E} \widetilde{s}(t') \widetilde{d}(t'-t)\,\rmd t'}
+ { \left[\int_{t_S}^{t_E} \widetilde{s}^2(t')\,\rmd t' \int_{t_S}^{t_E} \widetilde{d}^2(t'-t)\,\rmd t' \right]^{1/2}},
+ & \text{$-F(t_S-t_E) \le t \le F(t_S-t_E)$}, \\
+ 0, & \text{otherwise}.
+ \end{cases}
+\end{align}
+%
+where $F$ is the fraction of the window that we allow the data to slide passed the synthetics.
+Do you agree with this?
+%-----------------------------------------
+
+\subsubsection*{CC formula from PDF}
+
+I will assume that we are not aiming to list or implement (\ref{ccnorm}). So I dumped it all into this subsection.
+
+What would (2) look like for us? In our case, I suppose we have
+
$d$ is the image and the sum over $t$ under the window containing the feature $s$ positioned at $t'$
And then (2) looks like
@@ -51,26 +89,25 @@
\begin{eqnarray}
\gamma(t') &=& \frac{\sum_{t}[d(t) - \bar{d}_{t'} ][s(t-t') - \bar{s}] }
{\{\sum_{t} [d(t) - \bar{d}_{t'}]^2 \sum_{t} [s(t-t') - \bar{s} ]^2 \}^{1/2} }
-\\
+\nonumber \\
\bar{d}_{t'} &=&
-\\
+\nonumber \\
\bar{s} &=&
+\label{ccnorm}
\end{eqnarray}
%
This is only a start. It is still sloppy, considering the mix of continuous and discrete notation, and the vague summation (integration) limits. If this is the formula we are trying to implement, would you please adapt it to what we want to use, and also include explicit integration limits? Then we can see how it is implemented.
-{\em AM:
-The effective integration limits in (\ref{eq:paper}) are the start and end times of the window $t_S$ and $t_E$ (which are times that have been chosen on the synthetic seismogram). In the numerator, as $\widetilde{s}(t)$ is zero outside these times, then integrating over $-\infty +\infty$ is the same as integrating over $t_S$ to $t_E$. The same is true for the $\widetilde{s}(t)$ integral in the denominator. For the $\widetilde{d}(t)$ integral in the denominator, using the synthetic $t_S$ and $t_E$ as integration limits essentially extracts that portion of the windowed data that overlaps in time with the synthetic. I have added a line in the text of the paper to say what the integration limits are.
-}
+%-----------------------------------------
-\pagebreak
\subsection*{Limiting time-shift searches in CC calculation on the basis of maximum allowed time-shift}
-\begin{quote}
+
+{\em AM:
I do not think this is a good idea for window selection. The time-shift should be the shift at which the maximum cross-correlation is achieved. Say your true maximum cross-correlation is at 5 seconds, but you have limited your search to 4 seconds: you will have a value of cross-correlation which would be maximal at 4 seconds, and you would therefore say the time-shift between your two signals is 4 s, when in fact it is 5. I would much prefer taking the solution you have in the measurement code, i.e. using the entire window length as the time-shift search limits, then rejecting those windows that have larger time-shifts than we demand. The fact that we work essentially with box-car tapered windows means there is no danger of bringing in information from outside the windows themselves.
-\end{quote}
+}
%
-First, your are right in one thing: it is not the energy outside the window. It's the energy {\em inside} the window. This will be an issue when a window is longer than the allowable time shift, which is always the case, I think.
+First, I agree, it is not the energy outside the window, it's the energy {\em inside} the window. This will be an issue when a window is longer than the allowable time shift, which is always the case, I think.
Okay, consider this. If the real time shift is 5~s, but you think it should generally be no more than 4~s, then you pick 4~s as the cross-correlation. But the max CC value is probably not going to have a passable value to have this window picked. And if it does, then I still think the 4~s shift is much better than rejecting the window. It would be better to use this in the inversion than nothing at all. All it means is using a 4~s weight on the adjoint source rather than a 5~s weight.
@@ -86,12 +123,100 @@
The normalization scheme we have used up to now (where we essentially consider only the windowed data that overlaps with the windowed synthetic) gives approximately equal weight to the cross-correlation values at every time-shift (as the time-shifts get away from zero, the numerator and denominator decrease in equivalent ways, given that we consider less and less of the information within the data window). So the cross-correlation with a large time-shift (that only measures the correlation of the end portion of the data with the start of the synthetic, or the start of the data with the end of the synthetic) contributes to $\Gamma(t)$ with the same weight as the cross-correlation with small time-shift (that measures the correlation of most of the data window with most of the synthetic window). One simple way of emphasizing the small time-shift cross-correlations is to normalize by the entire synthetic and data windows (constant normalization). This solution has the consequence of giving systematically lower CC values for larger time-shifts, which may lead to the rejection of some otherwise OK windows.
}
-Please think about this. I am getting more convinced that this is the right way to go, but it is admittedly tricky, and I'm glad you are discussing this with me!
+I see the idea, but from the standpoint of convergence in the inversion, it is better to use a ``close to right'' measurement than to use no measurement at all. Believe me, I don't like the idea of it either, but my objective is to enhance convergence. I am actually interested in up-weighting the huge time-shifts, if they meet all the fitting criteria. Widening the allowable time shift is not an option, because $95\%$ of the paths do not need it. In essence, the problem is that some paths are extremely well fit, while others are not at all. So, yes, I could tune the code per event, but that's not the point.
-{\em AM:
-As far as I can make out, there is no easy answer to the question of the best way to evaluate the correlation between two waveforms. We have to understand the consequences of the choices we make, and make those that make the most sense given the problem at hand.
- }
+At some point, you said you liked what was implemented in \verb+measure_adj+ (code below). But it does what you strongly disagree with, right? I mean, it will not allow any time shift outside the allowable limit, and it will keep the time shift, even if it is a non-local maximum. (In practice, I put this in as a safety, because there should be no measurements larger than what was indicated in FLEXWIN.)
+Finally, you have made the choice to look no more than exactly half the window length outside $[t_S, t_E]$, as indicated by \verb+fout = 0.5+. So you believe it should be 0.5, and not 1.0? Probably using 0.5 avoids some of the cycle-skipping that can happen, which is the same reason I am advocating tying \verb+fout+ to \verb+DT_BASE+ in some way.
+
+Basically, I'm okay with whatever we choose, but it would be nice to have agreement between the implementations in \verb+flexwin+ and \verb+measure_adj+.
+
+%-----------------------------------------
+
+\pagebreak
+\subsection*{CC formula in FLEXWIN}
+
+\begin{verbatim}
+ ! 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.
+ ! If fout=0.5, then it looks outside by a time of 0.5*(window width).
+ ! Perhaps fout should be a parameter, or it should be tied to the max
+ ! allowable time-shift specified by the user.
+ ! However, it does not matter as much if the data and synthetics are
+ ! zeroed outside the windows, as currently done in calc_criteria.
+ fout = 0.5
+ i_left = -1*int(fout*nlen)
+ i_right = int(fout*nlen)
+
+ ! 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.
+ 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 <== ORIGINAL WINDOW
+ ! 5 25 46 25 46
+ ! 10 31 52 31 52
+\end{verbatim}
+
+%-----------------------------------------
+
+\subsection*{CC formula in measurement code}
+
+\begin{verbatim}
+ ! these choices will slide the entire windowed record past the other
+ cr_shift = nlen*dt
+ n_left = ceiling( -1.0 * cr_shift / dt )
+ n_right = floor( cr_shift / dt )
+
+ ! cross-correlation
+ ishift = 0
+ do i = n_left, n_right, 1
+
+ cc = 0.
+ do j = 1, nlen
+ if((j+i) > 1 .and. (j+i) < nlen) cc = cc + dzr(j) * dzr2(j+i)
+ enddo
+
+ !if(cc > cc_max) then
+ ! CHT, 07-Sept-2008: Do not allow time shifts larger than the specified input
+ if(cc > cc_max .and. abs(i*dt) <= BEFORE_TSHIFT ) then
+ cc_max = cc
+ ishift = i
+ endif
+
+ enddo
+ tshift_xc = ishift*dt
+\end{verbatim}
+
%----------------------------------
\bibliographystyle{plainnat}
More information about the CIG-COMMITS
mailing list