Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Statistical analyses of ambient seismic noise spectra with applications to detecting imperfectly diffuse wave field and deriving attenuation and phase velocity information
(USC Thesis Other)
Statistical analyses of ambient seismic noise spectra with applications to detecting imperfectly diffuse wave field and deriving attenuation and phase velocity information
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Statistical analyses of ambient seismic noise spectra with applications to detecting
imperfectly diffuse wave field and deriving attenuation and phase velocity information
Xin Liu ( 刘昕)
A Dissertation
Presented to the Faculty of the University of Southern California
Graduate School
In partial fulfillment of requirements for the degree of
Doctor of Philosophy
in
Geological Sciences
August 2016
2
Acknowledgments
I really enjoy my graduate study here in the Department of Earth Sciences, USC. The past 5
years have been phenomenal experiences to me and my future career. I’ve made some interesting
new discoveries in the tiny field of ambient seismic/acoustic noise studies and I’m hoping to
develop and improve new methods on better understanding of the physical properties of the
shallow subsurface earth structure as well as its temporal change that can be related to tremors,
water content, temperature and other physical processes.
The first person that I’d like to thank is my advisor Prof Yehuda Ben-Zion. Yehuda actually
visited my home town Hohhot in Inner Mongolia of China back in 2000, when I was a kid. So
curiously 5 years ago in 2011 he chose me as one of his new students. Yehuda is a very serious
scientist and his expertise span across many different sub-fields in seismology, fault zone studies
and fracture mechanics. He has a lot of new ideas and can quickly grasp the key points of new
research directions. He encourages students to work on new methodologies/theories and is very
strict on the accuracy/correctness of any scientific research he supervises. He educated me to use
most precise and clear language in research papers to better convey my ideas to the readers.
Yehuda also helped me making connections with the top experts in my field by taking me to
academic meetings around the world, asking external scholars to cooperate and inviting people to
give a talk here. He bought powerful clusters and laptops for students so that we don’t need to
spend our limited stipend on research computing devices. Apart from academic training, he also
cares about students’ life in graduate school. He routinely invites his group to lunch/dinner at
home or restaurant. He also likes Erhu (a fiddle-like musical instrument from my hometown)
music and sometimes invited me to play Erhu at group meetings or his home. Three years ago I
had a very difficult time due to family emergency and I couldn’t come back on time when the
semester began. He offered a lot of help and support during that time and he never blamed me
for missing classes and falling behind in research progress. Finally I recovered from that crisis
and published our first paper. That was the closest point when I considered dropout my PhD and
I am very grateful to be his student.
I would also like to thank Prof Tom Jordan for his excellent lectures in Advanced
Seismology class and his help on explaining clearly attenuation, scattering, anisotropy and other
physical processes behind seismic wave propagation and suggesting key references in these
topics. That was very helpful when I just started on deriving the analytical form of ambient
seismic noise cross-correlation in anelastic medium.
3
I want to thank Prof Robert Scholtz for teaching me the beautiful theory of random processes
in electrical engineering so that I can later apply it to analyzing ambient seismic noise in novel
perspectives. I enjoyed the pure curiosity-driven discussions with him and the valuable insights
and advices he provided from signal communication perspective.
I’d like to thank Prof Charles Sammis and Prof Stephan Haas for serving on my committee.
They both asked me a lot of good questions from different perspectives that are very
enlightening.
I want to thank Dimitri Zigone, who is my collaborator and helped me through a lot of
difficulties in the early years of my graduate study. I’m also grateful to have all the friends here
in the Earth Sciences department. I really appreciate the help from the staffs in my department
especially Cindy Waite, who is very kind to students and understands their concerns and worries.
Finally, I would like to thank my parents, S.L. Guo and Y. H. Liu, and my sister Sheng Liu,
who are always my source of support and love.
4
Table of contents
Abstract ...................................................................................................................................................................... 6
Introduction .............................................................................................................................................................. 7
1. Estimating correlations of neighboring frequencies in ambient seismic noise .................................. 10
1.1. Introduction ...................................................................................................................................... 10
1.2. Theory .............................................................................................................................................. 11
1.3. Data processing and results .............................................................................................................. 13
1.3.1 Correlation matrix estimation from stationary synthetic data using different window length ... 15
1.3.2 Correlations among neighboring frequencies between 0.03-0.3 Hz with T=500 s ......................... 16
1.4. Synthetic simulations of correlation features ................................................................................... 19
1.4.1 Single cross-frequency random component with square/Gaussian shape .......................................... 21
1.4.2 Multiple cross-frequency random components (vectors) ........................................................................... 23
1.5. Data analysis using two additional frequency bands ....................................................................... 25
1.5.1 Correlations with T=100 s and higher frequency range ............................................................................ 26
1.5.2 Correlations with T=1000 s and lower frequency range ........................................................................... 28
1.6. Discussion ........................................................................................................................................ 29
2. Frequency Domain Analysis of Errors in Cross-Correlations of Ambient Seismic Noise .............. 33
2.1. Introduction ...................................................................................................................................... 33
2.2. Theory .............................................................................................................................................. 35
2.2.1. Sample mean cross-spectrum and its mean and variance ........................................................................ 36
2.2.2. Errors of amplitude and phase measurements .............................................................................................. 37
2.2.3. Signal-to-Noise Ratio on narrow-band time domain cross-correlation ............................................. 40
2.2.4. Stacking of piecewise stationary cross-spectrum ......................................................................................... 41
2.3. Numerical simulation ....................................................................................................................... 43
2.4. Application to regional network stations ......................................................................................... 49
2.5. Application to densely spaced stations with high frequency data .................................................... 58
2.6. Discussion ........................................................................................................................................ 62
Appendix 2A: Derivation of variance and pseudo-variance of sample mean cross-spectrum ................ 65
Appendix 2B: Derivations for errors of amplitude and phase from the 2nd moment statistics of the
sample mean cross-spectrum .................................................................................................................. 67
Appendix 2C: Correlation between noise and its Hilbert transform in frequency domain ..................... 69
Appendix 2D: Estimating noise and SNR on narrow-band filtered cross-correlation in time domain ... 70
3. Theoretical and Numerical Results on Effects of Attenuation on Correlation Functions of
Ambient Seismic Noise ........................................................................................................................................ 72
3.1. Introduction ...................................................................................................................................... 72
3.2. Theoretical results for isotropic source distribution ......................................................................... 74
3.2.1. Uniform attenuation .................................................................................................................................................. 76
3.2.2. Attenuation only within the inter-station path
out in
Q Q Q , ......................................................... 79
3.2.3. General Case ............................................................................................................................................................... 80
3.2.4. Dispersion and frequency-dependent attenuation ........................................................................................ 81
3.3. Numerical simulation results ........................................................................................................... 83
3.3.1. Isotropic source distribution with and without attenuation ..................................................................... 85
5
3.3.2. Asymmetric amplitudes from non-isotropic source distribution and attenuation ........................... 86
3.3.3. Isotropic noise distribution and attenuation only within the inter-station region .......................... 87
3.4. Discussion ........................................................................................................................................ 89
Appendix 3A: Inverse Fourier transform of Eq. 9: Uniform Attenuation. ............................................. 92
Appendix 3B: Inverse Fourier transform of Eq. 14: Attenuation in the Inter-station Region. ............... 96
Appendix 3C: Inverse Fourier transform of Eq. 22: dispersion and frequency-dependent attenuation 100
Appendix 3D: Derivation of normalized cross spectra (coherency) for uniform attenuation and isotropic
source distribution ................................................................................................................................. 102
4. Extracting Seismic Attenuation Coefficients from Cross-Correlations of Ambient Noise at Linear
Triplets of Stations ............................................................................................................................................. 104
4.1. Introduction .................................................................................................................................... 104
4.2. Theory ............................................................................................................................................ 106
4.3. Numerical simulations of noise waveforms and inversion of Q values ......................................... 109
4.4. Inverting Q values for a regional array .......................................................................................... 115
4.5. Inverting Q values for a fault zone array ....................................................................................... 119
4.6. Discussion ...................................................................................................................................... 125
Appendix 4A: Derivation of cross-spectrum for isotropic noise sources and asymmetric attenuation
structure ................................................................................................................................................ 129
Discussion ............................................................................................................................................................. 132
References ............................................................................................................................................................. 134
6
Abstract
I study statistical properties of ambient seismic noise wave field and possible errors and
biases in the empirical Green’s function obtained from ambient seismic noise cross-correlations.
The techniques I developed can be applied to understanding the spatio-temporal non-diffuse
characteristics of ambient seismic noise field and estimating errors of phase and amplitude
information of empirical Green’s functions (noise cross-correlations) computed from finite
amount of noise data. In addition, I derive the theoretical form of noise cross-correlation
functions for different inter-station/background attenuation Q values and develop a methodology
for estimating inter-station attenuation Q values from linear array of three stations (triplet).
A technique is developed for estimating the non-diffuse characteristics of ambient seismic
noise by computing correlations of neighboring power spectral values at different frequencies.
This technique is applied to broadband stations in southern California and show distinct features
for primary and secondary ocean microseismic peaks. The observed features can be explained by
adding cross-frequency random components to the fully diffuse noise through numerical
simulation.
A method is proposed for estimating the error of stacked noise cross-spectrum from evenly
spaced noise windows and deriving phase and amplitude errors for tomography purposes. The
method is verified through numerical simulations and then applied to southern California
regional network stations across San Andreas fault and dense linear array stations across the San
Jacinto fault zone with different frequency bands. Outliers of the ambient seismic noise are
removed using a statistical approach based on probability distributions of noise cross-spectrum to
improve the measurements of phase and amplitude from the stacked cross-spectrum.
A theory on the effects of inter-station and background attenuation structure on the expected
noise cross-correlation amplitude is constructed to improve the understanding on the anelastic
decay of the noise amplitude cross-spectra. Based on the theory, a methodology is proposed to
estimated inter-station attenuation from linear array of three stations, which is validated by
numerical simulations and applied to array data in regional network and dense linear array in
southern California. The attenuation Q
-1
values with confidence intervals are estimated through
bootstrap resampling by days and the results for a dense linear array suggest the location of a
fault damage zone.
7
Introduction
The ambient seismic noise wave field has been an emerging topic in seismology over the past
10 years. It consists of continuous background noise recorded by seismometers between episodic
impulsive (e.g. earthquake) or transient (tremor) signals. If the noise wave field is fully diffuse,
the empirical Green’s function that contains structural information between two locations can be
retrieved by cross-correlating continuous seismic noise recordings on two receivers. This
empirical Green’s function doesn’t depend on earthquakes and can achieve better Signal-Noise
Ratio (SNR) by stacking noise cross-correlation of longer time durations, which are advantages
over traditional earthquake based methods.
The first theoretical framework on deriving sub-surface information from earth’s stochastic
noise was made by Aki (1957), who found that the normalized cross-spectra for stationary
random noise approximate an zero.th order Bessel function with velocity information between
two locations. Claerbout (1968) suggested that for layered medium, the autocorrelation of noise
recordings produces reflection response of the subsurface for vertical incidence of plane waves.
Similar methods have been applied to study the acoustic impulse responses of the sun (e.g.
Rickett and Claerbout, 2000).
Recent studies have confirmed with detailed theories, numerical simulations, laboratory
experiments and field observations that the empirical Green’s function can be obtained from
ensemble averaged correlation of fully diffuse noise wave field at two different locations (e.g.
Lobkis and Weaver, 2001; Shapiro and Campillo 2004; Snieder, 2004; Weaver and Lobkis,
2004). In practice, however, only finite time duration of noise recording data are available.
Therefore the finite time averaged cross-correlation between two receivers contains errors related
to the length of noise recording. The time domain error is defined as the variance of noise on
cross-correlation of finite amount of ambient seismic noise data and it is proportional to the
product of autocorrelation functions at two receivers and inversely proportional to the length of
noise recording and bandwidth (e.g. Snieder, 2004; Sabra et al., 2005; Weaver and Lobkis, 2005).
The ambient seismic noise field is assumed fully diffuse for the successful retrieval of
empirical Green’s function. A fully diffuse noise wave field implies that the noise field at any
location is stationary and different wave modes corresponding to different frequencies are
uncorrelated. This condition may not hold in practice and therefore it is important to quantify the
degree of diffuse for ambient seismic noise wave field. The stabilized ratio of P and S wave
energy in the ambient noise can be an indicator for the diffuse property of wave field (e.g.
8
Weaver, 1982; Hennino et al., 2001; Margerin et al., 2009). The spectral ratio of
horizontal/vertical components seismic noise data is also related to the P/S wave energy ratio and
can be applied for similar purpose (Kawase et al., 2015).
Numerous works have computed tomography maps based on ambient seismic noise cross-
correlations. Campillo and Paul (2003) first computed empirical Green’s function by stacking the
diffuse earthquake coda wave cross-correlations. Then this technique was applied to ambient
seismic noise recorded on broadband stations (Shapiro and Campillo, 2004; Shapiro et al., 2005)
for travel time tomography (e.g. Shapiro et al., 2005; Yao et al., 2006; Lin et al., 2008; Boué et
al., 2014; Zigone et al., 2015).
In addition to travel time tomography, several studies have focused on extracting information
from the amplitude of noise cross-correlation. Prieto et al (2009) conjectured a method to derive
inter-station attenuation by fitting the normalized cross-spectra curve (coherency) to a damped
Bessel function. Weaver (2013) derived attenuation coefficients and site factors from time
domain narrowband cross-correlation amplitudes on linear array of at least four stations by
providing a theoretical framework with numerical verifications based on radiative transfer theory.
Lin et al (2012) estimated site amplification and attenuation from USArray by taking the gradient
and Laplacian of travel time and amplitude measurements of surface waves in earthquake data.
This method accounts for focusing/defocusing effect due to heterogeneous velocity structure but
it doesn’t separate attenuation from site factor because it cannot solve for two unknowns from
one real-valued equation. Applying this technique to amplitude measurements on ambient noise
cross-correlations require further assumptions about the empirical Green’s functions.
This dissertation consists of four chapters. The 1st chapter is a study on non-diffuse
properties of the ambient noise wave field through correlations of power spectral values at
neighboring frequencies. A possible explanation for the features of correlated neighboring
frequencies observed on real data is proposed with simulations based on cross-frequency
components. In the 2nd chapter, I propose a new methodology to estimate the error of stacked
noise cross-spectrum from finite amount of data and derive phase and amplitude errors which
can be eventually linked to phase velocity error and attenuation Q value error. A statistical
preprocessing method is also introduced to better preserve the amplitude information with
improved SNR and to replace the traditional prewhiteing and onebit noise normalization methods.
In chapter 3, a theory for expected noise cross-correlation amplitude due to different inter-station
and background attenuation Q values is derived for isotropic far field noise source distribution.
9
The possibility of varying asymmetry of cross-correlation function amplitudes with frequency
due to asymmetric attenuation structure is predicted. In chapter 4, I propose a methodology for
estimating interstation attenuation Q values based on three stations on a line. By taking the
amplitude cross-spectral ratio of two pairs of stations, the inter-station Q values can be obtained
by fitting the amplitude ratio within a frequency band. Attenuation Q values with bootstrap
confidence intervals are estimated from both regional array scale (inter-station spacing ~25 km)
and dense linear array (inter-station spacing ~20 m) across a fault zone.
10
1. Estimating correlations of neighboring frequencies in ambient seismic noise
(Liu and Ben-Zion, 2016)
1.1. Introduction
Cross-correlations of the ambient seismic noise field can be used to derive the Green's
function between pairs of recording stations (e.g., Lobkis and Weaver 2001; Snieder 2004;
Wapenaar 2004). This allowed many studies to obtain information from the ambient noise field
on seismic velocities at local, regional and global scales (e.g., Sabra et al. 2005; Shapiro et al.
2005; Lin et al. 2008;Boué et al. 2014; Zigone et al. 2015). Some studies also attempted to
retrieve attenuation coefficients from amplitude information of noise cross-correlations (e.g.
Prieto 2009; Weaver 2013; Liu et al. 2015). Accurate retrieval of the Green's function and related
wave properties (arrival times and attenuations) from noise cross correlations requires the
wavefield to be fully diffuse (e.g., Weaver and Lobkis 2004; Campillo 2006). It is therefore
important to develop techniques to estimate the degree to which the ambient noise field is diffuse.
A fully diffuse wavefield is stationary and wave components with different frequencies are
not correlated (e.g., Weaver 1982; Sánchez-Sesma et al. 2008; Snieder et al. 2009). In addition,
the wavefield is equipartitioned among all possible modes, implying that waves propagate with
equal energy in all directions and that the ratio of P to S wave energy has a certain limit value
that depends on the velocity structure (e.g., Weaver 1982; Hennino et al. 2001; Margerin et al.
2009)). These properties can be used to assess the regime of a wavefield but they are not easy to
estimate. Separating the P and S waves requires, for example, calculating the divergence and curl
of three-component seismograms recorded by arrays of sensors (e.g. Margerin et al. 2009).
Testing that waves propagate with equal energy in all directions may be done with beamforming
analysis (Sens-Schonfelder et al. 2015), which again requires array data.The spectral ratio of
horizontal-to-vertical seismic data is related to the ratio of P to S wave energy and can also be
used to estimate if a wavefield is diffuse (e.g. Kawase et al. 2015), but this involves additional
simplifying assumptions.
In the present paper we propose a simple method for analyzing deviations from fully diffuse
seismic noise from correlations between different frequency components of the ambient noise
field. The technique is based on calculating a matrix of correlation coefficients of power spectra
at different frequencies. This can be easily computed from evenly spaced time windows of single
11
component data recorded by individual stations. Observed correlations between different
frequency components indicate that the wavefield is not fully diffuse and provide additional
information on the noise sources and propagation regime at that frequency range.
In the following sections we first provide basic theoretical results on the method and data
processing. We evaluate the resolution and other properties of the technique by applying it to
different types (fully diffuse and partially correlated) of synthetic wavefields. The method is then
applied to data recorded by three broadband stations of the southern California seismic networks.
The observations at all three stations show clear correlations of neighboring frequencies of the
secondary microseismic peak around 0.15 Hz, as well as frequencies below the primary peak
around 0.06 Hz. The power spectra recorded by a station near the edge of Los Angles basin is
higher compared to data recorded by stations outside the basin. This holds both in an absolute
sense and the level of frequencies higher than 0.2 Hz relative to the level of the primary
microseism peak. As discussed in the final section of the paper, these results have important
information on the physical origins of the ambient seismic noise in different frequency ranges
and propagation properties of the data.
1.2. Theory
A diffuse wave field can be expressed analytically with normal mode expansion as (Weaver
and Lobkis 2004),
n
n
n
n
t i u a t exp ,
) (
r r (1)
where u
(n)
(r) are normal mode eigenfunctions, a
n
are random complex modal amplitudes, r is
location vector, t is time and ω
n
is angular frequency corresponding to the n.th mode. In a fully
diffuse wavefield different modal amplitudes are not correlated. Under this assumptions, the
mean value, pseudo-covariance and covariance of the modal amplitudes are (Weaver and Lobkis
2004), respectively
nm n m n
m n
n
F a a
a a
a
*
E
0 E
0 E
(2)
12
where E[] denotes expectation, δ
nm
is the Kronecker delta function and superscript * denotes
complex conjugate. The Kronecker delta in the covariance function indicates that the modal
amplitudes at neighboring frequencies are delta-correlated.
In this paper we adopt first two expressions of the equations set (2), but assume that the
modal amplitudes at neighboring frequencies can be correlated. Specifically, we assume that the
covariance of modal amplitudes is given by
(3)
where the function g(ω
n
, ω
m
) is generally complex and the modal amplitudes are assumed
complex Gaussian random variables. Taking the covariance between the absolute squares of the
complex modal amplitudes we get,
(4)
The absolute square of the covariance function in (4) is related to the correlation matrix of power
spectral values studied in the following data processing section.
The correlation of the random wave field at two different locations and different times (field-
field correlation) can be derived base on Eqs. (1) and (3) in a similar fashion as in Weaver and
Lobkis (2004). The result is,
(5)
where Re denotes the real part. Eq. (5) is different from Eq. (6) of Weaver and Lobkis (2004) by
depending not only on the time difference t-t’, but on both t and t’ for non-zero covariance of
random modal amplitudes between neighboring frequencies. In other words, the random noise
field with correlated neighboring frequencies is not necessarily stationary. However, if the
covariance function g(ω
n
, ω
m
) is close to a constant times the Kronecker delta δ
nm
, the partial
13
time derivative of the correlation with respect to t-t’ can still approximate the Green’s function
between r and r’.
The above derivations are based on "global" normal modes for the entire system. For local
random noise field in an open system (e.g. Snieder 2004; Weaver and Lobkis 2004), if the
neighboring frequencies are not delta-correlated or the noise data are not stationary, the field is
not fully diffuse. Therefore, in the following sections we will not compute eigenmodes, but
instead correlations among power spectral samples at different frequencies.
1.3. Data processing and results
Before data processing, it is important to address limitations of continuous noise recordings
at seismic stations. Observed seismograms are recorded with a finite sampling frequency over a
finite time interval and data analyses are performed on windowed time segments. Assuming a
rectangle window of length T, the truncated (windowed) noise recording in frequency domain is,
. (6)
Power spectra are used for correlations between neighboring frequencies to avoid computing
complex correlation coefficients. For actual noise data we use ordinary frequency f instead of
angular frequency ω. Because the autocorrelation of the windowed noise has a length of 2T, the
frequency interval is df=1/2T. Based on Eq. (6), each spectrum sample contributes to the main
and side lobes of the sinc functions around neighboring frequencies. Therefore, the frequency
resolution is determined by two neighboring frequency samples separated at the first zero-
crossings of each other’s sinc functions with spacing df
R
=1/T, such that the main lobes of two
sinc functions of the two samples do not overlap. The sampling rate of the noise data should be
at least twice as the highest frequency analyzed.
There are two limit situations when comparing the frequency resolution df
R
with the interval
between neighboring eigenfrequencies df
n
=(ω
n
-ω
n-1
)/2π. If the frequency resolution is much finer
than the eigenfrequency interval (df
R
<< df
n
), different random modal amplitude terms can be
resolved and separated. In contrast, if the frequency resolution interval is much larger than the
eigenfrequency interval (df
R
>> df
n
), each frequency sample of the random wavefield data
14
contains contributions from several nearby modes, due to the sum of the sinc functions in Eq. (6),
and the random wave field spectra can be treated as continuous. The former case is relevant only
for global normal modes when the window length T is sufficiently long. The latter case is more
common and adopted in this study because of the frequency ranges of interest and the limitation
of employed window lengths.
The covariance of neighboring random modal amplitudes in Eq. (3) cannot be directly
derived from data. However, the correlation coefficient between random wave field power
spectral samples at different frequencies can be derived from the data using,
(7)
where ω
p
and ω
q
are two different frequency samples in the spectra of the truncated
(windowed) data ψ
tr
. If the frequency resolution is much finer than the eigenfrequency interval
(df
R
<< df
n
), Eq. (7) can be related to Eq. (4) when the covariance in Eq. (4) is normalized to the
correlation coefficient between two modal amplitudes. If the resolution is poorer than the
eigenfrequency interval (df
R
>> df
n
), the correlation coefficient in Eq. (7) reflects the correlation
between two groups of modal amplitudes, with each group centered at a different frequency.
Using correlation between power spectrum samples has the advantage that time shift does not
have any effect on the result of correlation coefficient:
(8)
where τ is any arbitrary time shift applied to the windowed noise.
If the power spectrum of noise data follows a certain probability distribution, the power
spectrum distribution at each frequency can be computed from numerous non-overlapping
windows that divide that data into segments of length T s. Assuming that each window represents
15
an observation of the random noise spectra, the correlation coefficient of two power spectral
values at different frequencies can be conveniently computed based on Eq. (7).
1.3.1 Correlation matrix estimation from stationary synthetic data using different window
length
To obtain basic guidelines for estimating correlation coefficient among neighboring
frequencies based on Eqs. (6) and (7), we perform synthetic tests on stationary noise data. We
first generate fully diffuse stationary synthetic noise data using the noise model of Patterson
(1993) for a total duration of one month. Then experiments are done using two different window
lengths, 500 s and 100 s, from the synthetic noise data. For each window lengthT, evenly spaced
windows with gaps of 10%T are selected from the continuous synthetic noise. The power
spectrum for each window is computed using the Fast Fourier Transform (FFT). Since the
autocorrelation of each windowed noise segment has a length of 2T (to avoid cyclicity), the
windowed segments are zero-padded to a length of 2T before using the FFT. Then power
spectrum values from many windows are combined to estimate correlation coefficient matrix
based on Eq. (7).
Figure 1.Correlation coefficient matrices of power spectra computed from one-month synthetic
stationary diffuse noise using (a) 500 s and (b) 100 s time windows. The correlation matrix
results are diagonal within the resolution limits. The stationary noise is simulated using the
model of Patterson (1993).
Fig. 1(a) and. 1(b) show the correlation coefficient matrices of 500 s and 100 s window
lengths estimated from the stationary synthetic data. The diagonal elements have correlation
16
coefficient values of one. They are surrounded by regions, with width inversely proportional to
the window length, having appreciable coefficient values due to the resolution limit of window
length discussed below Eq. (6). All the other elements show correlation values very close to zero
(below 0.03). These two examples suggest that the correlation matrix should be diagonal at its
resolution limit if the ambient noise time series is fully diffuse and stationary.
1.3.2 Correlations among neighboring frequencies between 0.03-0.3 Hz with T=500 s
−1 19 ˚ −1 18 ˚ −1 17 ˚
34˚
35˚
CHF
LMR2
OLI
Figure 2.Locations of three broadband stations of the southern California seismic network used
in this work.
In this section we study the correlation of neighboring frequencies in data recorded by three
broadband stations of the southern California network (Fig. 2). These stations are selected
because they have high-quality continuous data and can be used to show both common features
as well as differences at different locations. Station CHF is located in the San Gabriel mountains,
station LMR2 is in the Mojave desert and station OLI is close to the ocean between the Los
Angeles basin and Palos Verdes. Continuous waveform data recorded by the vertical component
(BHZ) of these stations between days 50 and 300 of the year 2014 are analyzed. The waveforms
are first decimated to 4 Hz and then high-pass filtered at 0.01 Hz and clipped at 4 times the
17
standard deviation within each four-hour long segment to reduce effects of small earthquakes.
Segments with larger earthquakes are simply removed.
This data analysis is centered around the primary (0.06 Hz) and secondary (0.15 Hz) ocean
microseismic peaks in the frequency band 0.03-0.3 Hz. The window length T is set to 500 s and
the gap between two windows is 50 s. Autocorrelation of length 2T is computed in the frequency
domain for each noise window, and the power spectral values of all windows are stored for each
analyzed frequency. Outlier values for each frequency are identified using the distribution of
power spectra. Windows with more than 5% outliers in their frequency content are excluded
from further analysis to minimize effects of small earthquakes and other transient noise sources
on neighboring frequency correlations and the stacked power spectrum. The remaining windows
are used for subsequent analyses discussed below.
Figure 3.Power spectra over 0.03-0.3 Hz and neighboring frequency correlation matrices
obtained for days 50-300 of year 2014 with window length of 500 s and frequency resolution of
0.002 Hz. (a) Results for station CHF: (top panel) Stacked power spectra (red line) and number
of power spectrum measurements with a given value in each power bin (gray scale) of each
discrete frequency value. The primary microseismic peak is around 0.06 Hz and secondary
microseismic peak is near 0.15 Hz. (bottom panel) Correlation coefficient matrix of the power
spectra. Three groups of correlated frequency features are observed: square (0.07-0.09 Hz),
18
diagonal halo (0.1-0.2 Hz) and stripes of correlated frequencies. (b) Corresponding results for
station LMR2.
Fig. 3(a) presents the stacked power spectra (top) and correlation coefficient matrix of
neighboring frequency (bottom) for station CHF. The correlation coefficient matrix is computed
based on Eq. (6) using power spectral values from numerous windows as done for the synthetic
test. The results show several features not present in the basic synthetic results of Fig. 1. The first
notable feature is that the secondary ocean microseismic peak around 0.15 Hz (top plot)
corresponds well with a significant correlated (~0.36) halo between 0.1-0.2 Hz in the correlation
coefficient matrix (bottom plot). The primary microseism peak near 0.06 Hz does not produce a
similar halo. However, the two troughs around the primary peak in the power spectrum are
associated with two correlated frequency bands at 0.03-0.05 Hz and 0.07-0.09 Hz that produce a
square-type shape in the bottom correlation plot. Additional notable features are the stripes in the
bottom plot associated with correlations over relatively wide frequency range (perhaps between
the primary and secondary ocean microseismic bands. Similar patterns are seen in Fig. 3(b) for
station LMR2 with slightly different spectral power and correlation coefficients.
19
Figure 4.Similar to Fig. (3) for station OLI with stacked power spectra and distribution of power
spectrum measurements at each frequency (top) and correlation coefficient matrix (bottom). The
correlation matrix for the data at OLI shows additional correlated features above 0.25 Hz
compared to Fig. (3).
The results for station OLI contain these and other features (Fig. 4). The power spectrum for
OLI (top plot) is amplified compared to those of CHF and LMR2, especially for the range above
0.15 Hz. The correlation matrix of OLI indicates another frequency band of noticeable
correlation (~0.3) among neighboring frequencies above 0.25 Hz (bottom plot). These
differences are seen more clearly in results for higher frequencies presented later.
1.4. Synthetic simulations of correlation features
To clarify the origin of the features seen in the correlation plots of the observed data (Figs. 3-
4), we perform additional synthetic tests. For N discrete frequency values, we simulate the noise
spectra vector d
N×1
of correlated frequency samples from fully uncorrelated random vector
m
(N+1)×1
of zero-mean standard complex Gaussian random variables,
20
(9a)
or equivalently
1 1 1 1 1 1 1
N N N N N N N N N N N N
m I s C m G C d , (9b)
where the real vector s
N×1
represents the cross-frequency component contributing to each
frequency sample, I
N×N
is an identity matrix, and C
N×N
rescales the variance of each resulting
random spectrum sample d(ω
i
) to unity. The matrix G
N×(N+1)
transforms an uncorrelated random
vector to a correlated random vector. More importantly, it adds cross-frequency random
component to the uncorrelated standard complex Gaussian random vector (corresponding to a
fully diffuse wave field) with predictable statistics and correlation of the resulting wave field.
The correlation matrix can be computed from ensemble average of product of random spectra
vectors,
N N N N
T
N N N N
T
N N N N N N N N
T
N N N N
C I s s C G C G C d d R
1 1 1 1
*
1 1
E (10)
which shows as expected that the diagonal entries of the matrix R
N×N
are equal to 1. The
correlation matrix of observed data is computed from correlations of power spectral samples
based on Eq. (7). Therefore each element of the power spectra correlation matrix is actually the
absolute square of the corresponding element in R
N×N
.
In the next two sub-sections, we first define one or multiple cross-frequency component
vectors and then compute the corresponding correlated noise spectra from uncorrelated random
vectors. Using the synthetic results, we produce correlation matrix plots with features of similar
shapes (square, halo, stripes) seen in the correlation results of the observed data.
21
1.4.1 Single cross-frequency random component with square/Gaussian shape
Figure 5.Synthetic results for non-stationary noise with a simple cross-frequency component
vector. (a) Cross-frequency component with values of 1.22 between 0.07-0.09 Hz and zero
elsewhere. (b) Correlation coefficient matrix of power spectra computed from 2000 realizations
of random noise vectors simulated usingthe boxcar cross-frequency component vector in (a). (c)
A realization of simulated noise using the cross-frequency component vector in (a).
22
Figure 6.Similar to Fig. (5) with a cross-frequency component vector shown in (a) consisting of
a boxcar between 0.07-0.09 Hz and a Gaussian function centered at 0.15 Hz elsewhere. The
maximum amplitudes are 1.22. (b) Correlation coefficient matrix of power spectra computed
from 2000 realizations of random noise vectors simulated using the cross-frequency component
vector in (a). The Gaussian and boxcar are correlated on the correlation matrix because they are
both affected by the same random variable m
0
in eq. (9a). (c) A realization of simulated noise
using the cross-frequency component vector in (a).
If there is only one cross-frequency component vector s
N×1
, the noise data simulation takes
exactly the same form as Eq. (9). Different forms of the vector s
N×1
produce different correlation
matrix features. As a first example, we consider a case with a cross-frequency component vector
defined by a simple boxcar function that equals 1.22 between 0.07-0.09 Hz and zero elsewhere
(Fig. 5a). Straightforward calculation based on Eqs. (9) and (10) show that all frequencies within
0.07-0.09 Hz are correlated with each other with a coefficient of 0.6 (except the diagonal and
surrounding small neighborhood that have higher values). The correlation coefficient matrix of
power spectra computed from 2000 random realizations has off-diagonal peak value of 0.36 (the
square of 0.6) in a square bounded between 0.07-0.09 Hz (Fig. 5b). A simulated noise of a single
random realization has significant non-stationary features around t = 0 s (Fig. 5c). The random
23
noise can be time-shifted arbitrarily and still preserve the same correlation matrix of power
spectrum samples according to Eq. (8).
In the second example the cross-frequency component vector (Fig. 6a) consists of the boxcar
function as before augmented by a Gaussian function (same peak value 1.22 centered at 0.15 Hz
and standard deviation is 0.17 Hz). Because the boxcar and the Gaussian function are on the
same cross-component vector with the same peak value, they are correlated with each other with
peak correlation coefficient of 0.6. The correlation matrix of power spectrum samples (Fig. 6b)
has three different features with peak value of 0.36: (1) the square that is the same as in previous
example; (2) a 2D Gaussian centered at 0.15 Hz and (3) vertical and horizontal stripes produced
by the correlation between the square and the Gaussian functions. A simulated noise of one
realization based on this cross-frequency component vector has again a non-stationary feature
(anomalous high amplitude) around t = 0 s (Fig. 6c).
1.4.2 Multiple cross-frequency random components (vectors)
The analysis of observed data in section 1.3.2 yields (Fig. 3 & 4), in addition to box and
stripe patterns, extended halo around the diagonal. This suggests that multiple cross-frequency
random component vectors may affect the data. To simulate random noise with multiple cross-
frequency components and corresponding power spectra correlation matrix, Eq. (9) can be
generalized by replacing the vector s
N×1
with a matrix s
N×K
, where K is the number of
independent cross-frequency component vectors. The length of uncorrelated complex Gaussian
random vector m
(N+K)×1
is N+K. The diagonal elements of the rescaling matrix C
N×N
in Eq. (9)
also needs to be modified to include squares of all the cross-frequency components at same
frequency.
24
Figure 7.Synthetic results for non-stationary noise with 20 cross-frequency component vectors.
(a) Assumed 20 cross-frequency component vectors, each with a boxcar between 0.07-0.09 Hz
and a Gaussian functionshifted at different center frequency with a variable variance. The
maximum amplitudes are indicated by the color bar (0.25 for boxcar and 0.82 for Gaussian
function). (b) Correlation coefficient matrix of power spectra computed from 2000 realizations
of random noise vectors simulated using the cross-frequency component vector in (a).The
Gaussian and boxcar are correlated on the correlation matrix because they are both affected by
the same random variable m
0
in eq. (9a). (c) A realization of the simulated noise using the cross-
frequency component vector in (a).
As an example, we assume 20 independent cross-frequency component vectors (Fig. 7a),
each with a boxcar equal to 0.25 between 0.07-0.09 Hz and a Gaussian function with peak value
of 0.82, variable variance and shifting center between 0.1-0.2 Hz. The correlation matrix of
power spectrum samples (Fig. 7b) shows a square between 0.07-0.09 Hz with correlation
coefficient of 0.3, an elongated halo between 0.1-0.2 Hz with peak correlation of ~0.36 (square
of 0.6) and variable width from the diagonal, and stripes with peak correlation of 0.08 due to the
correlation between the boxcar and Gaussian functions of the multiple cross-frequency
components. A simulated noise of one realization based on the multiple cross-frequency
25
component vectors produces as in the previous examples non-stationary high amplitude around t
= 0 s (Fig. 7c).
The similarity of the features in Fig. 7b to those seen in the correlation plots of the observed
data in Figs. (3) and (4) suggests that the ocean microseismic noise wave field is associated with
multiple sources at somewhat different frequencies (or nonlinear scattering involving frequency
shifts), that produce correlated rather than fully diffuse field. A quantitative comparison between
the correlated neighboring frequency features at the synthetics and observed data is not presented
here but may be performed in a future work involving inversion of the cross-frequency
components.
1.5. Data analysis using two additional frequency bands
Figure 8.Power spectra over 0.05-0.6 Hz and neighboring frequency correlation matrix obtained
for days 50-300 of year 2014 with window length of 100 s and frequency resolution of 0.01 Hz.
(a) Stacked power spectra with distribution of power spectrum measurements at each frequency
(top) and correlation coefficient matrix (bottom) for station CHF. The three correlated frequency
features in Fig. (3) are visible but with lower resolution. (b) Corresponding results for station
LMR2 with similar correlated frequency features.
26
To explore further non-diffuse characteristics of data recorded by the same three stations of
the CI network (Fig. 2), we use two additional window lengths, 100 s and 1000 s. The 100 s
window length corresponds to a resolution limit of 0.01 Hz and it can be used to study
frequencies between 0.05-0.6 Hz. On the other hand, the 1000 s window length is suitable for
much lower frequencies between 0.002-0.1 Hz with a resolution limit of 0.001 Hz.
1.5.1 Correlations with T=100 s and higher frequency range
In this section we use a shorter time window to analyze slightly higher frequencies between
0.05-0.6 Hz, using the same preprocessed data as in section 1.3.2. The window length is set to
100 s and the gap between consecutive windows is 10 s. As in section 1.3.2, the power spectral
values of all windows are calculated and stored for each frequency within the used band. Outlier
windows are again removed to minimize effects of small earthquakes and episodic noise sources
on the results.
Figure 9.Similar to Fig. (8) for station for OLI. Stacked power spectra and distribution of power
spectrum measurements at each frequency (top) and correlation coefficient matrix (bottom).
27
There is an additional correlated feature between 0.25-0.6 Hz compared to Fig. (8), which may
reflects basin effects.
Figs. (8a) and (8b) show stacked power spectra (top) and correlation matrix plots (bottom)
for stations CHF and LMR2, with broader frequency range and lower frequency resolution
compared to the results of Figs. (3a) and (3b). The neighboring frequency correlations near the
secondary ocean microseism peak at 0.15 Hz (halos) and the correlated frequency band between
0.07-0.09 Hz (box) are similar to those shown earlier for the low frequency range. The power
spectrum at 0.08 Hz also correlates with the secondary ocean microseismic peak at 0.15 Hz
(stripes). The frequency content above 0.3 Hz not examined in section 1.3.2 shows a broad zone
of weak correlations (between 0.1 and 0.18).
Figure 10.Power spectra over 0.002-0.1 Hz and neighboring frequency correlation matrix
obtained for days 50-300 of year 2014 with window length of 1000 s and frequency resolution of
0.001 Hz. Stacked power spectra and distribution of power spectrum measurements at each
frequency (top) and correlation coefficient matrix (bottom) for station CHF.
28
Station OLI again exhibits higher power spectra than observed at stations CHF and LMR2
and additional features in the correlation matrix (Fig. 9). There is a broad frequency range (up to
the analyzed 0.6 Hz) with higher power than the primary microseismic peak around 0.06 Hz. The
correlation coefficient matrix shows clearly a band of correlated neighboring frequencies
between 0.25-0.6 Hz, in addition to the lower frequency observations that are seen better in Fig.
(4). This high frequency band corresponds to anomalously high power relative to the primary
microseismic peak compared with the power spectra of CHF and LMR2.
1.5.2 Correlations with T=1000 s and lower frequency range
Figure 11.Similar to Fig. (10) for station for OLI. The correlated frequency feature below 0.05
Hz is stronger and broader compared with the results for station CHF in Fig. (10).
Finally, we use a longer time window to analyze lower frequencies between 0.002-0.1 Hz.
Here the waveforms are high-pass filtered at 0.001 Hz and the same preprocessing steps as in
section 1.3.2 are applied to suppress or remove earthquake signals. The window length is set to
29
1000 s and the gap between windows is 100 s. As before, power spectral values are calculated
and outliers are removed from subsequent analysis.
Fig. 10 shows the stacked power spectra and correlation matrix for station CHF. The finer
frequency resolution produces evidently stronger correlations (~0.6) among neighboring
frequencies between 0.005-0.05 Hz compared to the lower resolution results of sections 1.3.2
and 1.5.1. The square type correlations between 0.07-0.09 Hz are also shown more clearly. The
primary microseismic peak at ~0.06 Hz again corresponds to a less correlated frequency range,
suggesting more diffuse wavefield compared to those associated with lower frequency range
0.005-0.05 Hz, and the higher frequencies associated with the secondary microseismic peak (Fig.
3a).
The correlation matrix of station OLI (Fig. 11) shows stronger correlations compared to the
results at station CHF. The correlated frequencies between 0.005-0.05 Hz have correlation
coefficient of ~0.7 and the shape of the correlated zone approaches a square. This range is also
more strongly correlated with another square at 0.07-0.09 Hz. The overall power spectral
amplitude is similar to that of station CHF. The primary microseismic peak at ~0.06 Hz is also
less correlated with its neighboring frequencies as in station CHF.
1.6. Discussion
We propose a method to detect deviations in given frequency ranges of the ambient seismic
noise field from a fully diffuse regime, using a correlation matrix of power spectral samples. The
method is easy to use and specific signatures in the correlation matrix provide information on the
forms of correlations between different frequency components of the ambient noise field. The
existence of correlation in a given frequency bandwidth may produce a bias in the empirical
Green’s function computed from noise cross-correlation in that band. Variations of correlation
matrix characteristics in different frequency ranges observed at different stations imply
differences in the noise source mechanisms and/or general propagation/scattering aspects in
these frequency ranges. Observed variations at different stations provide additional information
on propagation/scattering properties specific to different sites.
30
Weaver and Lobkis (2004) showed that diffuse wavefield implies delta-correlation between
different frequency components of (normal mode expansion of) the wavefield. With this
assumption, they arrived at the conclusion that the time derivative of field-field correlation is the
imaginary part of the Green’s function. In our derivation (section 1.2) we relax the delta-
correlation condition and assume a covariance function for neighboring modal amplitudes at
different frequencies. We then derive the correlation of random noise field at two different
locations in a similar way as Weaver and Lobkis (2004), but our result implies non-stationary
random noise due to cross terms in neighboring frequencies. The correlation of random field
between two locations produces limitations to the ability of approximating the Green’s function
from cross correlations of the noise field at these locations.
The correlation coefficient matrix of power spectral samples is estimated in this work from
evenly spaced, non-overlapping windows of same length. Simple rectangle window is used on
the noise time series and the frequency resolution of the corresponding random spectra is
determined from the first zero-crossing of the sinc function. To verify the underlying idea of the
method, we apply it to one-month fully diffuse and stationary synthetic noise data simulated
using the Patterson (1993) noise power spectral model. Synthetic results with 500 s and 100 s
window lengths show diagonal correlation matrices of power spectrum values within the
resolution limits. These results demonstrate what should be obtained if the delta-correlated
condition of different frequency components is satisfied and the wavefield is fully diffuse.
Applying the technique to data recorded by three stations of the southern California seismic
network, we estimate the correlations of neighboring frequencies using different time window
lengths (100-1000 s) to study the properties of the ambient noise field in different frequency
ranges. The correlation matrix results for the three stations show significant common deviations
from the diagonal pattern expected for a fully diffuse field, along with some differences in the
data recorded by the different stations. Results for analysis with a time window length of 500 s,
focusing on the frequency range 0.05-0.3 Hz, show elevated correlations in a square between
0.07-0.09 Hz, a halo around the diagonal in the range 0.1-0.2 Hz around the secondary
microseismic peak (0.15 Hz) and correlation stripes between the square and halo at all three
stations. To probe the possible causes of these types of correlated features, we simulate synthetic
noise data where cross-frequency random components are added to an uncorrelated (fully diffuse)
noise spectral vector. The cross-frequency components may represent correlation of noise
31
sources in given frequency ranges and a wavefield that is not sufficiently scattered during
propagation from sources to receivers.
The synthetic simulations indicate that the square and stripes can be reproduced with one
cross-frequency component, while the halo pattern requires multiple cross-frequency
components. The correlation coefficients are controlled by the ratio of variances of cross-
frequency components and the uncorrelated noise (Eqs. 9 & 10). Moreover, the simulated noise
for each realization is non-stationary and the peak amplitude at t = 0 s is higher than at different
times. This suggests that the correlated frequency features may significantly increase the noise
amplitude variations with time, since time-shift does not change the correlation of power spectral
values. Eq. (9) may be applied in a future study to form an inversion of the observed correlation
coefficient matrix to determine the number and form of independent cross-frequency
components needed to match the off-diagonal energy. Such results derived from data at many
stations can provide insights on the mechanisms (source, structural or combined) producing the
cross-frequency components.
Results with a shorter time window of 100 s for the frequency range 0.05-0.6 Hz show again
similar correlated neighboring frequency features at 0.1-0.2 Hz (around the secondary
microseismic peak) and 0.07-0.09 Hz on all three stations, with lower resolution than in the case
using 500 s windows. In addition, the correlation matrix of station OLI has significant
correlations of neighboring frequencies between 0.25-0.6 Hz. The frequency range at that station
above 0.2 Hz also has higher power than the primary microseismic peak. The elevated power and
correlations at this frequency band recorded by station OLI may be related to reverberations at
the edge of the Los Angeles basin or near-coast propagation effects. This and additional possible
effects of different geological settings and proximity to the coast should be tested in a future
study with data at many more stations.
Analysis of lower frequencies in the range 0.002-0.1 Hz based on a window length of 1000 s
for stations CHF and OLI shows stronger correlations between 0.005-0.05 Hz than the correlated
features at higher frequencies (0.07-0.09 Hz and 0.1-0.2 Hz). These results indicate that the
ambient noise at such very low frequency is much less diffuse (or less equipartitioned (e.g. Sens-
Schönfelder et al. 2015). The frequency band 0.005-0.05 Hz partly overlaps with the range of the
earth hum (e.g., Tanimoto 2005; Webb 2008;) related to global normal modes. In contrast, the
32
primary microseismic peak (~0.06 Hz) corresponds to a more diffuse noise field among
neighboring frequencies. These differences suggest that the primary and secondary microseismic
peaks may have some different source mechanisms and/or propagation properties.
Finally, we note that the present study does not take into account seasonal effects of the
ambient noise field, which has been observed on wide ranges of frequencies (e.g. Hillers and
Ben-Zion 2011; Landès et al. 2010; Traer et al. 2012), and the three used stations reveal features
of the ambient noise in a relatively small region. Future studies should perform similar analyses
over broader frequency range using more stations in different regions and data covering multiple
years.
33
2. Frequency Domain Analysis of Errors in Cross-Correlations of Ambient Seismic Noise
(Liu et al., 2016)
2.1. Introduction
Extracting empirical Green’s function between multiple pairs of stations from ambient
seismic noise cross-correlation has been widely applied in various regions and scales (e.g.
Shapiro and Campillo 2004; Sabra et al. 2005a; Bensen et al., 2007; Lin et al., 2008; Hillers et al.
2014). Aki (1957) proposed that seismic properties of the subsurface material can be obtained
from the ensemble-averaged cross-spectra of two different locations in a stochastic wave field.
An empirical Green’s function can be derived from the expected cross-correlation of fully
diffused stationary noise wavefield (e.g. Lobkis and Weaver 2001; Roux et al., 2005; Sánchez-
Sesma and Campillo 2006;Gouédard et al., 2008). In practice, the ambient noise field may not be
fully diffuse at certain frequency ranges (e.g. Sens-Schönfelder et al. 2015; Liu and Ben-Zion
2016) due to imperfectly scattered ocean microseismic noise or other noise sources. The non-
diffuse ambient noise wave field may produce biased empirical Green’s function estimated from
cross-correlation (Liu and Ben-Zion 2016). Furthermore, the distribution of ambient seismic
noise sources can be non-isotropic (e.g. Campillo 2006; Weaver 2009; Tsai 2011), which can
also bias the empirical Green’s function. Despite those complications, ambient seismic
correlations have been applied in numerous studies of phase/group velocity tomography (Shapiro
et al., 2005; Lin et al., 2008; Zigone et al., 2015), attenuation coefficient/site factor estimation
(Prieto 2009; Lin et al., 2012; Liu et al., 2015), and body wave phases retrieval (e.g. Poli et al.,
2012; Lin and Tsai, 2013). In addition, noise correlations have been used to analyze temporal
changes of seismic velocities (e.g. Sens-Schonfelder and Wegler 2006; Brenguier et al. 2008;
Brenguier et al. 2014; Hillers et al. 2015).
The early work on convergence rate of diffused wavefield cross-correlation primarily
focused on the variance of time domain noise cross-correlation. Snieder (2004) derived the time
domain variance expression for isotropically distributed scatters (sources) embedded in
homogenous media. The time domain cross-correlation variance is proportional to the product of
two autocorrelation functions from two stations at zero lag time divided by the time-bandwidth
product. Weaver and Lobkis (2005) derived variance expression that decays like 1/(recording
time) for open systems based on modal expansion of diffused wavefield. Sabra et al. (2005a)
34
generalized the variance expression to heterogeneous medium and colored noise spectra
assuming isotropic noise source distribution and stationary noise process. They estimated
variance by windowing the coda of cross-correlation where the empirical Green’s function
equals zero. Based on the definition of variance for cross-correlation, the Signal/Noise Ratio
(SNR) of narrowband filtered cross-correlation can be conveniently defined as the ratio of peak
amplitude of wave packet envelope over the square root of variance (e.g., Lin et al., 2008; Sabra
et al., 2005b; Poli et al., 2013; Zigone et al., 2015). The square root of variance of coda noise is
simply the Root Mean Square (RMS) of the same noise segment.
In the present paper we develop a probabilistic description of the sample mean (stacked)
cross-spectrum based on assumptions on Gaussian noise model, stationary random process and N
independent, identically distributed (i.i.d) random cross-spectrum samples. We find that the
variance of sample mean cross-spectrum is equal to product of power spectrum values of two
stations at the same frequency divided by number of samples (section 2.1), which is linked to the
time domain variance expression mentioned before. The phase error and causal (or anti-causal)
part amplitude error based on standard errors of complex sample mean cross-spectrum are
derived in section 2.2. The SNR for narrowband Gaussian filtered cross-correlation is estimated
in time domain and linked to the variance and mean of sample mean cross-spectrum (section 2.3).
The mean and variance estimations for piecewise stationary ambient noise records are given in
section 2.4. We simulate synthetic noise data at two stations and study the distribution and
statistics of sample mean cross-spectrum at each frequency computed from numerous non-
overlapping windows (samples). The modified mean cross-spectrum distributions and causal
amplitude spectral decay curves based on two widely used preprocessing techniques, one-bit
normalization and pre-whitening, are compared to those associated with unprocessed data.
Synthetic data are also used to study phase and amplitude errors as well as the convergence
images of SNR. We apply similar error analysis procedures to three regional stations (~35 km
spacing) of the Southern California regional network crossing the San Andreas Fault. The
probability distribution of real cross-spectrum data at each frequency contains outliers that
should be removed. After removing outliers, we find SNR convergence images containing peaks
and troughs at different frequencies. Block bootstrap method is applied to the same dataset as an
independent way of estimating the empirical confidence intervals of sample mean cross-
spectrum compared with the results based on standard error of mean cross-spectrum. The results
35
suggest that the temporal correlation of cross-spectrum samples among different windows
increases the variance of sample mean cross-spectrum (between 0.05-0.2 Hz). In addition, we
analyze the errors in high frequency noise recorded by densely spaced (~20 m) arrays that cross
the San Jacinto Fault Zone and show that by removing outliers the SNR is increased, which
improves both amplitude and phase velocity estimations.
2.2. Theory
We consider ambient noise wave field on two stations a and b. The random wave field d
j
(ω)
at each station j (j = a or b) can be expressed as the sum of correlated noise u
j
(ω), which stands
out through cross-correlation, and uncorrelated noise v
j
(ω), which cancels out after cross-
correlation,
j j j
v u d . Assuming the correlated noise field at different sites have
common noise sources, the correlated noise spectrum of station j at angular frequency ω is
written as a sum of contributions from numerous noise sources s
k
(ω) (discrete version of Eq. 1 in
Liu and Ben-Zion 2013),
(1)
where Δ
jk
represents distance from noise source k to station j. The inverse complex phase
velocity is defined as c Q i c ) ( 2 / sgn 1
~
1 where Q is the spatially varying
(heterogeneous) attenuation quality factor and c(ω) is the real phase velocity.
The noise sources s
k
(ω) are assumed to be circular Gaussian random variables that are
mutually uncorrelated such that
kl k l k
B s s
*
E , where B
k
(ω) is the spectral density of
k.th noise source and δ
kl
is the Kronecker delta. The noise sources distribution can be non-
isotropic. As a result of the summation in Eq. (1), the correlated noise spectrum u
j
(ω) is a
Gaussian random variable. The uncorrelated noise spectrum is also assumed circular Gaussian
random variable and satisfies 0 E
*
l k
v u and
kl l l k
V v v
*
E , where V
l
(ω) is the
spectral density of uncorrelated noise at station l.
36
2.2.1. Sample mean cross-spectrum and its mean and variance
In frequency domain the cross-correlation operation can be expressed as multiplication of the
spectrums on two stations. In realistic cases there are only finite length of noise recordings
available. Here we denote the n.th sample of random spectrum on station j as
n
j
d , assuming N
independent and identically distributed (i.i.d) spectrum samples at the same frequency. Then the
sample mean of N i.i.d cross-spectrum sample values from stations a and b can be formulated as,
N
n
n
b
n
a
N
ab
d d
N
R
1
*
1
. (2)
The expectation of the sample mean cross-spectrum is
ab
N
n
n
b
n
a
N
n
n
b
n
a
N
ab
C u u
N
d d
N
R
1
*
1
*
E
1
E
1
E (3)
where C
ab
(ω) is the expected cross-spectrum (ensemble average) of correlated noise recordings
from two stations a and b (e.g. Weaver et al., 2009; Liu and Ben-Zion, 2013; Liu et al., 2015).
The uncorrelated noise terms v
a
(ω) and v
b
(ω) cancel out.
The sample mean cross-spectrum is a complex random variable that approaches a
Gaussian distribution for a large sample size N according to the Central Limit Theorem of
probability. The variance of is (Appendix A),
(4)
where A
aa
(ω) and A
bb
(ω) are expected (ensemble-averaged) power spectrum functions for
stations a and b, respectively.
The pseudo-variance of the sample mean cross-spectrum is (Appendix A),
ab ab b a b a
N
ab
C C
N
d d d d
N
R
1
E E
1
pVar
* *
(5)
37
which is generally not zero. A non-zero pseudo-variance indicates that the complex random
variable is not circular and therefore its real and imaginary parts are correlated
(Appendix B).
2.2.2. Errors of amplitude and phase measurements
For N i.i.d cross-spectrum samples (from finite amount of noise data), the standard errors for
both real and imaginary parts of the sample mean cross-spectrum can be measured directly from
the distributions of N cross-spectrum samples. The sample mean cross-spectrum can be
represented by a sum of its expected value C
ab
and a zero-mean noise term n
C
(ω) which is a
complex random variable,
I R ab C ab
N
ab
in n C n C R (6)
where n
R
(ω) and n
I
(ω) are real and imaginary parts of the random component n
C
(ω) of sample
mean cross-spectrum. The noise term n
C
(ω) approaches zero for infinite number of i.i.d cross-
spectrum realizations, which requires unlimited amount of data.
The expected cross-spectrum can be approximated with two wave packets on causal and anti-
causal parts propagating in opposite directions with the same speed (simplified from Eq. 1 of Liu
et al., 2015),
4
exp
4
exp
i
c
x
i i
c
x
i C
ab
(7)
where α(ω) and β(ω) are two real amplitude terms for causal and anticausal parts, respectively. A
phase shift of π/4 is equivalent to multiplying with i which is due to the approximation of
Bessel functions when the inter-station distance is greater than the wavelength (Liu et al., 2015).
In Eq. (7), the wave speed c(ω) is assumed finite and greater than zero.
Substituting Eq. (7) into Eq. (6), we have the sample mean cross-spectrum
(8)
38
where the real and imaginary parts have β(ω)+α(ω) and β(ω)-α(ω) as their expected envelope
functions, respectively. Applying frequency domain Hilbert transforms to the real and imaginary
parts of the sample mean cross-spectrum respectively, we get their analytic forms
(9)
where and are analytic forms of the real and imaginary parts of the expected
(mean) cross-spectrum, respectively. They are derived based on Bedrosian's theorem (e.g.
Boashash 1992) assuming the time domain representations of β(ω)±α(ω) are zero for
or the envelope of the real or imaginary part cross-spectrum is smoother than
the cross-spectrum itself. The analytic form of either real or imaginary part of cross-spectrum is
equivalent to multiplying its inverse Fourier transform by a step function (in time domain) and is
used for obtaining the envelope. The analytic form of random noises for the real and imaginary
parts are and , respectively
(10)
where Hn
R
(ω) and Hn
I
(ω) are, respectively, the Hilbert transforms of n
R
(ω) and n
I
(ω). The
Hilbert transforms Hn
R
(ω) and Hn
I
(ω) are uncorrelated with n
R
(ω) and n
I
(ω), respectively
(Appendix C).
The wave packets for causal and anticausal parts of the cross-spectrum can be separated
based on Eq. (9),
(11)
39
From Eq. (11) we can derive the mean values and uncertainties of the amplitude estimators
and as well as the phase angle estimator . The amplitude
uncertainties for causal and anticausal parts are, respectively (Appendix B)
(12)
where and are variances for the real and imaginary
parts of sample mean cross-spectrum, respectively. The frequency domain SNR on causal or
anti-causal part cross-spectrum can be conveniently defined as mean over square root of variance:
or .
The phase uncertainty can be derived based on error propagation in the non-linear case
(Appendix B),
(13)
where
2
|
ˆ
R
,
2
|
ˆ and
2
|
ˆ are the variances for real part, causal part and anticausal
part of cross-spectrum, respectively. Each variance term in Eq (13) can be understood as 1/SNR
2
(variance/mean
2
) of wave envelope of the corresponding analytic spectrum. Most studies derive
phase and group velocities from the symmetric component (real part) of the cross-spectrum (e.g.
Lin et al., 2008; Weaver et al., 2009), so the phase uncertainty derived for the real part cross-
spectrum should be used in such cases.
40
2.2.3. Signal-to-Noise Ratio on narrow-band time domain cross-correlation
Assuming phase velocity c(ω) and amplitude terms α(ω) and β(ω) on causal and anticausal
parts are constant in a narrow-band, the time domain analytic cross-correlation function of
narrow-band filtered Eq. (8) can be written as inverse Fourier Transform of the filtered spectrum,
(14)
where the Gaussian filter is ) exp(
~
2 2
a A with a the filter width parameter and 0 a .
Applying convolution theorem and assuming the unfiltered sample mean cross-spectrum satisfy
Hermitian symmetry, the time domain narrow-band analytic cross-correlation becomes,
(15)
where ) 2 /( )) 4 /( exp(
2 2
a a t t A is the inverse Fourier transform of the filter , and *
t
denotes convolution in time domain. Eq. (15) indicates that the causal part peak amplitude of the
analytic narrow-band filtered cross-correlation at frequency ω
0
is equal to ) 2 (
0
a , and
similarly for the anti-causal part the peak amplitude is .
The Root Mean-Square (RMS) estimation of noise in the real time domain cross-correlation
is (Appendix D),
2 2
1
E
2 4
1
RMS
0
2
0
2
0
*
0 0 ,
a
n n
a
I R
C C TD n
(16)
which is related to the filter width a. The narrower the Gaussian filter bandwidth is, the larger a
and consequently the smaller the RMS value become. The noise in the cross-correlation function
is assumed stationary.
Based on Eq. (15) & (16), the conventional SNR measurement in time domain on causal part
can be related to (Appendix A),
41
. (17)
This suggests that the SNR measured from time domain narrow-band filtered cross-correlation
depends on the filter width parameter a, such that if one increases the filter width (decreases a),
the SNR in time domain increases.
Theoretically, the time domain SNR is different from frequency domain SNR (defined below
Eq. 12) because it depends on the bandwidth of the narrowband Gaussian filter. In realistic cases,
however, the noise in the cross-correlation function is not stationary and the RMS noise
estimation is less than in Eq. (16). The RMS noise value also depends on the length/position of
the noise window, the tail of the Gaussian filter in time domain and the finite length of the cross-
correlation.
2.2.4. Stacking of piecewise stationary cross-spectrum
The actual seismic noise data are not stationary and the noise records during different time
periods can have different mean and variance values. Assuming that the noise recording within a
fixed length of time (day/month/etc, for simplicity we use “day” in the following) is stationary,
the mean and variance for either real or imaginary part cross-spectrum X
j
at a single frequency ω
of j.th day are,
. (18)
They can be estimated from the sample mean and variance at the same frequency and day. If
each day can be divided into M non-overlapping and uncorrelated time windows, the average
cross-spectrum of M i.i.d samples and the standard error from j.th day is, respectively
(19)
where t
k
is the time of the k.th time window in a day.
42
Let X
S
be the averaged stack of N days of cross-spectrum,
N
j
j S
X
N
X
1
1
. (20)
Then the mean and variance of X
S
are, respectively
2
1
2
1
2
, 2
2
1
1 1 1 1
Var
1
E
NM N NM N
X
m
N
X m
N
j
j
N
j
stde j S S
N
j
j S S
. (21)
where is the average variance of N different days and NM is the total number of
sample windows in N days. The mean value of the stacked cross-spectrum of N days is the
average of N daily mean values, and the stacked variance is the N-day averaged sample variance
divided by total number of samples NM. There are two ways to reduce the variance of
stacked cross-spectrum X
S
: (1) increase the number of samples NM or (2) remove or rescale the
outliers in samples of (j=1 to N) in order to minimize
2
.
43
2.3. Numerical simulation
#10
-3
-1 0 1
0
200
400
600
800
real part
0.14 Hz
(a)
(b)
(d) (e)
Frequency (Hz)
0.1 0.2 0.3 0.4 0.5 0.6 0.7
Cross-spectra
#10
-5
-4
-3
-2
-1
0
1
2
3
cross-spectra
3600 samples/5.00 days
Real part of averaged cross-spectra
Imaginary part of averaged cross-spectra
Error bar for real part
Error bar for imag part
window 100 s gap 20 s
Synthetic noise
time
(c)
Time (s)
-100 -50 0 50 100
# 10
-6
-2
0
2
(f)
Figure 1. (a) Geometry for numerical simulation of ambient noise on two stations. The inter-
station distance is 50 km and the inter-station attenuation is 30. The background attenuation is
500. The far field noise source distribution is isotropic and source-receiver distance is 8 degrees.
(b) Schematic plot for evenly spaced windows on noise recording. Window length is 100s and
gap is 20s. (c) 3600 independent cross-spectra curves computed from 3600 windows in 5 days.
(d) Stacked cross-spectra with error bars from 3600 sample windows. (e) Histograms of the real
part and joint real & imaginary parts cross-spectrum from 3600 samples. (f) Stacked cross-
correlation of 3600 windows from 5 days of data.
We perform numerical simulations of finite-length (25 days) stationary random noise records
at two stations (Fig. 1a) and compute stacked cross-spectrum of the two stations. The standard
errors of the real and imaginary parts of each cross-spectrum value are computed and the errors
on amplitude and phase are estimated based on the stacked cross-spectrum errors.
The synthetic ambient noise records are computed using Eq. (1) with added uncorrelated
noise. The percentage of uncorrelated noise increases with frequency following 1-exp(-0.35ω),
which is exaggerated compared with the realistic case to illustrate how the uncorrelated noise
44
increase the variance by increasing the power spectrum. The assumed inter-station distance is 50
km and the attenuation Q value between the two stations is 30, which is typical for regional fault
zone environments (Liu et al., 2015). The noise spectrum is computed by summing over 20,000
random noise sources that are 8 degrees away from the center of the two stations and are
assumed i.i.d stationary Gaussian random processes. The inverse Fourier Transform of the noise
cross-spectra produces the corresponding ambient noise in time domain with equal number of
data points. The Nyquist frequency is 1 Hz.
The cross-correlation can be obtained by directly cross-correlating the continuous noise
records of 25 days at the two stations since the synthetic noise data are stationary with no other
signals (e.g. earthquakes) included in realistic continuous waveform recordings. However, we
adopt an alternative method consistent with our theoretical results on stacked cross-spectrum
with error estimation. The continuous noise data on each station are divided by evenly spaced
windows with length of 100 s and gap of 20 s between windows (Fig.1b). The gap of 20 s is
sufficient to ensure that the cross-spectra curves computed from any two windows are
independent for frequencies above 0.1 Hz.
As a first test of the synthetic data and their connections with the theoretical errors, we select
5-day long noise data with 3600 non-overlapping windows, which produce 3600 independent
random cross-spectra curves (Fig. 1c). The stacked cross-spectra of the 3600 independent cross-
spectra curves are shown in Fig. (1d) with error bars for real and imaginary parts. The error bars
for the real and imaginary parts at the same frequency are very similar but the amplitude of the
real part is much higher than on the imaginary part due to the isotopic noise source distribution.
The real part and joint real-imaginary parts histograms of cross-spectrum at 0.14 Hz are shown in
Fig. (1e). The real part cross-spectrum histogram has sample mean of 2.26e-5, which is 10 times
more than its standard error 1.92e-6. The stacked cross-correlation function of the 5-day (3600
windows) data is shown in Fig. (1f).
45
Figure 2.1-σ phase error and errors for its derived quantities. (a) Phase error computed from the
symmetric component (real part) of cross-spectra in Fig. (1d). (b) Phase travel time error. (c)
Phase velocity dispersion curve with 1-σ confidence interval.(d) Relative error in phase velocity.
The phase error (Fig. 2a) can be derived from Eq. (13) using standard errors from the stacked
cross-spectrum errors. The phase error is inversely proportional to the square of the SNR so it's
larger when SNR gets lower. The phase error should be smaller than π/4, otherwise it may cause
phase unwrapping errors for travel time and phase velocity estimates. Assuming the phase error
is much smaller than π, the error in travel time can be conveniently derived from the linear
relation dt=ωdφ (Fig. 2b) and can be used for tomography based on travel time measurements on
surface waves (e.g. Lin et al., 2009). In addition, assuming that the ray path is straight, phase
velocity errors (Fig. 2c & d) can be derived from dc=-c
2
dφ/ωx and phase velocity can be derived
based on Eq. (7) in frequency domain or as in Lin et al. (2008) in time domain. Phase velocity
errors can potentially be used for tomography based on straight ray path approximation and
monitoring temporal velocity changes at different frequencies.
46
Figure 3.Comparison among different preprocessing techniques: no preprocessing (column a),
onebit normalization (column b) and prewhitening (column c). 18000 windows (25 days) are
used. (a), (b) and (c) are joint distributions for real and imaginary parts of cross-spectrum. (d), (e)
and (f) are causal part cross-spectra amplitude curves with error bars.
The effects of several commonly used preprocessing steps on distributions of errors and
amplitude sample mean cross-spectrum are evaluated using the 25-day long noise data with
18000 samples. Figs. (3a), (3b) and (3c) show results for no preprocessing, one-bit normalization
and pre-whitening, respectively. The top panels contain the joint distributions of the real and
imaginary parts of cross-spectrum at 0.14 Hz. The bottom panels show the causal part amplitude
decay curves computed from first line of Eq. (11) with different preprocessing steps. These
results suggest that both one-bit and pre-whitening modify the distribution of cross-spectrum
samples and produce larger errors relative to their corresponding amplitude mean values.
Moreover, the pre-whitening step also normalizes the absolute square of each cross-spectrum
value to one (Fig. 3c)
. (22)
Therefore, the joint distribution of real and imaginary parts with pre-whitening is non-zero only
on the unit circle. The pre-whitening step significantly modifies the amplitude decay curve (Fig.
3f) and does not allow for accurate amplitude estimations.
47
Figure 4. Causal part SNR convergence images for 25 days/18000 windows noise data. (a) and
(b) are SNR convergence images of causal part cross-correlation amplitude estimated in time and
frequency domains, respectively. (c)and (d) are stacked cross-correlations from 18000 samples
for onebit normalization and prewhitening, respectively. (e) Frequency domain SNR
convergence image of causal part cross-correlation amplitude with onebit normalization. (f)
Frequency domain SNR convergence image of causal part cross-correlation amplitude with
prewhitening.
The SNR of time domain narrowband cross-correlation (Fig. 4a) is defined like in theory
section 2.3 as the ratio between the peak of wave packet envelope over the RMS estimation of
48
coda noise . The narrow-band Gaussian filter width is a = 28.29 and the noise window is selected
between 50 s and 100 s. The frequency domain SNR for causal part amplitude cross-spectrum
(Fig. 4b) is defined below Eq. (12) as the ratio of causal amplitude over square root of variance.
Both time and frequency domain SNR plots show that higher SNR are obtained when more
samples are used for each frequency. For a fixed number of samples, the SNR also decreases as
the frequency increases because of the increasing high-frequency uncorrelated noise, relative to
correlated noise, added in the forward noise synthetics. Higher percentage of uncorrelated noise
in the ambient noise recordings results in larger power spectrum values and consequently slower
convergence rate according to Eq. (4). The time domain SNR is unstable and fluctuates with
increasing trend as the number of samples increases while the frequency domain SNR increases
steadily with more samples (Fig. 4a vs Fig. 4b). Moreover, the time domain SNR depends on the
width of the Gaussian filter and the definition of coda noise window, which make it non-unique.
For comparison, the frequency domain SNR is more stable and is a unique statistical estimation
for each discrete frequency value. The time domain cross-correlations of 25 days (18000
windows) are shown in Figs. (4c) & (4d) for one-bit normalization and pre-whitening,
respectively. The frequency domain SNR with one-bit normalization (Fig. 4e) or pre-whitening
(Fig. 4f) shows similar trend of increasing SNR with more samples, however, the maximum
SNR values (14.5 for onebit and 15 for pre-whitening) are much lower than in the case without
preprocessing (Fig. 4b, maximum SNR is 18) because of the nonlinear nature of the one-bit or
pre-whitening preprocessing steps.
Figure 5. (a) Map of three stations CHF, SBB2 and LMR2 in CI network in Southern California
(black triangles). The black lines indicate the fault traces. The green and brown colors show low
49
and high elevation respectively. (b) Cross-correlations of three station pairs computed from 250
days in 2014.
2.4. Application to regional network stations
To illustrate the utility of the theoretical results, we use data recorded by three stations CHF,
SBB2 and LMR2 selected from the regional CI seismic network in Southern California (Fig. 5a).
After removing instrument response, the continuous vertical component recordings for 250 days
in the year 2014 are divided into 4-hour long segments. Segments dominated by significant
earthquakes are removed and any spikes more than 4 times the standard deviation of the noise
segment are clipped to minimize the effects of smaller earthquakes (Poli et al., 2013; Zigone et
al., 2015). Neither one-bit normalization nor pre-whitening is applied in the preprocessing step to
preserve the original amplitude information. The cross-correlation functions of raw 250-day
stack for the three station pairs are shown in Fig. 5b.
50
(a) (b) (c)
(d)
(h)
(e)
(i)
cross-spectrum value
#10
9
-1 0 1 2
# 10
4
0
0.5
1
1.5
2
2.5
real part: 162250 samples.
mean/stderr=-25.2
at 0.24Hz
cross-spectrum value
# 10
8
-1 -0.5 0 0.5 1
# 10
4
0
0.5
1
1.5
2
2.5
real part: 154137 samples.
mean/stderr=-30.1
at 0.24Hz
(f) (g)
frequency (Hz)
0 0.1 0.2 0.3 0.4 0.5 0.6
stacked cross-spectra
# 10
7
-4
-2
0
2
4
6
8
Selected
Real part of xspectra
Imaginary part of xspectra
Error for real part
Error for imag part
frequency (Hz)
0 0.1 0.2 0.3 0.4 0.5 0.6
stacked cross-spectra
# 10
8
-2
-1
0
1
2
3
Outliers
Real part of xspectra
Imaginary part of xspectra
Error for real part
Error for imag part
Selected Outliers
frequency (Hz)
0 0.1 0.2 0.3 0.4 0.5 0.6
cross-spectra/stderr
-60
-40
-20
0
20
40
60
80
Outliers
Real part of xspectra/Std
Imaginary part of xspectra/Std
frequency (Hz)
0 0.1 0.2 0.3 0.4 0.5 0.6
cross-spectra/stderr
-100
-50
0
50
100
Selected
Real part of xspectra/Std
Imaginary part of xspectra/Std
Figure 6. (a) Raw distribution of 162250 real part cross-spectrum values from same number of
windows at 0.24 Hz for station pair CHF-SBB2. (b) Distribution of real part cross-spectrum
(154137 samples) for CHF-SBB2 after outlier removal. (c) Joint distribution of real and
imaginary parts of cross-spectrum for CHF-SBB2 after outlier removal. (d) First 2000 cross-
spectra curves from 115661 selected windows. (e) First 2000 cross-spectra curves from 46589
outlier (unselected) windows. (f) Stacked cross-spectra of CHF-SBB2 for selected windows
51
(115661 windows). Standard deviation curve for real part cross-spectra is plotted in black dashed
line. (g) Stacked cross-spectra of CHF-SBB2 for outlier windows (46589 windows). Note the
difference below 0.1 Hz. (h) Normalized cross-spectra of CHF-SBB2 for selected windows. The
variance of cross-spectrum at each frequency is normalized to one. (i) Normalized cross-spectra
of CHF-SBB2 for outlier windows with unite variance.
The error analysis primarily focuses on the station pair CHF-SBB2. As in the synthetic test,
the non-overlapping window length is 100 s with gap of 20 sec between windows (see sketch in
Fig. 1b). The 250-day data is equivalent to 162250 windows (samples) after earthquake segment
removal. The histogram of the real part cross-spectrum at 0.24 Hz is shown in Fig. (6a) and
suggests the existence of many outliers because of the long tails of the histogram. The samples
outside of three Median Absolute Deviation (MAD) are assumed to be outliers and removed.
After outlier removal, the new distribution of 154137 remaining samples of real part cross-
spectrum at the same frequency has larger ratio of mean over standard error (Fig. 6b). The
outliers may involve cross-spectrum of earthquakes coda or strong oceanic storms, which have
different spectral density compared with the regular noise energy. Therefore, it is better to
exclude outliers before estimating standard errors. The joint distribution of the real and
imaginary parts of the cross-spectrum at 0.24 Hz after outlier removal (Fig. 6c) indicates that the
real and the imaginary parts do not show visible correlation and have similar variance. These
results are consistent with the pseudo-variance results when comparing Eq. (5) with Eq. (B1).
Outliers for each frequency can be identified based on MAD in a similar way as in Figs. (6a)
& (6b). For each window, the percentage of outliers of all frequencies in the frequency passband
is an indicator of the amount of abnormal signal within the window. We define a threshold of
maximum 5% outliers of all frequencies in each time window selected for stacking. The 115661
selected windows produce the same number of random cross-spectra curves (Fig. 6d) while the
46589 unselected (outlier) windows give random cross-spectra curves (Fig. 6e). A simple stack
of the windows produces cross-spectra with associated standard errors computed from the
distribution at each frequency (Fig. 6f). Similarly, Fig. (6g) shows the stacked cross-spectra for
outlier windows, which presents a very different mean values and standard errors at all
frequencies within the passband since the outlier windows follow different statistics. The stacked
cross-spectra curves derived from both selected windows and outlier windows have
characteristic peak at 0.15 Hz due to the secondary ocean microseismic peak, and the amplitude
52
of the spectrum decays quickly above that peak frequency. The cross-spectra curve derived from
outlier windows has a smaller peak near 0.07 Hz, probably due to the primary microseismic peak.
Figure 7. (a) and (b) are SNR convergence images of causal part cross-correlation amplitude for
250 days estimated in time and frequency domains, respectively. Only selected windows
(defined in Fig. 6f) are used. (c) and (d) are SNR convergence curves sliced at different
53
frequencies (0.15, 0.20, 0.32 Hz) for time and frequency domain measurements, respectively. (e)
and (f) are SNR versus frequency curves sliced at different days (5, 125, 250 days) for time and
frequency domain measurements, respectively.
Dividing the stacked cross-spectrum by its standard error, the normalized cross-spectrum has
a unit variance. This procedure is connected to the approach of Aki (1957) because the standard
error is theoretically proportional to the square root of the product of power spectrums on two
stations (Eq. 4), which is the denominator in the normalization formula of Aki (1957). This
normalization produces cross-spectra curve with white additive noise and more evenly
distributed energy within the frequency band, which is helpful for frequency domain Hilbert
transform (better satisfy the condition for Bedrosian's theorem) or isolating fundamental mode
surface wave packet in time domain. Because of the strong variation of the ambient noise power
spectra due to ocean microseismic peaks, we normalize the stacked cross-spectrum by its
standard error before Hilbert transform and multiply back the standard error after Hilbert
transform. The normalized cross-spectra curve for selected windows (Fig. 6h) shows higher
SNRs than that for outlier windows (Fig. 6i).
The time and frequency domain SNR images of the causal part cross-correlation amplitude
(Fig. 7a and Fig. 7b, respectively) are computed in the same way as in the synthetic test. The
color scale represents the SNR values. Only the selected windows within the number of days are
used for stacking. The time domain SNR values are measured with the same Gaussian filter and
noise window parameters as in the synthetic test. The time domain SNR analysis shows
increasing SNR values at multiple frequencies with increasing number of days, but the values are
unstable and get saturated (stop increasing after certain number of days) at some frequencies. In
contrast, the frequency domain SNR analysis shows consistently increasing SNR with number of
days at multiple frequencies. Both images show generally consistent peaks and troughs of SNR
values at different frequencies, suggesting multiple noise sources at different frequencies other
than the main ocean microseism peak.
By slicing both SNR images at three different frequencies (0.15 Hz, 0.20 Hz and 0.32 Hz),
the corresponding SNR values are plotted as functions of number of days in Fig. (7c) and Fig.
(7d) for time and frequency domain SNR, respectively. At 0.15 Hz and 0.32 Hz, the SNR values
increase with number of days because they correspond to the noise peaks in panels 7a & 7b. At
0.20 Hz, however, the SNR stops increasing beyond 35 days probably because of a gap between
54
two noise sources in which the noise random process is highly non-stationary. By slicing the
SNR images in Fig. (7a) & (7b) at three different day counts, we plot in Fig. (7e) and Fig. (7f)
SNR functions of center frequency for time and frequency domain SNR measurements,
respectively. For 5 days stacking, the time domain SNR show anomalous peak near 0.14 Hz (also
in Fig. 7a), probably because the wave packet peak is not well defined due to the actual low SNR
value (the maximum after 5 days is only 40). For 125 and 250 days of stacking, the time domain
SNR curves show similar patterns and SNR values but at some frequencies (e.g. 0.14 Hz, 0.24
Hz) the SNR values get saturated while at some other frequencies (e.g. 0.32 Hz, 0.52 Hz) the
SNR values keep increasing up to 250 days. For frequency domain SNR (Fig. 7f), all three SNR
curves for 5, 125 and 250 days show similar patterns of peaks and troughs at different
frequencies and progressive increase in SNR values as the number of stacked days increases.
These results show more stable SNR estimates from frequency domain method probably because
the frequency domain method uses variance of cross-spectrum in addition to mean (the stacked)
cross-spectrum.
Figure 8. Block bootstrap (hourly block) estimation of empirical errors of 250-day cross-
spectrum for CHF-SBB2. (a) Temporal correlation of noise revealed by autocorrelation of real
55
part cross-spectrum from different windows. At 0.16 Hz there is weak correlation ~0.07; at 0.20
Hz there is nearly zero correlation (<0.01). (b) Distribution of bootstrap cross-spectrum at 0.16
Hz and 0.24 Hz, respectively. (c) Ratio of real/imaginary part block bootstrap error (more
realistic) over standard error (assuming i.i.d samples from different windows). Bootstrap method
suggests larger error between 0.05-0.2 Hz due to temporal correlation. The ratio approaches to
one for frequency between 0.2-0.6 Hz.(d) SNR convergence image of causal part cross-
correlation computed using hourly block bootstrap resampling. Because of temporal correlation,
the bootstrap SNR between 0.05-0.2 Hz are lower than the frequency domain SNR based on
standard error.
We estimate the correlation between different noise windows by computing the
autocorrelations of the real part cross-spectrum values for two different frequencies (Fig. 8a).
The temporal correlation among windows can add extra cross-window (non-zero covariance)
terms in calculating standard error for real and imaginary parts of stacked cross-spectrum. As a
result, temporal correlation of cross-spectrum values from different time windows will increase
the variance of stacked cross-spectrum. As shown on Fig 8a, the correlation coefficient is around
0.07 at 0.16 Hz from offset of 1 window to 49 windows, which suggest weak temporal
correlation longer than one and half hours. At 0.20 Hz, however, the real part of the cross-
spectrum is not correlated (correlation coefficient ±0.01 or less). For higher frequencies (between
0.2-0.6 Hz), the cross-spectrum value does not show temporal correlations among different
windows.
An hourly block bootstrap method is applied as an independent way to estimate the empirical
confidence interval for the cross-spectrum derived from the same selected windows (115661) for
pair CHF-SBB2 as in the previous SNR example. The purpose of hourly block bootstrap is to
capture some temporal correlation among the neighboring windows if they exist. For each one-
hour block, cross-spectra curves from all windows within the block are averaged. Bootstrap
resampling is applied to all hourly blocks with replacement. At each frequency, a bootstrap
cross-spectrum sample is computed by averaging randomly drawn block cross-spectrum values.
There are 4000 bootstrap cross-spectrum samples at each frequency (Fig. 8b).
The ratio of 250-day block bootstrap standard deviation over standard error (assuming i.i.d
windows) for the real and imaginary parts of cross-spectra are plotted in Fig. 8d. The ratio for
both real and imaginary parts are approaching to one above 0.2 Hz, suggesting that the standard
error based on the assumption of i.i.d cross-spectrum samples at each frequency from different
windows offers reliable estimation of the data variance between 0.2-0.6 Hz. The ratio for the real
56
and imaginary parts between 0.05-0.2 Hz show significant oscillation between 2.1 and 1 with
downward trending, consistent with the fact that the cross-spectrum within this frequency range
corresponds to weak temporal correlations among windows (Fig. 8a).
Applying block bootstrap to different number of days (1-250), the mean and standard
deviation of bootstrap distribution for each cross-spectrum value can be easily computed and the
SNR for the causal part cross-spectrum can be derived (Fig. 8d) in a similar way as frequency
domain SNR based on Eq. (12) by replacing standard error with bootstrap error. For frequency
above 0.2 Hz, the bootstrap SNRs (Fig. 8d) are similar (very close) to the frequency domain
SNR values (Fig. 7b) based on standard errors assuming i.i.d cross-spectrum samples at each
frequency from different windows. For frequency between 0.05-0.2 Hz, however, the SNR
values derived from bootstrap are evidently lower (~20% less) than the values in Fig. 7b based
on standard error, indicating temporal correlation among the windows that increase the variance
within this frequency range. Because the temporal correlation at some frequency (e.g. 0.16 Hz) is
greater than one hour, applying block bootstrap with longer block length can further increase the
variance and reduce the corresponding SNR value.
57
Figure 9.Phase and amplitude errors (1-σ) of cross-spectra based on block bootstrap resampling
of 115661 selected windows. The error estimations below 0.2 Hz are less accurate due to
temporal correlation. (a) Phase error curve of real part cross-spectra. (b) Phase travel time error.
(c) Phase velocity dispersion curve with 1-σ confidence intervals. (d) Relative error in phase
velocity. (e) Causal part amplitude cross-spectra with error bars. (f) Zoom-in version of causal
part amplitude cross-spectra in panel e.
The phase error curve (Fig. 9a) is derived directly from the bootstrap error of the real part
cross-spectrum according to Eq. (13). The maximum phase error (1-σ) is 0.11 radians, which is
much less than phase angle period 2π. This indicates that the phase angle can be reliably
unwrapped without 2π ambiguities. Following similar procedure as in synthetic test (section 3),
the phase travel time error can be estimated directly from the phase error. Assuming straight ray
path, the phase angle errors can be mapped into phase velocity uncertainties (Fig. 9 c&d) and the
largest uncertainty is 0.4 % velocity variation within one-sigma interval. The causal part
58
amplitude cross-spectrum (Fig. 9e & f) is estimated based on Eq. (11) and the corresponding
confidence interval is computed based on the variance values of the real and imaginary parts of
cross-spectrum in Eq. (12). The amplitude errors below 0.2 Hz are underestimated because of
temporal correlation of noise. By comparison with the bootstrap SNR figure (Fig. 8d), high SNR
indicate reliable estimation of causal part cross-spectrum amplitude.
2.5. Application to densely spaced stations with high frequency data
time (s)
-1 -0.5 0 0.5 1
-0.5
0
0.5
JF00-JFS1: dist=21.6 m
time (s)
-1 -0.5 0 0.5 1
-0.1
0
0.1
JF00-JFS2: dist=35.6 m
time (s)
-1 -0.5 0 0.5 1
-0.2
0
0.2
JFS1-JFS2: dist=14.5 m
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
JF00 Los$
Angeles$
frequency (Hz)
15 20 25 30 35 40 45 50
stacked cross-spectra
# 10
5
-6
-4
-2
0
2
4
6
JF00-JFS1: raw stack of 1 day
Real part of xspectra
Imaginary part of xspectra
Error for real part
Error for imag part
frequency (Hz)
15 20 25 30 35 40 45 50
stacked cross-spectra
# 10
5
-2
-1
0
1
2
JF00-JFS1: no outliers, 1 day
Real part of xspectra
Imaginary part of xspectra
Error for real part
Error for imag part
(a)
(d)
(c)
(g) (h)
6220 selected windows 7332 windows
window 10 s 1 s gap
High frequency
noise
time
(e)
1221 outliers removed
(f)
cross-spectrum value
# 10
8
-1 0 1 2
0
500
1000
1500
2000
2500
Real part: 7332 samples
mean/stderr=3.22
at 21.0 Hz
cross-spectrum value
# 10
6
-5 0 5
0
500
1000
1500
2000
Real part: 7185 samples
mean/stderr=13.14
at 21.0 Hz
(b)
Figure 10.Error analysis experiment using one-day high frequency noise (10-50 Hz) on JF array
with small station spacing (~20 m).(a) Location map of dense linear JF array across the San
Jacinto fault zone. (b) Cross-correlations of three stations JF00-JFS1-JFS2 computed from the
59
2nd days in 2012. (c) Schematic plot of evenly spaced windows with 10 s window length and 1 s
gap. (d) Raw distribution of 7332 real part cross-spectrum values from same number of windows
in one day for station pair JF00-JFS1. (e) Distribution of real part cross-spectrum (7185 samples)
for JF00-JFS1 after outlier removal. (f) First 2000 of 6220 selected cross-spectra curves from
6220 windows. (g) Raw stack of cross-spectra curves computed from all 7332 windows in the
2nd day of 2012. (h) Stacked cross-spectra of 6220 selected cross-spectra curves (panel e) from
same number of windows after removing outliers.
As an additional illustration, we compute noise cross-correlations using data recorded by
three stations JF00-JFS1-JFS2 of a dense linear array (JF) across the San Jacinto fault (Fig. 10a).
The cross-correlation functions of first 10 days in the year 2012 are plotted in Fig. (10b). Evenly
spaced windows with length of 10 s and gap of 1 s are used to compute cross-spectra curves
within the frequency range of 10-50 Hz. Because of the high variability of the noise intensity at
high frequency (Ben-Zion et al., 2015), we select only one day (the 2
nd
day) of 2012 and
compute cross-spectra for station pair JF00-JFS1. The histogram of the real part cross-spectrum
at 21 Hz with 7332 samples from same number of windows (Fig. 10c) has a ratio of mean over
standard error of 3.22. After removing the outliers according to 3 MAD value and zooming in
around the selected samples, the new distribution (Fig. 10d) of real part with 7185 samples has a
higher ratio of mean over standard error of 13.14. A selected window should contain no more
than 5% of outliers from all frequencies within 10-50 Hz. There are 6220 selected windows that
produce different cross-spectra curves (Fig. 10e). Therefore, the raw stacked cross-spectra before
outlier removal (Fig. 10f) have more fluctuations and larger error bars relative to spectral
amplitude than the stacked cross-spectra after outlier removal (Fig. 10e).
Figs. (11a & b) compare the amplitude SNR values based on raw stacked (with outliers)
cross-spectra and outlier-removed cross-spectraobtained with one-day causal part (and anticausal
part) cross-spectra for pair JF00-JFS1. After outlier removal, the SNR of causal part amplitude is
4 to 10 times greater than that before outlier removal between 17-45 Hz (Fig. 11a). For
anticausal part amplitude, the outlier-removed SNR is 4-10 times better than outlier-included
SNR between 10-23 Hz (Fig. 11b). In addition, after outlier removal the anticausal amplitude
SNR decays faster than the causal amplitude SNR, suggesting that different
source/scattering/propagation mechanisms for noise coming from two opposite directions may
affect the amplitude SNR differently.
60
Figure 11. Amplitude SNR and phase error based on one-day stack of cross-spectra of JF00-
JFS1. For each panel, the information derived from raw stack of cross-spectra with outliers are in
blue and those from stacked cross-spectra without outliers are in red. (a) Causal part amplitude
SNR. (b) Anticausal part amplitude SNR. (c) Phase error. (d) Phase travel time error derived
based phase error. (e). Phase velocity curves. (f) Phase velocity relative error.
The outliers also affect phase velocity retrieved from one-day stack of cross-spectra for the
pair JF00-JFS1. The phase errors (Fig. 11c) estimated from real part of stacked cross-spectra (Fig.
10g) with outliers are ~4 times greater than those without outliers and show spikes greater than
pi/4 at 27, 32-35 and 38 Hz, which may cause phase unwrapping errors. Based on phase error
information, the phase travel time error (Fig. 11d) and phase velocity error (Fig. 11f) can be
derived following the same way as in synthetic test (section 3). The phase velocity dispersion
curves based on stacked cross-spectra with and without outliers show 1-4 % discrepancy
61
between 11-20 Hz (Fig. 11e), which is probably due to systematic errors brought by outliers
across different frequencies.
Figure 12.Phase velocity and relative phase velocity error based on one-month stack of cross-
spectra for three station pairs in JF00-JFS1-JFS2. For each panel, the information derived from
raw stack of cross-spectra with outliers are in blue and those from stacked cross-spectra without
outliers are in red. (a) & (b) are phase velocity and corresponding relative error on JF00-JFS1. (c)
& (d) are phase velocity and corresponding relative error on JF00-JFS2. (e) & (f) are phase
velocity and corresponding relative error on JFS1-JFS2.
Phase velocity curves with confidence intervals (Fig. 12) are estimated from the stacked
cross-spectra of first month in year 2012. The one-month stack of the pair JF00-JFS1 yields
much smaller phase velocity errors (Figs. 12 a & b) compared to the phase velocity errors from
the 2nd day of 2012. The spikes in the phase velocity error (Fig. 12b) with outliers are also
reduced on the one-month stacked data. After outlier removal, however, the phase velocity
62
errors associated with different frequencies are reduced by a factor of 2-6. The phase velocity
error curves for two additional pairs, JF00-JFS2 (Figs. 12c & d) and JFS1-JFS2 (Figs. 12e & f)
show improved phase error estimations and less spikes after outlier removal.
2.6. Discussion
Error analysis of ambient noise cross-correlation can provide important information on the
deviation of the finite-length noise cross-correlation from the expected cross-correlation that is
linked, through time derivative, to the empirical Green’s function. We develop a formulation to
estimate the error of stacked cross-spectrum at each frequency by calculating standard or
bootstrap errors from cross-spectrum values of evenly spaced non-overlapping windows at the
corresponding frequency. The variance values for phase and amplitude cross-spectrum are
derived from variances of real and imaginary parts of stacked cross-spectrum. From phase angle
uncertainty, phase travel time and phase velocity errors can be easily derived, which can help to
better constrain uncertainties in phase velocity temporal changes and in phase velocity
tomography studies. In addition, these theoretical results can be used for uncertainty estimations
of amplitude spectra for either causal or anti-causal part of cross-correlation. Such amplitude
uncertainties may help inferring attenuation coefficients or site amplifications with confidence
intervals, which is a key ingredient in future noise-based attenuation tomography (Liu et al.
2015).
This work complements previous studies on time domain noise cross-correlation. These
studies (e.g. Snieder 2004; Sabra et al. 2005a) state that the variance of noise is proportional to
the product of two autocorrelations at two stations and inversely proportional to the length of
noise recording. However, real ambient noise data are not stationary and contain earthquakes and
other transient signals, leading to complex pre-processing and segment cuttings steps (e.g.
Bensen et al. 2007; Seats et al. 2012; Zigone et al. 2015), which make the time domain variance
formulae less practical. The time domain noise variance is commonly estimated from coda of
noise cross-correlation, while the frequency domain noise variance can be derived directly from
the standard error of numerous i.i.d cross-spectrum samples at the same frequency in different
windows. Therefore, the error of stacked cross-spectrum at each frequency is easier to derive and
can provide more information than the time domain variance, as its a statistical quantity that can
be derived independently from power spectrum functions in the Eq. (4).
63
At the same frequency, the cross-spectrum values from non-overlapping time windows are
assumed i.i.d, which is a valid assumption for stationary noise data provided that the gap
between two windows is greater than the noise autocorrelation width. The independence of those
values can be verified by computing autocorrelation of either real or imaginary part of the cross-
spectrum sequence from different windows. As mentioned, real ambient noise is not a stationary
random process and non-stationary noise can increase the stacked cross-spectrum variance,
thereby reducing the SNR.
Synthetic tests with 25-day stationary noise produced from Gaussian noise model are used to
study statistical properties of sample mean cross-spectrum under ideal conditions. As can be
inferred from the histograms of real and imaginary parts cross-spectrum samples at the same
frequency from non-overlapping windows (Fig. 1e), the distribution of cross-spectrum is not
Gaussian although it is a product of two Gaussian random variables. By increasing the number of
samples, the distribution remains the same (due to the stationary property) but the standard error
decreases.
Effects of two preprocessing techniques, onebit normalization and pre-whitening, are
compared to the case without preprocessing normalization in terms of distribution of random
cross-spectrum values and causal part amplitude decay. The results show that both pre-
processing methods modify the distribution of cross-spectrum samples and increase the errors by
~15% (decrease the SNRs) relative to their corresponding amplitude mean values. In addition to
larger amplitude errors, the pre-whitening also modifies the shape of the amplitude decay curve
(Fig. 3f), preventing accurate amplitude estimations. Together, these results highlight the strong
influence of the preprocessing procedures on the estimation of any amplitude related quantities.
The proposed methodology allows assessing these errors, which is useful for computing
confidence intervals for any amplitude measurements. Such information can provide guidelines
for attenuation tomography and site amplification factor estimations.
In ambient noise cross-correlation study, the SNR convergence provides important
information on the data quality of cross-spectrum at different frequency. The stacked cross-
spectrum function also has information about SNR at its peaks but not at other frequencies. The
time domain SNR estimates are less stable than their frequency domain equivalents in our
64
synthetic tests, and the time domain SNR is non-unique depending on the narrowband filter
width and coda noise window definition.
The distribution of cross-spectrum from real data contains outliers that should be removed to
enhance SNR. The outliers may correspond to clipped amplitudes, small earthquakes, earthquake
coda waves or other transient signals with different statistics. By removing the outlier windows,
the SNR can potentially be improved according to Eq. (21) because the outliers have very
different means and standard errors than the collection of selected windows. Converged cross-
spectra at two stations with lower SNR can be obtained from outlier windows of regional array
CHF-SBB2, as most of the outliers are likely from scattered earthquake coda waves or
tremor/ocean storm which may correspond to higher noise amplitude. This method can
potentially help to identify those tremor-like activities.
The SNR convergence image for causal part amplitude of station pair CHF-SBB2 shows
peaks and troughs of SNR values at different frequencies. The troughs are probably due to gaps
among different noise sources and should be excluded from cross-correlation amplitude
estimations. The time domain SNR convergence image is again less stable than the frequency
domain SNR image. However, for frequency below 0.2 Hz, the time domain SNR image
suggests the SNR should be relatively lower than predicted by the frequency domain SNR image
because of the temporal correlation at these frequencies.
One potential problem of the error analysis assuming i.i.d cross-spectrum values from
different windows is temporal correlation of the noise between 0.05-0.2 Hz, probably related to
the ocean microseismic peaks (primary at ~0.06 Hz and secondary at ~ 0.15 Hz). We apply
hourly block bootstrap method (end of section 4) to account for temporal correlation of the noise
and achieve more realistic estimation for the cross-spectrum error that increases due to temporal
correlation between 0.05-0.2 Hz. We note that the temporal correlation time scale can be several
hours for cross-spectrum values around the ocean microseismic peaks, and longer block for
bootstrap should be used for more accurate error estimation at those frequencies. For cross-
spectrum between 0.2-0.6 Hz, the temporal correlation is close to zero (<0.01) and the bootstrap
error converges to the standard error assuming i.i.d cross-spectrum samples from different
windows. The frequency band of temporal correlation coincides with correlated neighboring
frequency band (0.07-0.2 Hz) near the secondary microseismic peak (Liu and Ben-Zion 2016).
65
Based on the bootstrap cross-spectrum error, the phase velocity error for 250-day stack on
selected windows is less than 0.23% for frequency between 0.15-0.6 Hz. To understand temporal
variations of phase velocity at different frequencies by stacking fewer number of days, the
variation of phase velocity should be statistically significant compared to phase velocity error to
claim temporal change of phase velocity at certain frequency. The phase error can be reduced in
time domain by applying a narrow-band Gaussian filter because wider filter bandwidth can
increase the time domain SNR (e.g. Snieder 2004; Sabra et al., 2005a; Eq. 17 in this paper).
Error analysis is also applied to high frequency (10-50 Hz) noise cross-correlations in San
Jacinto fault zone with inter-station distance of ~ 20 m. A shorter window length (10 s) is used
due to the high frequency content and small distance. The results of one-day stack suggest that
after outlier removal, the amplitude SNRs of cross-spectra are improved by a factor of 4-10 and
phase velocity errors are reduced to ~1/4 with less fluctuations. Because of the high variability of
high frequency noise intensity over time, the error estimation for multiple days stacking is more
difficult than for regional arrays and we only show phase velocity results with confidence
intervals for one-month stack of noise cross-correlation of data from JF00-JFS1-JFS2. The phase
velocity errors before outlier removal are ~2-6 times greater than those after removing outliers.
The presented error analysis method is promising for analyzing the error of stacked noise
cross-correlation data at various scales and frequencies with applications to amplitude
(attenuation) and phase (velocity). More work is needed to account for the strong variability of
the high frequency noise intensity with time. One possibility of addressing this is using a
weighted average of daily/hourly cross-spectrum assuming piecewise stationary noise as
suggested in section 2.4. This and additional related work is left for a future study.
Appendix 2A: Derivation of variance and pseudo-variance of sample mean cross-spectrum
The variance of the sample mean cross-spectrum is by definition
(A1)
Combined with Eq.(2), the correlation of sample mean is,
66
N
p n
p
b
n
b
p
a
n
a
N
n
p
b
p
a
N
n
n
b
n
a
N
p
p
b
p
a
N
n
n
b
n
a
N
ab
N
ab
d d d d
N
d d d d
N
d d
N
d d
N
R R
1
* *
2
1
*
1
*
2
1
*
1
* *
E E
1
E E
1
1 1
E E
(A2)
Substituting Eq. (A2) back into Eq. (A1), we get the variance of the sample mean cross-spectrum
bb aa ab ab bb aa ab ab
N
ab
A A
N
C C A A
N
C C
N
N
R
1 1
Var
* *
2
2
(A3)
where the expected power spectrum on station a is and it takes a similar
form on station b.
The pseudo-variance of the sample mean cross-spectrum is defined as
(A4)
and similarly, when combined with Eq. (2), the pseudo-correlation is
(A5)
Substituting Eq. (A5) back into Eq. (A4), we get the pseudo-variance of the sample mean cross-
spectrum
ab ab ab ab ab ab ab ab
N
ab
C C
N
C C C C
N
C C
N
N
R
1 1
pVar
2
2
(A6)
67
Appendix 2B: Derivations for errors of amplitude and phase from the 2nd moment
statistics of the sample mean cross-spectrum
Based on Eq. (6), variance and pseudo-variance of sample mean cross-spectrum can be
related to the standard errors of its real and imaginary parts,
(B1)
where and are variances for the real and imaginary parts of sample mean cross-
spectrum, respectively. Combining Eq. (B1) with Eq. (5), we have
ab ab I R
C C
N
n n
1
Im E 2 (B2)
which suggests that n
R
(ω) and n
I
(ω) are weakly correlated.
Next the variance of causal part estimator is
(B3)
where the cross terms have no correlations due to the Hilbert transform (see Appendix C):
0 2 E
I R
Hn n and . The Hilbert transform also preserves the
distribution of the zero-mean random variables n
R
(ω) and n
I
(ω) if they are white noises, which
can be achieved by normalizing the stacked cross-spectrum at each frequency by its standard
error (e.g. Aki 1957). If the noise terms on the stacked cross-spectrum function are not white, the
distributions of Hilbert transform of n
R
(ω) and n
I
(ω) are only approximately the same as those of
n
R
(ω) and n
I
(ω). The variance for anticausal part estimator can be derived in a similar way,
68
(B4)
The phase angle information can be derived from the analytic form of the real part of sample
mean cross-spectrum in Eq. (9),
(B5)
The random variables on the right-hand side are n
R
(ω) and its Hilbert transform Hn
R
(ω), which
are not correlated (according to Appendix C).
The phase variance can be estimated with error propagation equation,
(B6)
Evaluating partial derivatives in Eq. (B6) and knowing n
R
(ω) has the same variance with its
Hilbert transform Hn
R
(ω), the phase variance becomes,
(B7)
where is simply the envelope of the real part of sample mean cross-
spectrum. Eq. (B7) can be viewed intuitively as inverse square of SNR of the real part sample
mean cross-spectrum.
Following similar approaches, the phase variances for causal and anticausal parts are,
respectively,
69
(B8)
Appendix 2C: Correlation between noise and its Hilbert transform in frequency domain
Suppose there are two zero-mean real random functions X(f) and Y(f) in frequency domain.
The covariance function of them is a delta function,
, (C1)
which is equivalent to a diagonal covariance matrix in discrete frequency case. The amplitude of
the delta function is scaled by a real function D
XY
(f
1
).
The covariance function is simply defined as
, (C2)
where the new variable w is defined as w=t-t’. Here x(t) and y(t) are inverse Fourier Transforms
of X(f) and Y(f), respectively. Combining Eq. (C1) and (C2), we have
, (C3)
which is a real function of frequency f
1
.
The Hilbert transform of Y(f) is,
HY f
( )
= H Y f
( )
é
ë
ù
û
= H F y t
( )
é
ë
ù
û
é
ë
ù
û
= F -isgn t
( )
y t
( )
é
ë
ù
û
, (C4)
70
Therefore the covariance function of X(f) and HY(f) is,
. (C5)
Therefore X(f) and HY(f) are not correlated, which still holds if the noise covariance function
is a real symmetric function of .
Appendix 2D: Estimating noise and SNR on narrow-band filtered cross-correlation in time
domain
Assuming the real part of the narrow-band filtered noise on the time domain cross-correlation
is widely-sense stationary, its correlation function is,
(D1)
where S
n,flt
(ω;ω
0
) is the power spectral density of the noise narrow-band filtered at center
frequency ω
0
.
The real part of the narrow-band filtered noise term on time domain cross-correlation
according to Eq. (15) is,
. (D2)
Applying inverse Fourier Transform to Eq. (D2), the narrow-band filtered noise spectra is
. (D3)
Its power spectral density is,
(D4)
Substituting Eq. (D4) back into Eq. (D1), the correlation function can be evaluated as,
71
(D5)
Assuming the noise on real time domain cross-correlation in Eq. (D2) is an ergodic random
process, the Mean-Square of it is,
(D6)
Based on Eq. (D6), the Root Mean-Square (RMS) of the time domain noise is,
(D7)
By taking the ratio of causal part envelope amplitude (from Eq. 15) and RMS of real noise
(Eq. D7) on time domain cross-correlation, we arrive at the SNR of the causal part cross-
correlation in time domain,
0
2
0
2
0
0
2 1
SNR
I R
a
. (D8)
Similarly the SNR of the anticausal part cross-correlation in time domain is,
(D9)
Beware that in realistic cases the cross-correlation between two receivers only have finite
length (usually a little longer than the latest arrival; so the coda noise is not stationary and it
decays towards both ends) and therefore the RMS noise estimation based on coda of cross-
correlation can be much less than the theoretical RMS value in Eq. (D7).
72
3. Theoretical and Numerical Results on Effects of Attenuation on Correlation Functions of
Ambient Seismic Noise
(Liu and Ben-Zion, 2013)
3.1. Introduction
Seismic imaging techniques using the ambient noise field have developed considerably in
the last decade (e.g., Gouédard et al. 2008; Campillo et al. 2011, and references therein) and
proved successful in various seismological applications. The origin of these techniques can be
traced to the pioneering studies by Aki (1957) and Claerbout (1968). Based on spectral analysis
of stochastic waves from stationary random processes, Aki (1957) noted that the normalized
cross-spectra at two locations resemble a 0-order Bessel function that depends on the seismic
velocity between the stations. Claerbout (1968) suggested that the Green’s function can be
obtained from the temporal average of the correlation function at different positions. Numerous
recent studies demonstrated that the Green’s function between stations can indeed be retrieved
from the time derivative of the noise cross-correlation function (e.g. Lobkis and Weaver 2001;
Shapiro and Campillo 2004; Campillo 2006). A complete retrieval of the Green function requires
that the noise wavefield is fully diffuse and the noise sources are stationary (e.g., Lobkis and
Weaver 2001; Weaver and Lobkis 2004; Sánchez-Sesma et al. 2008). However, even when these
assumptions are not satisfied, noise cross-correlation can still be used to obtain apparent velocity
information associated with various scales and seismic phases (e.g. Shapiro et al. 2005; Roux et
al. 2005; Bensen et al. 2007).
The ambient noise wavefield in the Earth is generated by sources having strong spatio-
temporal variations (e.g. Stehly et al. 2006; Kimman and Trampert 2010; Landès et al. 2010;
Hillers and Ben-Zion 2011). The lack of isotropic far field noise distribution produces
asymmetric cross correlation functions, with different amplitudes of the causal and anti-causal
parts corresponding to opposite propagation directions. As demonstrated in this paper, seismic
attenuation can amplify and modify the asymmetry of cross correlation functions in the presence
of non uniform source distribution and/or non uniform material properties. Intrinsic attenuation
is related to the fluctuation-dissipation theorem which is fundamental to the noise cross-
correlation of random wavefield (e.g. van Tiggelen 2003). However, retrieving attenuation from
73
noise cross-correlation functions is considerably more difficult than obtaining velocity
information.
Several recent works attempted to analyze the amplitude of noise cross-correlation
functions in relation to attenuation. Prieto et al. (2009) proposed a generalization of the inference
of Aki (1957) where the normalized cross-spectra, referred to as coherency, follows the 0-order
Bessel function scaled by x exp with being an attenuation coefficient and x denoting the
inter-station distance. Several studies (e.g. Lawrence and Prieto 2011; Prieto et al. 2011; Lin et al.
2012) attempted to extract attenuation from empirical coherencies based on this idea. Tsai (2011)
examined the conjecture of Prieto et al. (2009) assuming homogenous attenuation, and showed
that it is valid for uniform sources everywhere including the near field, but is not an appropriate
solution for noise generated by far-field isotropic sources. Weaver (2011b) demonstrated that
strong directional noise distribution can mask the effects of attenuation on the coherency.
In this paper we examine with analytical and numerical results several effects of seismic
attenuation in a fully diffuse wavefield on properties of cross-correlation functions. We consider
a situation where the region between two stations can have a different seismic attenuation factor
than the outside region. The analytical derivations (section 2) assume the Fresnel approximation
and limited attenuation (cases with stronger attenuation are included in several appendices). The
results indicate that attenuation plays a number of different roles even in the simple examined
cases. The attenuation between the noise sources and receivers produces the following effects on
the amplitudes and phases of the cross correlation function: (i) overall exponential amplitude
decay scaled by frequency and source-receiver propagation time, (ii) asymmetric amplitude
decay of the causal and anti-causal parts of the correlation function, in the presence of some non
uniformities, which increases with frequency, and (iii) phase shifts in the correlation function.
On the other hand, attenuation in the inter-station regions leads to phase shifts and non-
exponential symmetric decay of the correlation function amplitude. These effects are also present,
in a somewhat modified form, in numerical simulation results independent of the Fresnel
approximation and with non-isotropic source distributions (section 3). The results indicate that
the coherency method requires several modifications to be used for extracting information on
attenuation between pairs of stations. This issue and suggestions for observational studies are
discussed in section 4.
74
3.2. Theoretical results for isotropic source distribution
Figure 1. (a) Illustration of basic quantities used in equation (1) withstars representing noise
sources and the triangle denoting the receiver. (b) Geometry of 2D random point sources in a
ring and two receivers, one located at the origin and the other at (x, 0). The inter-station region
(small green circle) and outer region (white) are characterized by quality factors
in
Q and
out
Q ,
respectively.
The ambient seismic noise at location r
can be written for a single component of motion
using the scalar wavefield
) (
~
exp ) , (
'
1
, r dV
c
x
i r s
x
r u
V
, (1)
where x is the distance between the receiver and source ) , ( r s
at location r
(Figure 1a) and
' / 1 x is the geometric spreading factor. The term ) (
~
c is a complex phase velocity that
accounts for attenuation and is given by
) ( 2
sgn
1
~
Q
i c c , (2)
75
where Q(w) is the quality factor, c is a real phase velocity and sgn ensures that the
integral in (1) maintains the required Hermitian symmetry to have a real-valued wavefield in
time domain. The expected cross-correlation of wavefields at two receiver positions
1
r
and
2
r
is
V
u u
r r dV
c
x
c
x
i r s r s
x x
r u r u r r
) , (
~ ~
ω exp ω , ω , E
1
ω , ω , E ω ; , C E
*
*
2 1
*
2 1
2 1
(3)
where * denotes complex conjugate, r
and r
are source location vectors for receivers 1 and 2,
respectively, and x and x are the distances between source-receiver pairs 1 and 2, respectively.
We assume that sources at different locations are not correlated and each is a widely-sense
stationary random process
) ( δ ) ω ; ( B ω , ω , E
*
r r r r s r s
, (4)
where ) ω ; ( B r
is the source spectral density at location r
. Substituting this into (3) gives
V
u u
r dV
c
x
c
x
i r
x x
r r ) ( )]
~ ~
( ω exp[ ) ω ; ( B
1
ω ; , C E
*
2 1
2 1
. (5)
We simplify the problem to two dimensions (2D) by assuming that all sources and receivers
are at the surface (z=0 km), and convert r
to 2D polar coordinates ) , ( R . We put receiver 1 at
the origin (0,0) and assume that receiver 2 is at (x,0). We also assume in this section that the
random noise sources are uniformly distributed in a ring with radii R
min
and R
max
(Figure 1b).
With these assumptions, the expectation value of the cross-correlation between two receivers
becomes
d
c c
R
i R
R
dR
x x
R
R
R
u u
~
) (
~
ω exp ) ω ; , ( B
1
ω ; x , 0 C E
*
2
0
max
min
2 1
, (6)
where xR x R ) cos( 2
2 2
. For conditions corresponding to the Fresnel approximation
x R (source-receiver distance much greater than the receivers spacing), we have
x R R x R cos / cos 2 1
.
76
We consider situations where the region that includes the inter-station path (green circle in
Figure 1b) has a certain quality factor
in
Q , while the outer region may have a different quality
factor
out
Q . With such separation of Q values, the Fresnel approximation, and source spectral
density B , ; ω B( ω) R that is only a function of frequency, Eq. (6) becomes
d
c Q
x
i
c Q
x x R
i
c
x
i dR
in out
R
R
u u
2
cos sgn
2
cos cos 2 sgn cos
ω exp ) ω ( B ω ; x , 0 C E
2
0
max
min
2 1
(7)
In the following subsections we provide analytical results based on (7) for several basic spatial
distributions and properties of Q values.
3.2.1. Uniform attenuation
Assuming first that Q
in
= Q
out
= Q with no frequency-dependency, and collapsing Q back into
) (
~
c as in Eq. (2), the angular integral in (7) can be evaluated analytically to be
c
x
J d
c
x
i
~
2 ]
~
cos
ω exp[
0
2
0
. (8)
To convert (7) with uniform attenuation to time domain, we apply band-pass filter and inverse
Fourier transform
d t i
c
x
J
Qc
R
H dR
R
R
u u ) exp(
~
2
sgn
exp ) ( B
2
1
t ; x , 0 C
0
2
max
min
2 1
. (9)
The term
2
B( ω) B ( ) H is source spectral density that is narrow band-limited (which is
equivalent to applying a filter () H to the waveforms). We assume that
2
) ( H is a Gaussian
of the form (Weaver et al. 2009)
0 0
2 ~ ~
) ( A A H , (10)
where ) exp(
~
2 2
a A and 0 a .
77
The inverse Fourier transform of (9) is derived in Appendix A, assuming that the attenuation
is not overly strong so that Qc >> x , Q >R/(ac), and that the wave velocity is nearly
constant for a narrow band around
0
so
0
~ ~
c c . The last assumption of no dispersion
will be relaxed in section 2.4, along with the assumption on the frequency-independence of Q.
For these conditions, the expected value of the time domain cross-correlation can be written as
t t A
dt
d
t
t
rect
Q
Bt
t t A
t
t
rect
B
Q c
R
R R
t t
u u
0
2
0
2
0
0
min max
sin 2 *
) / ( 1
2
cos 4 *
) / ( 1
2
) (
exp t ; x , 0 C
2 1
, (11)
where
0
/ xc > 0 and a a t t A 2 / ) 4 /( exp
2 2
is a Gaussian envelope that is the
inverse Fourier transform of A
~
in (10). The two convolution terms in the bracket of (11) stem
from the inverse Fourier Transform of the Bessel function
0
Js with complex parameter in
equation (A3). The first convolution involves purely elastic waves, is symmetric in time (Figure
A1a) and appears in a similar form also in Weaver et al. (2009). In contrast, the second
convolution with the (1/Q) attenuation factor has an asymmetric term acting as two delta
functions shifted at ±τ with opposite signs (Figure A1b), and the wave function
t t A t t A t t A
dt
d
0 0 0 0
cos sin 2 sin 2
.
78
Figure 2. Schematic illustration of theoretical attenuation effects on cross-correlation functions.
(a) A cosine function with a Gaussian envelope centered at zero. (b) Exponential decay with
frequency scaling the overall amplitude due to attenuation between the receivers and far-field
sources. (c) Phase delays ofthe cross-correlation function in a case with attenuation (solid blue
line) compared to purely elastic case (dashed black line). (d) Asymmetry in the correlation
function in a case with uniform attenuation due to the asymmetric source-receiver propagation
distances. (e) Symmetric amplitude reduction of the correlation function governed by Struve
functions due to attenuation in theinter-station region. (f) Cross correlation function with the
combined theoretical effects of (b)-(e).
The narrow-band Gaussian filter leads to an elementary wave packet with a Gaussian
envelope (Figure 2a). The existence of attenuation modifies the cross-correlation function in
several ways. One basic effect is exponential amplitude decay with propagation time normalized
by the wave period (Figure 2b). The cosine term in the wave packet has the same phase as in the
elastic convolution and its amplitude increases with frequency, while the sine term in the wave
packet has a modified envelope t A and it shifts the phases in the cross correlation function
(Figure 2c). Carrying the convolution of the second term in (11) leads to asymmetric function in
time (Figure 2d) that becomes more prominent with increasing frequency. The discussed effects
can impact the amplitudes of cross-correlation functions as well as measurements of phase and
group velocities.
79
In Appendix A we consider the situation when the assumption Qc >> x is not satisfied,
so the 1
st
order Taylor expansion in (A2) is not a good approximation. This case corresponds to
conditions involving very strong attenuation, such as in fault zone environments with Q ≤ 15-30
(e.g., Ben-Zion et al. 2003; Peng et al. 2003), or when the central frequency of the band-pass
filter is very high. In section 3 we demonstrate further various effects of attenuation with
numerical simulations not limited to the assumptions associated with the analytical results.
3.2.2. Attenuation only within the inter-station path
out in
Q Q Q ,
In this case, Eq. (7) becomes,
d
c
x
i d
c
x
i dR
R
R
u u
~
cos
ω exp 2
~
cos
ω exp 2 ) ω ( B ω ; x , 0 C E
2 /
*
2 /
0
max
min
2 1
(12)
The inner integrals of Eq. (12) can be evaluated analytically to be
2 /
~
H
~ ~
cos
ω exp
2 /
~
H
~ ~
cos
ω exp
* 0 * 0 *
2 /
0
0 0
2 /
c
x
i
c
x
J d
c
x
i
c
x
i
c
x
J d
c
x
i
(13)
where
0
H is the zero
th
order Struve function.
We apply again band-pass filter and inverse Fourier transform to convert equation (12) to
time domain
d t i
c
x
i
c
x
i
c
x
J
c
x
J H dR
R
R
u u ) exp(
~
H
~
H
~ ~
) ( B
2
1
t ; x , 0 C
*
0 0
*
0 0
2
max
min
2 1
(14)
The inverse Fourier transform of (14) is given in Appendix B, and we make the same
assumptions on ) ω ( B ,
0
~ ~
c c , weak attenuation,narrow-band filter
2
) ( H as in section
2.1. Under these conditions, the expected value of time domain cross-correlation is
80
t t A
dt
d
t
t
rect t
Q
B
t t A
t
t
rect
B
R R
t t
u u
0
2
0
2
min max
sin 2 *
1
2 1
cos 4 *
1
2
t ; x , 0 C
2 1
(15)
where the phase delay time is same as in section 3.2.1. Eq. (15) shows that local attenuation in
a circle including the two receivers (Figure 1b) reduces the amplitudes of the signal
symmetrically (Figure 2e) and shifts the cross-correlation function phases (Figure 2c). These
effects can again impact the amplitudes of cross-correlation functions and measurements of
phase and group velocity. In Appendix B, we also discuss and illustrate with numerical
simulations the situation when the assumption Qc >> x is not satisfied, corresponding to
conditions when attenuation is very strong (small Q value) or the central frequency of the band-
pass filter is very high.
3.2.3. General Case
Combining results from sections 3.2.1 and 3.2.2, we can write the expected value of the cross
correlation function for the general case with different
out
Q and
in
Q as
d
c
x
i d
c
x
i
Qc
R
dR
R
R
u u
cos
ω exp 2
cos
ω exp 2
sgn
exp ω B ω ; x , 0 C E
2 /
2 /
0
max
min
2 1
(16)
where
out in
Q Q
i
2 1
2
sgn
1
and
in
Q
i
2
sgn
1
. The integrals in Eq. (16) can be evaluated to
be
2 / H
cos
ω exp
2 / H
cos
ω exp
0 0
2 /
0
0 0
2 /
c
x
i
c
x
J d
c
x
i
c
x
i
c
x
J d
c
x
i
(17)
where
0
J is the zero
th
order Bessel function and
0
H is the zero
th
order Struve function.
81
We make the same assumptions on ) ω ( B ,
0
~ ~
c c , weak attenuation and narrow-band
filter
2
) ( H as in sections 3.2.1 and 3.2.2. Based on a similar Taylor expansion as equation
(A2), we have
1 0 0 0
sgn
2 J
Q
i J
c
x
J
c
x
J
out
(18)
Likewise, we expand the Struve functions as in Eq. (B2) and get
2
H
1 1
H H
1 0 0
out in
Q Q
sgn i
c
x
c
x
(19)
Substituting Eq. (17-19) into Eq. (16), we apply band-pass filter and inverse Fourier
transform and get
]
sgn
exp[ *
2
H
1 1
sgn
sgn
2
* ) ( ) ( t ; x , 0 C
0
1
1 1 0
1
2
1
min max 2 1
c Q
R
Q Q
J
Q
i J
H R R B
out
t
t
out in out
t
u u
F F
F
(20)
Following similar steps as in the previous cases, we evaluate the inverse Fourier transforms
and obtain,
2
0
2
0 0
2
0
0
min max
1
2 1
* sin 2
1
2 1
* sin 2
1 1
cos 4 *
1
2
exp t ; x , 0 C
2 1
t
t
rect t
t t A
dt
d
Q
B
t
t
rect t
t t A
dt
d
Q Q
B t t A
t
t
rect
B
c Q
R
R R
t
out
t
out in
out
u u
(21)
The solution in Eq. (21) for the general case includes a superposition of the various effects
discussed in the previous sections (Figure 2f). If
out in
Q Q , Eq. (21) reduces to Eq. (11) of
section 3.2.1. If
out
Q , Eq. (21) reduces to Eq. (15) of section 3.2.2.
3.2.4. Dispersion and frequency-dependent attenuation
Here we also take into consideration dispersion and frequency dependence of Q. Assuming
variations of frequency within the narrow band of the filter, Eq. (20) becomes,
82
2
H
) (
1
) (
1
sgn
) (
sgn
2
*
sgn
exp * ) ( ) ( t ; x , 0 C
1 1 0
1
1
2
1
min max 2 1
x k
Q Q
x k x k J
Q
x k
i x k J
c Q
R
H R R B
out in out
t
out
t
u u
F
F F
. (22)
where c k / is the wave number.
Because the bandwidth of the filter is narrow, we linearize the wave number and attenuation
as a function of frequency. We use the same assumptions on ) ω ( B andnarrow-band filter
2
) ( H as in section 2.1. The inverse Fourier transform of Eq. (22) is derived in Appendix C,
leading to expected time domain cross-correlation function that is given in Eq. (C9). The result
can be written approximately as
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
min max
2 2
exp ) ( t ; x , 0 C
2 1
c
x
t sin t A
c
x
t sin t A
dt
d
p p x
c
x
t cos t A
c
x
t cos t A q q x
c
x
t sin t A
c
x
t sin t A
dt
d
xp
c
x
t cos t A
c
x
t cos t A xq
c
x
t cos t A
c
x
t cos t A
Qc
R
R R B
g g out in
g g out in
g g out
g g out
g g
u u
, (23)
where
g
is the group delay time and
0
/ c x is the phase delay time. Eq. (23) is comparable to
Eq. (21), but the incorporation of dispersion and frequency dependence of Q lead to a more
complicated form with constants
in
q ,
out
q ,
in
p , and
in
q that are defined in Appendix C near Eq.
(C5). The dispersion and frequency-dependant attenuation generally amplify the symmetric and
asymmetric effects of attenuation. If there is no dispersion or frequency dependant attenuation,
in
q and
out
q will be zero. As a result, the terms which contain
in
q and/or
out
q are caused entirely
by dispersion and frequency dependence of Q.
83
3.3. Numerical simulation results
Figure 3. Simulation results for isotropic source distribution without attenuation. Here f
0
denotes
the central frequency of a narrow-band filter in the range 0.03-0.15 Hz.
To validate the theoretical results of section 2, we perform numerical simulations based
directly on equations (1) and (2) for frequency-independent Q, without any of the other
assumptions or approximations made in the theory section. The simulation procedure includes
the following steps: (i) Use identical source spectral density (implemented as a FIR filter) for
each independent uncorrelated noise source. (ii) For each source-receiver pair apply travel time
delays to the filtered source spectrum based on a simple 1D linear dispersion curve of phase
velocity. (iii) For each source-receiver pair attenuate the amplitude of the spectrum following the
exponential decay with frequency, distance and Q in equations (1) and (2). (iv) For each receiver,
add the contributions in frequency domain from a distribution of sources in a ring. (v) Check for
Hermitian symmetry and convert the summed noise spectra to time domain.
84
In the results below, we use a ring of radii between 10 and 12 spherical degrees that contains
uniformly-distributed 30000 random noise sources. One station is located at the origin and the
other is 141.5 km away. We assume a linear dispersion curve defined with a maximum phase
velocity 3.1 km/s corresponding to 0.01 Hz and a minimum velocity 2.2 km/s corresponding to
0.18 Hz. For each station, the simulated continuous noise waveform is 10-day long with 2 Hz
sampling frequency and 1 Hz Nyquist frequency. All parameters can be easily modified to
accommodate different situations.
Figure 4. (a) Similar to Figure 3 for a case with homogeneous attenuation Q
in
=Q
out
=80. The red
horizontal bars mark portions of the correlation functions used to measure amplitudes of the
causal and acausal parts. (b) Amplitude asymmetry ratios A
R
versus frequency based on peak
amplitudes at the causal and acausal parts of the correlation functions. See text for more
explanation.
85
3.3.1. Isotropic source distribution with and without attenuation
The case for elastic wave propagation and isotropic source distribution provides a simple
reference model. Figure 3 shows simulated cross-correlation functions for the basic receiver
configuration of Figure 1. The amplitudes of the correlation function are generally symmetric
and reduce as the frequency increases, which is an inherent property of the zero
th
order Bessel
function. For the causal or anti-causal parts of the correlation function, the peak envelope
corresponds to the group delay time, which is different from the phase delay time due to the
surface wave dispersion. Because the two stations are not symmetric about the origin, there is a
slight asymmetry (less than 5% difference in relative amplitudes) in the cross-correlation
functions of Figure 3. This phenomenon is not predicted by our theoretical results since we use
the Fresnel approximation in the theory section which assumes source-receiver distance much
larger than the inter-station separation.
As a simple comparison example to the elastic case with isotropic sources, we consider a
situation with uniform attenuation associated with Q
in
=Q
out
=80. Figure 4 demonstrates that
homogenous attenuation produces two additional effects on the results: 1) exponential decay
with source-receiver distance and frequency of the overall amplitudes, and 2) larger asymmetry
of the cross-correlation functions at higher frequencies. We define the asymmetry of the peak
amplitudes in the causal and anti-causal parts as
) /( ) (
c ac c ac R
A A A A A , (24)
where
c
A and
ac
A are peak amplitudes measured from the causal and anti-causal parts,
respectively. The range of
R
A is [-1, 1]:
R
A equals 1 or -1 represents cross-correlation function
dominated by anti-causal or causal part, respectively, while 0
R
A corresponds to a symmetric
cross-correlation function. In the next section we show that the increasing asymmetry with
frequency can be separated from the asymmetry produced by non-isotropic source distribution.
86
Figure 5. Comparison between effects of nonisotropic source distributions and homogenous
attenuation.(a) Geometry for simulations with nonisotropic sources. The blue circle is a polar
plot of the azimuth noise source intensity with denoting theangle between the inter-station
direction and the maximum source intensity. Theangle gives the direction between an arbitrary
source and the inter-station line. (b) Amplitude asymmetry ratios A
R
for different cases versus
frequency. The values of A
R
are measured as discussed in the text and illustrated in Figure 4. The
max/min values in the legend denote ratios of the maximum noise intensity over the minimum
value for purely elastic cases with non isotropic sources distribution. The A
R
values do not change
in elastic cases systematically with frequency andwhen is 90 degrees the nonisotropic source
distribution produces an almost symmetric correlation function.The squares represent results of
simulation associated with isotropic source distribution and homogenous attenuation. The A
R
values in this case increase linearly with frequency.
3.3.2. Asymmetric amplitudes from non-isotropic source distribution and attenuation
Assuming the ambient noise sources have the same azimuthal distribution at all examined
frequencies, non-isotropic source distribution will generally produce the same asymmetry of
correlation function amplitudes for different frequencies. To simulate cases with non-isotropic
sources, we assume that the noise source intensity depends on the azimuth angle
) γ c os( 2 = B ,where θ is the azimuth angle, f is the angle of maximum intensity with
respect to the inter-station line, and γ < 2 is a positive number (Figure 5a). As in the isotropic
source distribution, the noise source locations are uniformly generated in the ring, but the
intensity varies with the azimuth angle.
Figure 5b displays simulation results for various elastic cases with non-isotropic source
intensities, along with a case having isotropic source distribution and attenuation. The red solid
line corresponds to a case where the maximum noise intensity is in the direction as inter-station
path and the maximum/minimum noise intensity ratio is 1.5. The green dashed line represents a
87
similar case with maximum/minimum noise intensity ratio of 1.2. The four dash-dot lines
represent cases with maximum/minimum noise intensity ratio of 1.5 and maximum intensity
direction rotated at different angles. The example with 90 degree rotation angle produces as
expected nearly symmetric results. In all the non-isotropic source cases, the amplitude
asymmetry fluctuates closely around a constant value (with the fluctuations related to numerical
simulation errors). In contrast, the case associated with isotropic source distribution and
attenuation (square symbols) shows clearly increasing amplitude asymmetry with frequency with
small fluctuations around a linear trend. If we can estimate the effect of non-isotropic source
distribution in a given region on the amplitude asymmetry, measurements of
R
A values versus
frequency may be used to estimate the averaged attenuation between the noise sources and
receivers pair.
3.3.3. Isotropic noise distribution and attenuation only within the inter-station region
Assuming the attenuation quality factor is finite only within the inter-station region, we study
the amplitudes of the resulting cross-correlation functions (Figure 6) in comparison with the
elastic case (Figure 3). As expected from the theoretical results, the inter-station attenuation
reduces the correlation function amplitudes symmetrically with larger amplitude reduction for
higher frequencies. This effect is potentially important for application of attenuation tomography.
If the contributions of other effects (e.g. asymmetric source distribution and attenuation between
sources and receivers) are carefully addressed, the inter-station attenuation can be recovered
using the frequency domain representation in equations (20) and (22).
88
Figure 6. (a) Simulation results for isotropic source distribution with inter-station attenuation
Q
in
=50. Here f
0
denotes the central frequency of a narrow-band filter in the range 0.03-0.15 Hz.
(b) comparison of peak amplitudes with the elastic case in Figure 3.
89
3.4. Discussion
We provide 2D analytical and numerical results on effects of attenuation on amplitudes and
phases of cross-correlation functions of seismic noise generated by sources in a far-field ring.
We analyze effects of attenuation in two complementary regions: a circle around the station pair
and the outer region between the stations and the noise sources. The incorporation of different
attenuation factors in the regions between and outside the stations is important for realistic
situations, and especially cases where receivers are in highly attenuating environments such as
fault zones (e.g. Ben-Zion et al. 2003; Lewis and Ben-Zion 2010) and volcanoes (e.g. Brenguier
et al. 2011) that are at considerable distance from the dominant noise sources. The analytical
derivations assume fully diffuse wavefield, the Fresnel approximation and relatively low
frequency and/or relatively large Q value. Cases with stronger attenuation are included in the
appendices, and more general situations including non-isotropic source distributions are
addressed in the numerical simulations section.
The analytical results indicate that the amplitude of noise cross-correlation functions decays
due to attenuation differently from what is assumed by Prieto et al. (2009) in several distinct
ways. The coherency expression of Prieto et al. (2009) is assumed to decay exponentially with
the inter-station distance. While our model results include exponential decay function, this is
associated with the propagation time between the sources and stations (sections 2.1 and 2.3). Our
analytical results indicate that the inter-station attenuation causes additional amplitude reduction,
which is non-exponential and related to properties of complex Struve functions in the frequency
domain (section 2.2). Our results indicate further that the attenuation also produces delay times
that can impact measurements of phase velocities (Figure 2). With isotropic far-field sources but
some non uniformity of properties between the sources and receivers or propagation distances,
the time shifts on the causal and anti-causal lags for uniform attenuation are asymmetric, while
the time shifts for inter-station attenuation are symmetric. In more realistic cases with different
regions for different frequencies and modes of noise sources (e.g. Stehly et al. 2006; Kimman
and Trampert 2010; Hillers et al. 2113) and heterogeneous velocity structures (e.g. Hauksson and
Shearer 2006; Tape et al. 2010; Allam and Ben-Zion 2012), more complex results are expected.
The amplitude reduction produced by attenuation outside the inter-station region can
generally be asymmetric because of differences in the regions between the sources and each
receiver. Such differences are represented in our work simply by different propagation distances.
90
If all properties between the sources and receivers are the same (propagation distances,
attenuation coefficients, directional properties of the sources, etc.), the attenuation in the region
outside the two stations leads to symmetric exponential amplitude decay. This would correspond
to results of Tsai (2011) on coherency with far-field sources and homogenous attenuation. Some
cases involving high frequencies and/or strong attenuation discussed in Appendices A and B
show that the correlation function amplitude reduces nonlinearly (increasingly fast) with
frequency, since the truncated Taylor expansion is not valid under these conditions.
As shown in Appendix D, for a case with uniform attenuation and isotropic source
distribution, the coherency between a pair of stations using the normalization of Weaver (2012)
is given from our results by
) (
~
/
0
) 2 (
c x J C E . (25)
The real part of equation (25) is the same as the real part of the coherency expression of Weaver
(2012), if we convert the quality factor Q to the attenuation coefficient α. However, the
imaginary parts of (25) and the coherency of Weaver (2012) have opposite signs. We note that
the imaginary part of the coherency produces asymmetry in the causal and anti-causal parts of
the noise cross-spectra. In the special case of uniform attenuation and isotropic sources, the
exponential decay term in equation (9) cancels out. However, for the more general case
discussed in section 2.3, the different attenuation in the inter-station region produces Struve and
Bessel functions in the normalization factor of the coherency which will make the form of the
coherency more complex.
The inter-station attenuation is a key parameter that is the target of attenuation tomography
studies. Equation (21) is a linear combination of first order Struve function and 0-order Bessel
function. This suggests that the quality factor Q
in
, which accounts for the ratio between these two
functions, can be retrieved together with an amplitude scaling constant that accounts for the
noise intensity and site effects. The asymmetric cross-correlation function due to different
propagation distances between the far-field sources and two receivers may be used to derive
information on the attenuation factor Q
out
outside the inter-station region. For example, in
California with dominant noise sources at the coast (e.g. Schulte-Pelkum et al. 2004; Hillers et al.
2013), using pairs of stations that are a few degrees off from the coast-parallel direction would
produce different propagation distances between the sources and receivers. An asymmetric
91
distribution of attenuation factors, as might exist for example for station pairs located on one side
of a fault zone, may also produce a similar effect on the cross-correlation functions.
If the inter-station attenuation is heterogeneous and asymmetric about the two stations (e.g. a
fault zone separates with unequal distance the two stations), the correlation function amplitude
reduction due to inter-station attenuation might be asymmetric (with variable asymmetry ratio).
This may be seen from the amplitude decay term c Q x
in
2 cos sgn exp in Eq. (7),
evaluated in the current context with
out
Q at θ and θ+π,. We note that the amplitude
reduction is most sensitive to the attenuation in the along-path directions (θ=0 or π) where the
decay is symmetric. In the perpendicular directions (θ=±π/2), the attenuation structure makes no
contribution to the amplitude of the correlation function. Such asymmetric amplitude decay due
to inter-station attenuation can be minimized if the inter-station distance is small and the Fresnel
approximation is satisfied. In addition, the inter-station formula may represent an averaged
attenuation effect of the inter-station region (primarily along path direction) as in travel time
tomography. These inferences should be substantiated with additional numerical experiments.
The numerical simulation method of section 3 is based on a simple ray theoretical summation
of the contributions of all noise sources with possible non-isotropic distribution to the waveforms.
The results demonstrate that different propagation distances in a uniform attenuation case
produce increasing asymmetry of the causal and anti-causal parts with increasing frequency
(section 3.1), while non-isotropic source distribution produces an asymmetric amplitude ratio
which fluctuates around a constant value for different frequencies (section 3.2). The results
suggest that these two effects on amplitude asymmetry can be separated. For a case dominated
by uniform inter-station attenuation (Q
out
very high) and non-isotropic source distribution, the
correlation function amplitude will also be asymmetric but the amplitude asymmetry ratio will
not vary with frequency. More complicated cases with heterogeneous attenuation between the
stations will lead to more complex results.
Effects of dispersion are examined briefly in the theory section 2.4 and the numerical
simulations. The results indicate that dispersion will amplify the predicted amplitude reduction
and amplitude asymmetry. The effects of dispersion should be accounted for in algorithms that
attempt to invert data for attenuation coefficients. We also note that it is easier to invert for
attenuation parameters from amplitude information in the frequency domain. This is because
92
fitting the combination of complex Bessel and Struve functions in frequency domain is more
straightforward, and not affected by the narrow-band filtering which is necessary to measure the
amplitude in time domain.
Finally, since the assumption on isotropic noise source distribution does not typically hold, it
may be possible to expand the non-isotropic noise source spectral density as a Fourier series
(Weaver et al., 2009; Tsai 2011) and evaluate the integrals for the cases examined in this paper
analytically for different Fourier terms. This and other important generalities of the simple cases
considered in this work (e.g. heterogeneous velocity structures) are left for future work.
Appendix 3A: Inverse Fourier transform of Eq. 9: Uniform Attenuation.
To simplify the inverse Fourier transform, we apply the convolution theorem
0
0
1
0
1
2
1
min max
~
*
sgn
exp * ) ( ) ( 2 t ; x , 0 C
2 1
c
x
J
Qc
R
H R R B
t t
u u F F F
(A1)
Let
Q
i
c
x
c
x
s
2
sgn
1
~
denote the complex phase delay time. Assuming that
Qc >> x
0
, we can apply a 1
st
order Taylor expansion to the Bessel function with complex
parameter (Yousif and Melka, 1997),
c
x
J
Qc
x
i
c
x
J s J
1 0 0
2
sgn
. (A2)
For narrow band analysis and large Q value, we assume small variations of phase velocity
c so the derivative of phase velocity can be ignored. Then we have
1
1 1
0
1
0
1
*
2
sgn
J
Q
i J s J
t
F F F F , (A3)
where
0
/ c x is the travel time between the two receivers for waves in a narrow frequency
band around
0
. The inverse Fourier transforms of the Bessel functions are
93
2
1
1
2
0
1
) / ( 1
2
2
2
1
) / ( 1
2
2
2
1
t
t
rect
t
i
J
t
t
rect
J
F
F
(A4)
The inverse Fourier transform of equation (10) is
t i t i t A
t A t i t A t i H
0 0
0 0
2
1
exp exp ) (
exp exp
F
(A5)
where ) 2 /( )) 4 /( exp(
2 2
a a t t A and hence t A t A .
Assuming 0
0
and applying the narrow band filter to (A1) gives
0
1
1 1
0
1
0
1
0
1
1 1
0
1
0
1
0
0
1
0
1
0
0
1
0
1
0
1
0
1
2
1
* *
2
sgn
*
~
* *
2
sgn
*
~
* *
~
* *
~
sgn
exp * *
0
0
0
0
0
0
0 0
Qc
R
i t J
Q
i s J A
Qc
R
i t J
Q
i s J A
Qc
R
i t s J A
Qc
R
i t s J A
Qc
R
s J H
D t t t
D t t t
D t t D t t
t t
F F F F
F F F F
F F F F
F F F
(A6)
The attenuation term is given by
iQ
t
Q Q
sgn
iQ
t
Q Q
D
D
2 2
sgn
2
2 2
sgn
2
sgn
0 0
0 0
1
0
1
1
0
1
F F
F F
(A7)
Substituting equations (A4), (A5) and (A7) into equation (A6), we have
0
1
0
1
2
1
sgn
exp * *
Qc
R
s J H
t t
F F F
94
2
0
2
0
2
0
2
0
2
0
2
2
0
0
) / ( 1
2
2
2
1
*
2
sin 2 *
2
1
2
2
2
1
*
2
cos 2
4
exp exp
t
t
rect
t
i
Qc a
R
t t A
Q
t i
t
t
rect
Qc a
R
t t A
Qc a
R
Qc
R
t t
D
t
(A8)
For the narrow-band filter
2
() H the width parameter a is very large. Therefore, when
attenuation is not very significant (i.e. Q >R/(ac)), we have 0 4 /
0
2 2 2 2
c Q a R and
0 ) 2 /(
0
2
Qc a R . Under these conditions, the time domain cross-correlation (expected value)
can be evaluated to be equation (12) in section 2.1.
-1 -0.5 0 0.5 1
1
1.5
2
2.5
3
3.5
4
4.5
5
t/ τ
-1 -0.5 0 0.5 1
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
t/ τ
(a) (b)
Figure A1. The inverse Fourier transform of Bessel function with complex parameter in
equation (A3). The time is scaled by
0
/ xc denoting the predicted travel time between the
two receivers for wave at frequency
0
.(a) The real part of (A3) yields a symmetric
function
2
) / ( 1 2 / t t rect for the elastic case. (b) The imaginary part gives an asymmetric
function
2
) / ( 1 2 / ) ( t Q t rect t scaled by attenuation assumed here to be Q = 50.
Figure A1 shows the inverse Fourier transform of the real and imaginary components of
s J
0
in equation (A3) for an example with Q = 50. The real part corresponding to the elastic
case produces a symmetric function in time (Figure A1a). The imaginary part associated with
uniform attenuation produces an asymmetric function in time (Figure A1b). The attenuation part
95
is small at lower frequency compared with the elastic part but becomes more significant as
frequency increases and/or Q decreases.
-1.5 -1 -0.5 0 0.5 1 1.5
-0.5
0
0.5
Frequency (Hz)
(a) Spectra representation of attenuation term: Q=50,station distance:50 km
-100 -80 -60 -40 -20 0 20 40 60 80 100
-5
0
5
Time (sec)
(b) Time domain delta function for frequencies up to 10 Hz.
-100 -80 -60 -40 -20 0 20 40 60 80 100
-0.2
0
0.2
Time (sec)
(c) Low frequency approximate delta function in time domain.
Accurate
Approximate
Figure A2. Comparison between approximate solution to equation (A2) associated with 1
st
order
Taylor expansion and full numerical solutionassuming uniform Q= 50, inter-station distance of
50 km, and phase velocity of 2 km/sec. (a) The approximate solution (red line) is accurate for the
assumed parameters only up to 0.7 Hz. (b) Thefull attenuation term in time domain computed
numerically from the blue curve in (a) with FFT. (c) The approximate attenuation term in time
domaincomputed from the red curve in (a) with FFT.
To supplement the analytical results associated with the assumption Qc >> x
0
we made
in the truncated Taylor expansion of equation (A2), we also simulate a case without truncation of
the Taylor series. This is done by computing s J
0
numerically with a Nyquist frequency of 10
Hz. and then using inverse Discrete Fourier transform to convert the results to time domain.
Figure A2 compares the imaginary part of s J
0
calculated numerically with the one computed
from the 1
st
order Taylor expansion in equation (A2) for Q=50 and 25 sec. For these
96
parameters, the 1st order Taylor expansion is a good approximation up to approximately 0.7 Hz
(Figure A2a). The amplitude of the pulses derived from the 1
st
order Taylor expansion is much
smaller (Figure A2b) than the full numerical solution (Figure A2c) because the Bessel function
with complex parameter increases very rapidly with frequency for f > 1 Hz. The results of Figure
A2 imply that cases with high frequency (or low Q values) should be examined with the full
numerical solution rather than the 1
st
order Taylor expansion.
Appendix 3B: Inverse Fourier transform of Eq. 14: Attenuation in the Inter-station Region.
We apply equation (A2) again to expand the Bessel function with complex parameter to the
1
st
order of Taylor expansion, assuming Qc >> x , and obtain
0
*
0 0
*
0 0
2
~ ~
J s J s J
c
x
J
c
x
J
(B1)
where the phase delay time is defined in equation (A3).
For the Struve function in equation (14) with an imaginary part of the parameter much less
than 1 and a real part much greater than the imaginary part, we can apply a 1
st
order Taylor
expansion,
2
H
2
sgn
H
~
H
2
H
2
sgn
H
~
H
1 0 * 0
1 0 0
Q
i
c
x
Q
i
c
x
(B2)
where
1
H is the first order Struve function. We therefore have
2
H
sgn
~
H
~
H
1 0
*
0
Q
i
c
x
c
x
. (B3)
Substituting equation (B1) and (B3) into equation (14) and applying the convolution theorem
lead to
t
t
u u
Q
J H R R B
2
H
sgn
2 * ) ( ) ( t ; x , 0 C
1 0
1
2
1
min max 2 1
F F
(B4)
97
Solving for the inverse Fourier transforms of the Bessel function and Struve function, we get
2
2
1
1
2
0
1
) / ( 1
2
2
2
1 2
H
) / ( 1
2
2
2
1
t
t
rect t
t
t
rect
J
F
F
(B5)
We now apply the narrow band filter to equation (B4), assuming
0
0 , and evaluate the
inverse Fourier transform in the neighborhoods of
0
and
0
,
2
H
sgn
2 * ) (
1 0
1
2
1
Q
J H
t
F F
0
~
1 0
1
0
1
2
H
sgn
2 *
~
Q
J A
t
F F
0
~
1 0
1
0
1
2
H
sgn
2 *
~
Q
J A
t
F F (B6)
For the attenuation coefficient term we have
iQ
t
Q Q
iQ
t
Q Q
D
D
'
1
0
1
'
1
0
1
0 0
0 0
sgn
sgn
sgn
sgn
F F
F F
(B7)
Substituting equation (B5) and (B7) into equation (B6) we get
2
2
0
'
2
0
1 0
1
2
1
1
2
2
2
1
* sin 2 *
1
2
2
2
1
* cos 4
2
H
sgn
2 * ) (
t
t
rect t
t t A
Q
t
t
t
rect
t t A
Q
J H
t t
D
t
t
F F
(B8)
98
The time domain cross-correlation (expected value) can be evaluated to be equation (15) in
section 2.2.
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
-0.012
-0.01
-0.008
-0.006
-0.004
-0.002
0
t/ τ
Figure B1. The inverse Fourier transform of the attenuation component (Struve function) in (B4).
The attenuation part
2
) / ( 1 2 / t Q t rect t reduces the amplitude of the signal
symmetricallyand is scaled by the inter-station attenuation Q =50.The time is scaled by
0
/ xc denoting the predicted travel time between the two receivers for wave at frequency
0
.
Figure B1 shows the attenuation (Struve function) component of the inverse Fourier
transforms in (B4) for a case with Q = 50. The elastic part is the same as in Figure A1a. The
attenuation part associated with the inter-station region produces a symmetric function in time
(Figure B1). To supplement the analytical results associated with the assumption Qc >> x
0
we made in the truncated Taylor expansion of equation (B2), we again simulate numerically with
Nyquist frequency of 10 Hz a case without truncation of the Taylor series. This is done by
computing s s
0
*
0
H H numerically and then using inverse Discrete Fourier transform to
convert the results to time domain. Figure B2 compares the s s
0
*
0
H H calculated
numerically with the one computed from the 1
st
order Taylor expansion in equation (B3) for
99
Q=50 and 25 sec. For these parameters, the 1st order Taylor expansion is a good
approximation up to 0.7 Hz (Figure B2a). The amplitude of the approximate delta functions
derived from 1
st
order Taylor expansion is much smaller than the full numerical solution (Figures
B2b and B2c) because s s
0
*
0
H H increases very rapidly with frequency when frequency
is greater than 1 Hz.
-1.5 -1 -0.5 0 0.5 1 1.5
-0.5
0
0.5
Frequency (Hz)
(a) Spectra representation of attenuation term: Q=50,station distance:50 km
-100 -80 -60 -40 -20 0 20 40 60 80 100
-6
-4
-2
0
2
Time (sec)
(b) Time domain delta function for frequencies up to 10 Hz.
-100 -80 -60 -40 -20 0 20 40 60 80 100
-0.15
-0.1
-0.05
0
0.05
Time (sec)
(c) Low frequency approximate delta function in time domain.
Accurate
Approximate
Figure B2. Comparison between approximate solution to equation (B3) associated with 1st order
Taylor expansion (RHS) and full numerical solution (LHS) of 0-order Struve function assuming
Q = 50 in the inter-station region, inter-station distance of 50 km, and phase velocity of 2 km/sec.
(a) The approximate solution (red line) is accurate for the assumed parameters only up to 0.7 Hz.
(b) The full attenuation term in time domain computed numerically from the blue curve in (a)
with FFT. (c) The approximate attenuation term in time domain computed from the red curve in
(a) with FFT.
100
Appendix 3C: Inverse Fourier transform of Eq. 22: dispersion and frequency-dependent
attenuation
We apply a Taylor expansion on k at
0
~ ,
d
dk
k k
0
0 0
) ( (C1)
We then solve the inverse Fourier transform of Bessel functions,
2
2
0
0
0 1
1
2
0
0
0 0
1
0
0 0 0
1
0
1
) / ( 1
2
) ( 2
2
1
1 exp
) / ( 1
2
2
2
1
1 exp
) (
0
0
g
g
g
g
g
g
g
g
t
t
t
rect t
c
v
it
c
x
J
t
t
rect
c
v
it
c
x
J
x
d
dk
x k J x k J
F
F
F F
(C2)
where
0
g
v is group velocity given by
0
0
dk
d
v
g
and
g
is group delay time given by
x
d
dk
v
x
g
g
0
0
.
Similarly we derive the inverse Fourier transform for the 1
st
order Struve function near
0
~
2
2
0
0
0 1
1
) / ( 1
2
2
2
1
1 exp
2
H
0
g
g
g
t
t
rect t
c
v
it
c
x
F . (C3)
We define
out
out
Q
k
and expand it at
0
as 1
st
order truncated Taylor series
101
d
d
out
out out
0
0 0
.
(C4)
We then define two related constants
out
out
out
out out
out
out
out
out
out
p
Q
k
d
d
q
Q
k
d
dQ
Q
d
dk
d
d
p
0
0
0 0
0 0
0
2
0
0
0
0
0
(C5)
Similarly, we define
in
in
Q
k
, expand it as 1
st
order Taylor series, and define constants
in
p and
in
q as in equation (C4) and (C5) by replacing
out
Q with
in
Q . With these definitions we
have
i t p t q x x sgn
c Q
x sgn
D out D out
out
/
' 1
0
1
0
0
F (C6)
For the narrow-band filtering at
0
~ , we note that the wave number is an odd function
0 0
k k (the velocities and attenuation are even functions of frequency). Therefore, for
equations (C1) to (C4) we only need to replace
0
by
0
and a similar version of (C6) can be
written as
i t p t q x x sgn
c Q
x sgn
D out D out
out
/
' 1
0
1
0
0
F (C7)
We now narrow-band filter the exponential decay part and convert it back to time domain
t i t A
c Q
R
t i t A
c Q
R
Rq iRp t A Rq iRp t A
c Q
R
H
out out
out out D t out out D t
out
t
0
0
0
0
0
0
0
1
0
1
1
2
1
exp exp exp exp
exp *
~
exp *
~
sgn
exp *
F F
F F
(C8)
102
where we make the same assumptions about the width parameter a (i.e. Q >R/(ac)) as in
Appendix A. Combining equations (C2) to (C8), we can solve the inverse Fourier transforms in
equation (22),
t h t sin t cos t A t h t cos t sin t A t xp
t h t sin t sin t A t h t cos t cos t A xq
t h t g t sin t cos t A t h t g t cos t sin t A t xp
t h t g t sin t sin t A t h t g t cos t cos t A xq
t f t sin t sin t A t f t cos t cos t A
Qc
R
R R B
t t t D in
t t in
t t t D out
t t out
t t
u u
* * * 2
) ( * * 2
* * * 2
* * 2
) ( * ) ( 4 * 4
exp ) ( t ; x , 0 C
0 0
'
0 0
0 0
'
0 0
0 0
0
0
min max 2 1
(C9)
where
0
0
0
1
g
v
c
and the three functions t f , t g and t h are defined as
g D g D
g
g
g
g
g D g D
g
g
g
g
g D g D
g
g
g
t t
t
t
t
rect t
t h
t t
t
t
t
rect t
t g
t t
t
t
rect
t f
2
) / ( 1
2
1
2
) / ( 1
2
) (
1
2
1
) / ( 1
2
1
2
2
2
2
2
(C10)
In equation (C10), we approximate the three functions with Dirac delta functions which are
the RHS of the approximately equal sign. With this approximation, we substitute (C10) into (C9)
and arrive at equation (24) in section 2.4.
Appendix 3D: Derivation of normalized cross spectra (coherency) for uniform attenuation
and isotropic source distribution
The normalized cross spectra of noise correlation functions is referred to as coherency (e.g.
Prieto et al. 2009). There are several definitions of coherency that differ by the normalization
103
factor. Here we derive the noise coherency for uniform attenuation and isotropic source
distribution using the definition of Weaver (2012). In our notation this is expressed as
ω ; 0 , 0 C E
ω ; x , 0 C E
ω
1 1
2 1
) 2 (
u u
u u
C E
, (D1)
where ω ; x , 0 C E
2 1
u u
is the expected cross-spectrum given by equation (7) and simplified for
uniform attenuation in equation (9). The normalization factor ω ; 0 , 0 C E
1 1
u u
is the power
spectral density (autocorrelation in time domain) of the noise recorded by the station at the origin,
dR
Qc
R
R
R
u u
sgn
exp B 2 ω ; 0 , 0 C E
max
min
1 1
, (D2)
and the integral is evaluated using the Mean Value Theorem. Combining equations (7), (9) and
(D2), we evaluate the coherency in (D1) as,
c
x
J C E
~
0
) 2 (
. (D3)
104
4. Extracting Seismic Attenuation Coefficients from Cross-Correlations of Ambient Noise
at Linear Triplets of Stations
(Liu et al., 2015)
4.1. Introduction
The ambient seismic noise contains rich information on structural properties between the
recording stations (e.g., Weaver and Lobkis, 2004; Campillo, 2006). Empirical Green’s functions
for surface waves can be retrieved by cross-correlating continuous noise data at different sites
(Shapiro and Campillo, 2004). The phase and amplitude information contained in the empirical
Green’s function have various applications in seismology (e.g., Campillo et al., 2011; Lin and
Ritzwoller, 2011; Boué et al., 2014). For successful retrieval of empirical Green’s function, the
noise wave field should be diffuse and the noise sources should be isotropic and uncorrelated. In
practical cases the noise field is not fully diffuse and often has strong directionality, leading to
asymmetric cross-correlations and ambiguity in the derivation and interpretation of results (e.g.,
Weaver 2011). Despite these difficulties, there have been many successful applications of
retrieving information on seismic velocities from cross-correlation of ambient noise (e.g.,
Shapiro et al., 2005; Zhan et al. 2014; Zigone et al., 2015). Group velocities can be derived with
Multiple Filter Analysis (MFA) techniques (e.g., Dziewonski et al., 1969; Herrmann 1973;
Pedersen et al., 2003) while phase velocities can be measured from the phases of narrow-band
filtered noise cross-correlation functions (e.g., Yao et al., 2006; Lin et al., 2008; Boschi et al.,
2013).
Getting amplitude information and inferring attenuation coefficients (Q values) from noise
cross-correlation data is considerably more difficult. Prieto et al. (2009) derived inter-station
attenuation by fitting ensemble averaged coherency (Aki 1957) as a function of inter-station
distance with a damped Bessel function. However, the method ignores the attenuation between
the sources and stations and involves additional assumptions that limit its utility (e.g. Tsai 2011;
Weaver 2013; Liu and Ben-Zion 2013). Weaver (2013) proposed a method for deriving
attenuation from linear array of at least five stations based on the radiative transfer equation and
stationary phase approximation (Snieder 2004). Applying the radiative transfer equation to both
sides (causal and anticausal) of synthetic linear arrays, Weaver (2013) estimated the attenuation
105
coefficients together with site amplification factors and confidence intervals from synthetic data.
Lin et al (2012) derived site amplification factors by applying the eikonal equation with damping
factors to earthquake data recorded by a 2D array. To extend that study to noise cross-correlation
data may require examination of amplitudes across the entire 2D array, accounting for anisotropy
of the noise sources and attenuation between the numerous sources and receivers (Galetti and
Curtis, 2012).
In the current study, we develop an inversion algorithm based on previous theoretical results
of Liu and Ben-Zion (2013) on effects of attenuation on noise cross-correlation between pairs of
stations in a model with separate inter-station Q value from that of the background attenuation.
We apply the inversion formula to linear arrays of three stations (triplet), accounting for inter-
station attenuation values, background attenuation and site amplification factors. We perform
numerical simulations with noise cross-correlation data described in terms of matrix
multiplication and data covariance matrix. Synthetic noise cross-correlation data are generated
from thousands of random noise sources, each of which is a sampling of the noise source
probability space. Applying the inversion procedure to synthetic data shows that the three Q
values between three stations can be reliably estimated for non-isotropic sources when the
maximum noise intensity is oriented along the inter-station direction. In contrast, only the Q
value between the end stations can be reliably recovered for isotropic source distributions while
the other two Q values are subjected to trade-offs. Applying the procedure to observed data, we
analyze two linear arrays having different length scales. The first has 3 stations of the southern
California regional network separated by about 35 km, while the second consists of 6 stations
across the San Jacinto fault zone with inter-station separation of about 25 m. For the former array
we use a frequency band between 0.2-0.5 Hz and obtain Q
R
values associated with Rayleigh
waves around 25. For the latter array we use several frequency bands between 15-25 Hz and
obtain Q
R
values of about 6-30 with the lower values characterizing the fault damage zone.
Bootstrapping the results with data of single days yields for most station pairs small confidence
intervals of Q
-1
values indicative of reliable estimations.
106
4.2. Theory
Figure 1. (a) Illustration of the 2D geometry for random point noise sources in a ring (gray),
attenuation structure (blue, orange and green areas described below), and station location (black
triangles). The two stations are symmetric about the center of the geometry. The interstation
region (green) and the west and east outer regions are characterized by quality factor Q
in
, Q
outW
and Q
outE
, respectively. (b) Under the stationary phase assumption, only the noise wavefields
propagating near the inter-station direction (end-fire lobes)have constructive interference on the
cross-correlation.
We consider a pair of stations in a solid with 3 different attenuation coefficients
characterizing the medium between the stations and the bounding regions to the left (west) and
right (east) of the stations (Fig. 1a). Noise sources are assumed in the far field (denoted by the
ring). Under stationary phase approximation (Snieder 2004), the cross-correlation amplitude is
affected mostly by noise propagation along the inter-station direction and noise sources in other
directions have little contribution (Fig. 1b). This has been generally referred to as end-fire lobes
in studies of noise correlations without attenuation (e.g., Roux et al., 2004, Goué dard et al .,
2008). Below we combine this well-known phenomenon with the expression of noise cross-
spectrum in dissipative medium (Liu and Ben-Zion, 2013) to develop an inversion algorithm.
Based on Liu and Ben-Zion (2013) and derivations in Appendix A, under the assumption that
interstation distance x is greater than the wavelength λ (2πx>>λ), the expected cross-spectrum
between two stations can be approximated as
107
. (1)
where R is the distance between far-field sources and origin, ω is angular frequency and c(ω) is
phase velocity. Here β
1
and β
2
are real-value frequency-dependent site amplification factors at
stations 1 and 2, respectively. The site amplification factors are assumed real; otherwise the
phase velocity dispersion derived from cross-correlation will not be accurate. The two complex
numbers ε
w
and ε
E
accounting for rotation of phase angle due to attenuation are defined below Eq.
(A2) of Appendix A.
Figure 2.Schematic diagram for inverting attenuation from a linear array of three stations. Black
triangles show the location of stations 1, 2 and 3 with interstation distances x
1
, x
2
, x
3
and
interstation quality factors Q
1
, Q
2
, Q
3
, respectively. Here Q
3
correspond to the apparent
attenuation between stations 1 and 3. The distance between station 1 and 3, x
3
, is the sum of x
1
and x
2
. The background attenuation can be canceled with this type of array as explained below
Eq. (2).
The causal part (first line) of Eq. (1) contains noise source spectral density B(ω), amplitude
correction factor
W
x c i 2 , an exponential decay term (attenuation term) and a phase
term (the imaginary complex exponential). An intuitive interpretation of the attenuation term in
Eq. (1) is that the attenuation of the causal part of the cross-correlation function is only relevant
108
to noise wave field propagating from the west along (nearly) the inter-station direction, and vice
versa. As a result, the exponential decay of the causal part is determined by the Q
outW
and Q
in
values between the western noise sources and two receivers along the end-fire lobe direction (Fig.
1b). The phase velocity can be extracted from the phase term multiplied by i (giving π/4 phase
shift) with non-uniqueness because the phase angle is wrapped on an interval between (-π, π).
A straightforward linear inversion can be formed by taking the ratios of amplitude cross-
spectra functions in a given frequency band of different station pairs of a linear array. The
simplest situation involves 3 stations on a line (Fig. 2), because the background attenuation term
and noise source power spectral density B(ω) can be canceled in this situation. The noise
waveforms at three stations are denoted u
1
, u
2
and u
3
, respectively. The amplitudes may be
measured on either the causal or anti-causal parts of the amplitude cross-spectra, and then
corrected by dividing with the amplitude correction term
W
x c 2 (We assume ε
w
≈1 in
realistic cases). Let
2 1
C
u u
,
3 2
C
u u
,
3 1
C
u u
be the corrected amplitude decay curves of the
causal cross-correlation functions of three station pairs. We form three linear least-square
equations for the three attenuation factors
(2)
where the relevant distances x and quality factors Q are defined in Fig. 2. Here γ
2/1
, γ
3/2
, γ
3/1
account for the natural logarithms of the ratios of corresponding site amplification factors (
values in Eq. 1) and they are assumed to be constant in the frequency band used for inversion. Eq.
(2) still holds if the site amplification factors at three stations are correlated with each other and
can be removed by taking the ratio of them in Eq. (2).
A Q value and a ratio of corresponding site amplification factors can be inverted from each
expression of Eq. (2) with simple linear least-square inversion in a chosen frequency band. The
frequency band should be wide enough to allow for sufficient amplitude decay that exceeds
109
random amplitude fluctuation. The inverted Q value only reflects the average amplitude decay
within that frequency band. A similar version to Eq. (2) defined for a noise wave field
propagating from the west (assumed to be the causal part of the cross-correlation function), can
be derived for noise wave field propagating from the east (anti-causal part of cross-correlation).
We note that if the attenuation is strong enough, it may shift the phase velocity retrieved from the
cross-correlation function. This can be observed from the complex coefficients
W
and
E
(Eq. 1)
which contain attenuation quality factors and rotate the phase angle in the causal and anti-causal
parts. However, the attenuation has to be very strong (e.g., Q ~ 2-3) to cause a noticeable phase
shift in the cross-correlation function. Therefore it is safe to assume ε
w
= ε
E
=1 in most practical
cases.
4.3. Numerical simulations of noise waveforms and inversion of Q values
The discrete version of noise spectrum produced by far-field noise sources at k receivers (Eq.
1 in Liu and Ben-Zion 2013) is
, (3a)
or
PS U , (3b)
where P is the propagation matrix, U is a vector that contains the spectrum at each station and S
is a random vector of noise sources at different azimuth angle θ. In (3a), θ
1
to θ
n
are equally
spaced azimuth angles (used to approximate an integral) of n assumed sources with respect to the
center of the array. The source-receiver distance from source m to receiver l is
l m
x . The
inverse complex phase velocity is defined as c Q i c ) ( 2 / sgn 1
~
1 where Q is the
(potentially space-dependent) attenuation quality factor. The geometric spreading factor is
approximated by 1/ R where R is the average source-receiver distance (Fig. 1).
110
Assuming a stationary ambient noise process, the covariance matrix of the random spectrum
data from all receivers is,
CovU = PCov S
( )
P
+
, (4)
where Cov(S) is the covariance matrix of the noise sources and the dagger operator denotes
conjugate transpose. For equipartitioned noise wavefield, any two different noise sources are
uncorrelated, which implies that Cov(S) is a diagonal matrix. For isotropic source distribution,
Cov(S) is the identity matrix. For the data covariance matrix at frequency ω, each non-diagonal
entry represents the cross-spectrum of two stations at that frequency while each diagonal entry
represents auto-spectrum. Synthetic data can be generated from Eqs. (3) and (4) by either (i)
repeatedly compute data covariance matrix
UU from random noise source vector S and then
average the results, or (ii) compute CovU directly from a known source covariance matrix
Cov(S). The former can simulate finite length of data while the latter is directly computing the
ensemble average of cross-spectrum. The propagation matrix P can be pre-computed to save
time.
30’
117
o
W
30’
116
o
W
30’
33
o
N
20’
40’
34
o
N
20’
40’
35
o
N
Longitude
1 2 3
Q1 Q2
Latitude (a)
(b)
Figure 3. (a) Forward simulation geometry with heterogeneous attenuation. The black triangles
show the station locations. The region between stations 1 and 2 has Q
1
=80 and x
1
=24.2 km; the
region between stations 2 and 3 has Q
2
=30 and x
2
=32.3 km. Their weighted average attenuation
between stations 1 and 3 isQ
3
=41 with distance x
3
=56.5 km. The background Q (white area) is
500.(b) Along-path attenuation and grid cells between two stations (black triangle). Each grid
cell is defined with an attenuation Q value.Exponential decay of amplitude is computed by
taking into account each grid cell along the ray path.
111
The forward wave propagation calculations are based on ray theory. We first discretize the
region of interest. In the following we use grid spacing of 4.04 km and overall grid size of 202
km x 222 km (Fig. 3a). We focus on Q
R
values associated with attenuation of Rayleigh waves,
although for brevity we sometime just use Q to denote Q
R
. The region between stations 1 and 2 is
defined with attenuation Q
1
= 80 and distance x
1
= 24.2 km while the rectangle between stations
2 and 3 is given Q
2
= 30 and x
2
= 32.3 km. The averaged Q
3
value between stations 1 and 3 is 41
with distance x
3
= 56.5 km. The background attenuation is assumed to be Q
out
= 500. We apply a
UTM projection to convert spherical coordinates to local Cartesian coordinates for the region.
For a ray path between a source-receiver pair, only the portion within the stations is discretized
to account for heterogeneous attenuation by adding up the contribution of each grid cell (Fig. 3b).
The portion of ray path outside the local grid is assumed to travel in homogeneous dissipative
media with time delay and amplitude loss computed along the great circle path.
112
Figure 4.Geometry and results for the synthetic example. (a) Non-isotropic source distribution
around 2 stations (black triangles): maximum intensity rotation angle ϕ and max/min ratio of
noise intensity distribution.(b) Phase velocity inversion for three station pairs. Black line: model
phase velocity; dashed line: inverted phase velocity. There isgood agreement between model and
inverted results. (c) Cross-correlation for stations 1 and 3 before windowing (blue) and the causal
part of the cross-correlation after applying a Turkey window (red) based on lower and upper
bounds of phase velocity (1-6 km/s). (d) Corrected amplitude measurement (log scale) as a
113
function of frequency for the causal parts of the 3 cross-correlation functions in a non-isotropic
source case with max/min=1.5 and ϕ =180 (maximum intensity in the West).The colors indicate
the station (see legend on the panel). (e) Relative decay curveas a function of frequency for the
non-isotropic source distribution defined in panel d. The curves are differences between the
logarithmic amplitude curves in panel d.The inverted Qvalues indicated in the panel are very
close to the model values in Fig. 3b. Best-fitting curves are plotted as solid lines with the same
color as data. (f) Same as (e) for an isotropic source distribution. The invertedQ
3
is very close to
the assumed value, while the inverted Q
2
and Q
1
are less accurate because oftradeoffs.
We simulate 20000 noise sources assumed to be uncorrelated stochastic sources located 11
degrees (distance in spherical degrees) away from the center of the grid. We consider both non-
isotropic sources with maximum intensity direction in the west along the inter-station direction
(Fig. 4a, max/min=1.5) and isotropic source distributions. The estimated cross-spectra and auto-
spectra (power spectra) are computed by 8000 iterations. The frequency band is 0-0.86 Hz with
sample interval of 0.005 Hz. The computation takes about 11 hours on a single CPU.
As a prior information for the attenuation inversion, phase velocity curves for the three
stations in Fig. 3a are first inverted from real part of the cross-spectra. Fig. 4b shows the inverted
phase velocity for all three station-pairs for the case with non-isotropic source distribution. The
isotropic case yields very similar phase velocity results. The inverted phase velocity curves are
very similar to the model phase velocity except that the map projection produces small errors in
distance calculations. The causal part of the cross-correlation (red waveform in Fig. 4c) is
isolated by applying a Turkey window defined by velocity bounds between 1 km/s and 6 km/s on
the cross-correlation (blue waveform in Fig. 4c). Amplitude decay curves are measured from the
amplitude cross-spectra of the causal correlation function, which represent noise wave
propagating from west to east. Then the decay curves are corrected according to Eq. (1) by
removing the amplitude correction term (Fig. 4d, in logarithmic scale).
Attenuation Q values are inverted according to Eq. (2). For non-isotropic noise source
distribution, Fig. 4e shows results for the LHS of Eq. (2) as three relative decay curves giving the
differences between the logarithms of the measured-corrected amplitudes of the three cross-
correlation functions. Each expression in Eq. (2) is solved with linear least-square inversion for a
Q value in the frequency range from 0.2 to 0.5 Hz. The inverted Q values for this case agree well
with the forward model parameters in Fig. 3a. For isotropic source distribution, the averaged
attenuation Q
3
between stations 1 and 3 is recovered with good precision (Fig. 4f). However,
there are some tradeoffs between the inverted Q
1
and Q
2
values. One explanation is that the
114
cross-correlation function between stations 1 and 3, t
u u
3 1
C
, should be symmetric as predicted
by the stationary phase approximation under symmetric background attenuation. However, the
actual simulated data have slight asymmetry since the noise waveforms propagating outside the
end-fire lobes do not cancel out completely because of the asymmetric inter-station attenuation
(Q
1
, Q
2
) between stations 1 and 3. An inaccurate t
u u
3 1
C
produces errors in Q
1
and Q
2
derived
using Eq. (2).
The synthetic tests show that the best scenario for this inversion method is non-isotropic
noise source distribution (at least for the parameters used in this section), with maximum
intensity aligned with the inter-station direction and smoothly varying noise intensity. For
isotropic source distribution, only the inverted average attenuation between stations 1 and 3 has
good precision. In addition, because the aperture angle for the end-fire lobes decreases as
interstation distance increases (Snieder 2004; Roux et al., 2004), the spacing between stations in
the triplet should not be significantly different. Although the synthetic test is set up for a
relatively low frequency range between 0.2-0.5 Hz (and corresponding relatively large station
separation), the conclusions hold for higher frequencies and densely spaced arrays if 2πx>>λ is
satisfied.
115
4.4. Inverting Q values for a regional array
Figure 5.A map for a triplet of stations, CHF-SBB2-LMR2 (black triangles), of the regional
southern California network. Black lines indicate fault traces. The green and yellow show low
and high elevation respectively.
Using a similar technique, we derive Q
R
values for the paths between a triplet of stations of
the southern California seismic network. The three used stations, CHF-SBB2-LMR2, cross the
Mojave section of the San Andreas fault (Fig. 5). The analysis is based on vertical component
seismograms recorded from Julian day 50 to day 350 of year 2014. The chosen three stations are
located in relatively simple crust outside large sedimentary basins (although station CHF is in a
mountainous area). The inter-station distances for CHF-SBB2, CHF-LMR2 and SBB2-LMR2
are 43.6 km, 73.3 km and 29.7 km, respectively. Before preprocessing step, the instrument
responses are removed, the raw waveforms are highpass filtered with a lower cutoff frequency of
0.02 Hz and the sampling rate is decimated to 4 Hz. The pre-processing of the noise data
includes glitch correction and removal of segments contaminated by earthquakes (see Zigone et
al. 2015 for more details). However, in contrast to typical noise-based imaging studies no
frequency whitening is applied to avoid amplitude distortion (e.g. Weemstra et al. 2014). The
computed daily cross-correlations are stacked to obtain one combined cross-correlation for each
station pair (Fig. 6a).
116
As a prior information for amplitude correction and attenuation inversion, phase velocity
dispersion curves are computed from the symmetric component of the cross-spectra, which is the
real part of Eq. (1). As a post-correlation whitening step for phase velocity measurements, we
normalize the stacked cross-spectra by the power spectra functions of both receivers (e.g. Aki
1957; Prieto et al. 2009) before phase velocity measurements. We adopt two methods for phase
velocity analysis: (i) estimate and unwrap phase angle data measured on narrow-band filtered
cross-correlation in time domain (e.g. Lin et al. 2008), and (ii) estimate phase angle data in
frequency domain from the inverse Fourier transform of the causal part of cross-correlation
(phase term in Eq. 1). These two methods produce consistent results. The measured phase
velocity curves are shown in Fig. 6b; the dispersion curves intersect at two points since the phase
velocity on the longest path (CHF-LMR2) is the weighted average of those on the two shorter
paths (CHF-SBB2 and SBB2-LMR2) in the triplet.
In a second analysis step, the cross-spectral amplitude is measured by separating the causal
and anti-causal parts of un-whitened complex cross-spectra according to Eq. (1). There are
generally two ways of separating the two opposite propagating terms: (i) frequency domain
Hilbert Transform of the bandpass filtered complex cross-spectra based on Eq. (1), and (ii) time
domain windowing on the causal part cross-correlation function followed by Fourier transform to
frequency domain. The former has better resolution in the frequency domain but does not
remove the correlation coda and requires that the cross-correlation function only contain the
fundamental mode Rayleigh wave. The latter assumes that signals other than the fundamental
Rayleigh wave travel at different speeds and therefore the fundamental Rayleigh wave mode can
be isolated by windowing at predicted time interval (although different modes may overlap and
can be hard to separate). However, the energy of the traveling wave in the causal part may suffer
from leakage at the anti-causal part due to fast decaying ambient noise energy with frequency,
small inter-station distances and P wave incidence below the array. Therefore, we use the former
method for triplets with much stronger surface wave packets than body waves and the latter
method for cases with weaker surface wave packets that need to be separated from body waves.
117
Figure 6.Cross-correlation functions, phase velocity curves, amplitude measurements and
inverse Q distribution for CHF-SBB2-LMR2 (see station locations in Fig. 5). (a) Cross-
correlations for CHF-SBB2-LMR2. Distances between each pair of stations are shown in the
titles. (b) Phase velocity dispersion curves. (c) Corrected amplitude decay curves for the three
correlation functions. The frequency range is 0.2-0.36 Hz and is determined in relation to the
inter-station distance. (d) Relative amplitude decay curve and inverted Q
R
values from a
randomresampling of 300 days in year 2012.Best-fitting curves are plotted as solid lines with the
118
same color as data. (e) Distribution of inverse Q values for the 3 station pairs based on
bootstrapping statistics (1/Q
1
=0.049±0.011; 1/Q
3
=0.040±0.006; 1/Q
2
=0.031±0.013).
For amplitude measurements with data of the regional triplet, we apply Wiener filter to
reduce the noise level in the cross-correlations. We first define a signal window by phase
velocity bounds between 1 km/s and 6 km/s. The noise window starts 5 sec after the signal
window and has a total duration of 50 sec. The Wiener filter is a minimum mean-square
estimator of the surface wave packet signal. Because of the strong ocean microseism peak
around 0.16 Hz, the amplitude of cross-spectra above 0.2 Hz decays very fast and the cross-
correlation is dominated by 7 sec waves. Therefore we deconvolve the autocorrelation of the two
stations from the cross-correlation following the approach of Aki (1957) before Wiener filtering
and convolve them back after the Wiener filtering.
The corrected amplitude decay curves for the three station pairs, based on the stacked cross-
spectra for random resampling of 300 days, are shown in Fig. 6c. We smooth the amplitude
decay curves with a running average smoothing parameter of 0.05 Hz. The frequency band used
for the linear least-square inversion is between 0.2 Hz and 0.36 Hz. The relative amplitude decay
curves and inverted Q values are shown in Fig. 6d. To estimate the ensemble mean and
confidence interval of the inverted Q values, we bootstrap the 300 days and get distributions of
three inverse Q values (Fig. 6e). The inverted results suggest that CHF-SBB2 has slightly
stronger attenuation (Q
1
≈21, 1/Q
1
=0.049±0.011) than SBB2-LMR2 (Q
2
≈32, 1/Q
2
=0.031±0.013),
possibly due to damaged rocks around the San Andreas Fault, but the difference is not highly
significant statistically. The average attenuation between CHF and LMR2 is Q
3
≈25
(1/Q
3
=0.040±0.006) associated with a narrow confidence interval.
119
4.5. Inverting Q values for a fault zone array
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
JF00 Los
Angeles
Figure 7.A map for the JF array with 9 broadband stations (black triangles) across the San
Jacinto Fault Zone.The black lines mark surface fault traces. The inversion for Q values is done
for data by the 6 stations in the red square forming 4 station triplets.
Here we derive Q
R
values for the paths between 6 stations of a dense linear array that crosses
the southern San Jacinto Fault zone at site JF (Fig. 7). As in the preprocessing stage for the three
stations of the regional network, we avoid pre-whitening, perform glitch corrections and remove
segments affected by earthquakes. The cross-correlations are computed on the radial RR
component instead of the vertical ZZ component in order to minimize body P wave phases (e.g.
Hillers et al., 2013) which could affect the amplitude measurements. The phase velocity
dispersion curves are derived from the stacked cross-correlations as discussed in the previous
section. We compute Signal-to-Noise Ratio (SNR) for daily cross-correlations from a set of 80
days for the year 2012 and select 34 days with highest SNRs for all 6 stations. The amplitude
decay curves are measured as in the previous section, using either frequency domain Hilbert
transform or time domain windowing, except that Wiener filter is not applied to triplets with very
high SNR (e.g. JFS2-JFS1-JF00). We assume that each one of the 34 days is an independent
observation of the amplitude cross-spectra and that the ambient noise field is relatively stationary
within each day (so that the cross-correlation can converge). A bootstrapping technique is used
to estimate empirical mean value and confidence interval of the attenuation Q
-1
by randomly
resampling the 34 days of cross-correlations.
120
Figure 8.Cross-correlation functions, phase velocity curves, amplitude measurements and
inverse Q distribution for JFS2-JFS1-JF00 (station locations are in Fig. 7). (a) Cross-correlations
for JFS2-JFS1-JF00. Distances between each pair of stations are shown in the titles. (b) Phase
velocity dispersion curves. (c) Corrected amplitude decay curves for three correlation functions.
The frequency range is 20-25 Hz related to the small inter-station distance. (d) Relative
amplitude decay curve and inverted Q
R
values from a random resampling of 34 days.Best-fitting
curves are plotted as solid lines with the same color as data.(e) Distribution of inverse Q values
121
for three station pairs based on bootstrapping statistics (1/Q
1
=0.071±0.026; 1/Q
3
=0.120±0.005;
1/Q
2
= 0.171±0.016).
The cross-correlations are one-sided for triplet JFS2-JFS1-JF00 (Fig. 8a), showing a strong
non-isotropic source distribution or asymmetric background attenuation. The inverted phase
velocity on the longest path JFS2-JF00 is about the average of those for the two shorter paths
JFS2-JFS1 and JFS1-JF00 (Fig. 8b). The inter-station distances are small, so the amplitude
measurements are done on relatively high frequencies (20-25 Hz) to approximate the attenuation
between stations. The amplitude curves (Fig. 8c) are measured by the frequency domain Hilbert
transform approach and are corrected based on Eq. (1). These curves are smoothed using a
moving average method with a length of 10 Hz. The amplitude curves before smoothing have
noticeable fluctuations on the semi-logarithm scale. This is probably because there are correlated
noise sources and/or strong variation due to random scattering and random excitation of noise,
and/or remaining small earthquakes (e.g. Ben-Zion et al. 2015) that cannot be completely
removed by the preprocessing. The Q values are inverted from the relative decay curves (Fig. 8d)
using Eq. (2) and a random resampling of 34 days of the cross-correlation results. Then the 34
days are randomly resampled 1000 times using bootstrap technique to create a set of Q
-1
values
from which the mean and standard deviation are calculated (Fig. 8e). The bootstrap results
indicate that the segment between JFS2 and JF00 has strong attenuation (Q
3
≈8,
1/Q
3
=0.120±0.005), which probably results from the fault damage zone. The three cross-
correlation functions show consistent asymmetry, suggesting that the Q values of two short
segments JFS2-JFS1 (Q
1
≈14, 1/Q
1
=0.071±0.026) and JFS1-JF00 (Q
2
≈6, 1/Q
2
=0.170±0.016) are
reliable as found in the synthetic test with non-isotropic source distribution.
122
Figure 9.Cross-correlation functions, phase velocity curves, amplitude measurements and
inverse Q distribution for JFS2-JF00-JFN2 (see Fig. 7 for station locations). (a) Cross-
correlations for JFS2-JF00-JFN2. (b) Phase velocity dispersion curves. (c) Corrected amplitude
decay curves on three correlation functions. The frequency range is 15-20 Hz related to the
slightly larger inter-station distance compare to the triplet JFS2-JFS1-JF00 in Fig. 8. (d) Relative
amplitude decay curves and inverted attenuation Q values from random resampling of 34
days.Best-fitting curves are plotted as solid lines with the same color as data. (e) Distribution of
123
inverse Q
R
values for three station pairs based on bootstrapping statistics (1/Q
1
=0.094±0.017;
1/Q
3
=0.041±0.006; 1/Q
2
=-0.005±0.011).
The second triplet JFS2-JF00-JFN2 has larger inter-station distance (Fig. 9a), the associated
analysis uses a lower frequency range (15-20 Hz), and the phase velocity dispersion curves are
more separated apart (Fig. 9b). The cross-correlation (Fig. 9a) is again asymmetric due to non-
isotropic source distribution and attenuation/scattering. As a result, the inverted Q values should
provide reliable estimates for the attenuation structure along the two shorter segments (JFS2-
JF00 and JF00-JFN2). However, the Signal-to-Noise Ratio is not as good as for the first triplet
JFS2-JFS1-JF00 due to the larger distance. The amplitude decay curve (Fig. 9c) is again
measured in frequency domain Hilbert transform. The southern pair JFS2-JF00 has strong
attenuation (Q
1
≈11, 1/Q
1
= 0.094±0.017), while the northern pair JF00-JFN2 has higher phase
velocity and its attenuation is not reliably determined (1/Q
2
= -0.005±0.011) with standard
deviation twice as large as the absolute value of mean in the bootstrap distribution. The
attenuation on the longest path JFS2-JFN2 is Q
3
≈24 (1/Q
3
=0.041±0.006). The pair JFS2-JF00
also corresponds to the longest path (Q
3
) in triplet JFS2-JFS1-JF00 (Fig. 8) and the inverted Q
values for this segment in the two triplets are in good agreement. The minor difference comes
from the shorter interstation distances and higher frequency band used in JFS2-JFS1-JF00 (Fig. 8)
leading to Q=8, compared to those used here for JFS2-JF00-JFN2 leading to Q=11. These results
suggest that the fault damage zone exists primarily to the SW of station JF00 (Figs. 7-9). This is
consistent with observations of fault zone trapped waves in that fault section (Qiu et al. 2015).
We also note that the attenuation corresponding to the longest path is more accurate (smaller
error of 1/Q
3
) than those of the two shorter segments.
124
Figure 10. Cross-correlation functions, amplitude measurements and inverse Q distribution for
triplet JFS3-JFS1-JFN1 (see locations in Fig. 7). (a) Cross-correlations for JFS3-JFS1-JFN1. (b)
Relative amplitude decay curve.Best-fitting curves are plotted as solid lines with the same color
as data.The frequency range is 15-25 Hz. (c) Distribution of inverse Q values for the 3 station
pairs in JFS3-JFS1-JFN1 based on bootstrapping statistics (1/Q
1
=0.198±0.037;
1/Q
3
=0.094±0.012; 1/Q
2
=0.049±0.018).
For the two additional triplets, JFS3-JFS1-JFN1 and JFN2-JFN1-JFS1, the cross-correlations
(panel a), relative decay curves (panel b) and bootstrap inverse Q distributions (panel c) are
shown in Fig. 10 (frequency band use for the inversion is 15-25 Hz) and Fig. 11 (frequency band
used for the inversion is 18-25 Hz), respectively. Because of the relatively large distance and
strong body phase associated with JFS3-JFN1, we apply time domain windowing to localize the
surface wave packet for amplitude measurements. The inverted Q values show that JFS3-JFS1
has stronger attenuation (Q
1
≈5, 1/Q
1
= 0.198±0.037) than the northern segment JFS1-JFN1
(Q
2
≈20, 1/Q
2
= 0.049±0.018), which is consistent with the previous two triplets. Interestingly, the
Q values for JFN2-JFN1-JFS1 show also that the northern segment JFN2-JFN1 has weaker
attenuation (Q
1
≈31, 1/Q
1
= 0.032±0.016) compared to the southern segment JFN1-JFS1 (Q
2
≈16,
1/Q
2
= 0.064±0.009). However, there are noticeable differences in the inverted Q values of JFN1-
125
JFS1 from these two triplets, which is probably due to relatively low SNR and large inter-station
distances.
Figure 11. Cross-correlation functions, amplitude measurements and inverse Q distribution for
JFN2-JFN1-JFS1. (a) Cross-correlations for JFN2-JFN1-JFS1. Interstation distances are
indicated in titles.(b) Relative amplitude decay curve.Best-fitting curves are plotted as solid lines
with the same color as data.The frequency range is 18-25 Hz. (c) Distribution of inverse Q values
for the 3 station pairs in JFN2-JFN1-JFS1 based on bootstrapping statistics (1/Q
1
=0.032±0.016;
1/Q
3
=0.049±0.003; 1/Q
2
=0.064±0.009).
4.6. Discussion
Deriving attenuation from ambient noise cross-correlation can provide important information
on structural and wave propagation properties between recording stations. Based on theoretical
results on amplitude decay of ensemble-averaged cross-correlations due to attenuation (Liu and
Ben-Zion 2013), we develop an inversion algorithm to retrieve inter-station Q
R
values from
cross-spectra of ambient noise data by using linear triplets of stations. For this particular
geometry, the three cross-correlation functions share common noise sources and background
attenuation in the causal (or anti-causal) direction. Using redundant information in the triplet and
assuming that the coherent noise comes mainly from the end-fire lobe direction (with interstation
126
distance > wavelength), the background attenuation, noise source spectra and site amplification
can be removed from the analysis. This allows us to form a linear least-square inversion (Eq. 2)
for inter-station Q values using cross-spectra in the same frequency band at the three station pairs.
The site amplification factors of the three stations are assumed constant or following similar
pattern in the used frequency range.
The proposed method complements previous studies on deriving attenuation information
between pairs of stations. In particular, Weaver (2013) developed a formulation for obtaining
attenuation based on radiative transfer equation for linear array of at least four stations. The
method is applied to narrow-band filtered cross-correlation functions in time domain and
involves measuring the peak amplitude on both causal or anticausal part across the array to invert
for attenuation coefficient and site amplification. The formulation of Weaver (2013) was proven
theoretically and applied to real data, but it requires more stations and data (at least six good
cross-correlation functions with both causal and anti-causal parts) than the inversion proposed
here. The basic equations in this paper and Weaver (2013) are consistent although they are
obtained by different approaches; our results on amplitude and phase terms are derived in the
frequency domain while Weaver (2013) provides the amplitude of narrow-band filtered noise
cross-correlation in time. The causal part of Eq. (1) agrees with Eq. (1) of Weaver (2013) in that
they both have exponential decay term, amplitude correction term (geometrical decay term in
Weaver (2013)) and noise source intensity term.
Our inversion formula is validated with synthetic noise cross-correlations and is then applied
to estimate Q values along linear arrays of two different scales in southern California. The
synthetic test results suggest that the proposed method can better resolve Q values of the two
inner station pairs for non-isotropic source distribution case with maximum intensity parallel
with inter-station direction than isotropic source distribution. This is likely because the noise
waveforms propagating in the latter case outside the end-fire lobes do not cancel out each other
sufficiently in heterogeneous dissipative medium. In application of the method to a station triplet
within the regional southern California network with inter-station distances of ~35 km (Fig. 5),
the frequency band for attenuation inversion is relatively low (0.20 to 0.36 Hz; slightly higher
than the ocean microseism peak). The obtained Q
R
values are ~25 along the used triplet and
likely reflect attenuation in the top 5 km of the crust. The depth is based on 1/3 of wavelength
that the Rayleigh wave is most sensitive to (Nolet 1987). The Q
R
values are significantly lower
127
than the body wave Q
S
~ 400 at 5 km depth estimated by Hauksson and Shearer (2006) based on
S wave spectra of earthquakes associated with a much higher frequency band.
Applying the inversion algorithm to 6 stations (4 triplets) of a linear array with smaller
interstation spacing (~25 m) that crosses the San Jacinto fault zone (Fig. 7), we use higher
frequencies in the range 10-25 Hz. The analysis involves four triplets (from JFS3 to JFN2),
excluding several long-distance triplets with weak surface wave signal and low SNR. The results
point generally to a significant fault damage zone at the southern part of JF array, between
stations JFS3 and JF00, with Q
R
values between 6 and 11. This is consistent with observations of
fault zone trapped waves at these stations (Qiu et al. 2015). The northern two stations, JFN2 and
JFN1, have higher Q
R
values of ~20 to 30. The results likely reflect attenuation in the top 30 m
of the crust (again roughly estimated based on 1/3 of the dominant used wavelength). The
derived Q
R
values are in the same range of Q
P
and Q
S
results obtained from earthquake data
recorded in two shallow boreholes in the area (Aster and Shearer 1991). The results between
stations JFS3 and JF00 are slightly lower than Q
S
values obtained by waveform inversions of
trapped waves at the JF and nearby site (Qiu et al. 2015; Lewis et al 2005), probably because the
trapped waves results reflect average values in the top ~3.5 km of the crust.
A bootstrap technique enables us to quantify the variability of the inverted Q
-1
values over
different combinations of days. The ambient noise wavefield has seasonal variations (e.g. Stehly
et al. 2006; Hiller and Ben-Zion 2011; Zhan et al. 2013). This produces variations in amplitude
spectra of ambient noise cross-correlation over 2-3 months for the ocean microseism band
(~0.05-0.2 Hz). There can be additional short-term variations for the amplitude of high frequency
noise cross-correlation (~15-25 Hz) included in our study. Assuming that the attenuation/velocity
structure do not vary significantly within the time range of our data, the inverted attenuation Q
-1
values should be stable for different resampling sets of days. The bootstrap results for both the
high frequency linear array and the low frequency regional triplet suggest that most of the
inverted Q
-1
values have relatively narrow confidence intervals.
The requirement of constant or similar site amplification factors in the frequency range used
for inversion for the three stations limits the choice of triplets because the three stations should
not be located in very different geological units (e.g. one in thick sedimentary basin and the
others in hard rock or on a mountain). Otherwise the effect of site amplification may not be
128
canceled properly in Eq. (2). The linear triplet geometry is a limitation on the number of
available paths from which we can get attenuation (compared to all paths among randomly
located 3 stations). Nevertheless, as shown in the paper with examples from regional and fault
zone arrays, the method can be applied to many stations over wide range of scales. In particular,
the method is suitable for regularly spaced arrays on local or regional (e.g. USArray) scales
where many triplets can be formed to cover the region of interest.
The dense linear array requires high frequency noise that samples the upper tens of meters of
crust. The advantage of high frequency noise is that it requires less number of days for the cross-
correlation to converge. One disadvantage is that body wave phases can have comparable
amplitudes to fundamental Rayleigh waves at larger distances due to attenuation of surface wave
noise, which could potentially bias amplitude measurements. For example, a body P wave
appears on some of the cross-correlations with high apparent velocities and can affect the phase
and amplitude measurements of the surface wave. We derive both ZZ and RR cross-correlation
components for the high frequency noise and then select the RR component for amplitude
measurements to minimize the effect of body P wave phases.
The essential signal preprocessing steps should be improved in future. It is desirable to have
a processing scheme that preserves most information in the noise cross-correlation. Pre-
whitening is not included in our analysis, since it modifies the amplitude decay curve (Weemstra
et al. 2014). Instead, a time domain normalization is performed by dividing each trace segment
of the ambient noise by its total power before cross-correlation. This strengthens the stationary
noise assumption which could enhance the convergence of cross-correlation. Earthquakes and
other unwanted signals (e.g. human activity) should be carefully removed (here we use
amplitude clipping for small events and remove segments with larger earthquakes). The clipping
method cannot completely remove small earthquakes that can pollute the noise data. Numerous
tiny events with small SNR cannot be detected and may cause bias on the amplitude estimates.
Removing many earthquake segments will reduce the length of available data and make the
cross-correlations slower to converge. Different steps of preprocessing procedures should be
experimented to determine improved processing flow in the future.
In the inversion step, the quality of results suffer from artifacts due to filtering and
windowing in the process of isolating the fundamental mode surface wave packet in the causal or
129
anti-causal parts. The filter and windowing function should be optimized to minimize these
artifacts. Generally, the goal is to select only the fundamental mode surface wave packet based
on frequency and velocity constraints without introducing distortion of amplitude spectra due to
the transition band of filter. It is challenging to select the desired signals from cross-correlations
due to potential overlaps with other surface wave modes, body waves and other unwanted energy.
Improvements in various signal processing steps will be attempted in a future work.
Appendix 4A: Derivation of cross-spectrum for isotropic noise sources and asymmetric
attenuation structure
We assume the two stations are symmetric about the origin of a coordinate system, which is
the center of a ring containing the isotropic noise sources in the far field. The background
attenuation quality factor values in the West (Q
outW
) and East (Q
outE
) half planes, are slightly
different (Fig. 1a). Based on Eq. (5) in Liu and Ben-Zion (2013), we have
d
c c
i R RdR
R
R
u u
~
) (
~
) (
ω exp ) ω ; , ( B
) ( ) (
1
ω ;
2
x
,
2
x
C E
2
*
1
2 1
2
0
max
min
2 1
(A1)
where
ω ;
2
x
,
2
x
C E
2 1
u u
is the expected cross-spectrum function, ) ω ; , ( B R is the far-field
noise source spectral density at polar coordinate location ) , ( R with central frequency ω, and
) (
1
and ) (
2
are distances from source at angle θ to receivers 1 and 2, respectively. R is the
distance between the far field noise source and the origin. The inter-station distance is x. The
inverse complex phase velocity is defined the in a similar way as in Liu and Ben-Zion (2013):
c Q i c ) ( 2 / sgn 1
~
1 , where Q is the attenuation quality factor.
We assign different attenuation Q values to different regions: Q
in
for inter-station region,
Q
outW
for background attenuation in the west half plane and Q
outE
for background attenuation in
the east half plane (see Fig. 1a). With stationary phase approximation (Fig. 1b), only the noise
propagating within the end-fire lobes (Roux et al. 2004, Snieder et al. 2004) contribute
constructively to the cross-correlation of which the causal part decay is determined by Q
in
and
Q
outW
and the anti-causal part decay is determined by Q
in
and Q
outE
. We also assume the Fresnel
approximation that source-receiver distances are much greater than the inter-station distance. As
a result, we have 2 cos ) (
1
x R and 2 cos ) (
2
x R , where R is the average
130
distance from noise sources to origin. By evaluating the integrals in Eq. (A1) for two half planes,
we arrive at
c
x
i
c
x
J
Q c
R
c
x
i
c
x
J
Q c
R
E E
outE
W W
outW
u u
0 0
0 0
H ω exp ) ω ( B
H ω exp ) ω ( B ω ;
2
x
,
2
x
C E
2 1
(A2)
where
outW in
W
Q Q
i
1 1
2
sgn
1
and
outE in
E
Q Q
i
1 1
2
sgn
1
are complex
coefficients containing the attenuation Q factors.
0
J represents zero.th order Bessel function of
1
st
kind and
0
H represents zero.th order Struve function of 1
st
kind. This equation, when narrow-
band filtered, shows two wave packets propagating from west and east directions, respectively.
We simplify Eq. (A2) assuming that the interstation distance is greater than the wavelength:
x 2 . As a result, we have 1 c x c x
E
. The Bessel function terms can be
simplified
c
x
H
c
x
iY
c
x
J
c
x
i
c
x
J
E E E E E ) 1 (
0 0 0 0 0
H , (A3))
where
) 1 (
0
H represents zero.th order Hankel function of the 1
st
kind and
0
Y represents zero.th
order Bessel function of 2
nd
kind. The Hankel function can be further simplified (e.g.,
Mathematics Handbook; Wikipedia; Boschi et al., 2013)
c
x
i
c
x
x
c i
c
x
i
x
c i
c
x
H
E
E
E
E
E
exp
Im
ω exp
2
exp
2
) 1 (
0
, (A4)
where the three factors on the right-hand-side are amplitude correction factor, attenuation
exponential decay factor and wave propagation factor, respectively. Similarly, the Bessel
function term representing wave traveling from west can be simplified as
(A5)
131
Substituting Eq. (A3)-(A5) back to Eq. (A2), we have
(A6)
We arrive at Eq. (1) by multiplying site amplification factors β
1
β
2
on Eq. (A6).
132
Discussion
I develop tools for analyzing statistically the ambient seismic noise data and extracting
attenuation Q values from amplitude information of ambient seismic noise cross-correlations.
The statistical tools include methods for measuring the non-diffuse characteristics of ambient
seismic noise wave field and estimating the error of stacked ambient noise cross-spectrum.
Together these tools provide important information on the origin and propagation/scattering of
ambient seismic noise, the quality of non-diffuse noise wave field for imaging purposes and the
error of phase and amplitude information that can be applied to constructing accurate travel time
and attenuation/site amplification tomographic maps.
The method on correlations of neighboring power spectral values at different frequencies
has been applied to broadband stations in southern California. The results show similar
correlated frequency patterns for mountain (CHF) and desert (LMR2) stations, but a different
pattern for a Los Angeles basin station (OLI) with more non-diffuse wave components at
frequencies above 0.25 Hz. This method will be applied in a future study to many broadband
stations in southern California to better analyze the spatio-temporal distributions of non-diffuse
ambient seismic wave field. Furthermore, the non-diffuse wave field will be decomposed into
cross-frequency noise components and fully diffuse noise as introduced in the simulation test of
chapter 1.4. There is also a potential to apply this method to long-period (>10 s) coda wave of
large earthquakes with magnitude greater than 6.5. This method requires sufficiently long noise
data recordings with enough number of evenly spaced time windows to converge. Insufficient
number of time windows will reduce the stability of the correlation coefficient matrix of power
spectral values.
The error analysis of stacked noise cross-spectrum has been applied to stations in both the
regional southern California seismic network across the San Andreas fault and a dense linear
array across the San Jacinto fault zone. The continuous noise recordings contain outliers (e.g.
earthquakes and other transient signals) that can be removed statistically based on distribution of
cross-spectrum at each frequency. For the regional network stations, the stacked cross-spectra
curve from outlier windows has 3 times greater amplitude with more variations than that from
selected windows, which can be problematic for extracting site factor and attenuation from
amplitude information if the outliers are not properly removed. For the dense array across the
San Jacinto fault, the outliers can significantly reduce the SNR of stacked cross-spectra curve by
a factor of ~4, probably due to more impulsive signals at high frequencies (10-50 Hz). After
133
removing outliers, the phase and amplitude errors can be used to derive errors in phase velocity
and attenuation/site factor, which can be applied to mapping temporal changes in subsurface
structure. The different noise windows are assumed independent and having similar statistics in
order to estimate errors. However, the autocorrelation test suggests temporal correlations for
frequency below 0.2 Hz. The error bars for stacked cross-spectra within that frequency band can
be estimated through a block bootstrap method.
The method of estimating attenuation from linear array of three stations requires accurate
amplitude measurements and isolation of fundamental mode surface wave. Two commonly used
noise preprocessing techniques, onebit and prewhitening, are non-linear operations and can bring
in strong bias in estimating attenuation and site amplifications. Therefore it is better to follow the
approach in chapter 2 to remove outliers statistically. Combined with the error analysis of
stacked cross-spectral amplitude, the confidence intervals for Q
-1
values can be estimated
through error propagation and this can be verified by bootstrap resampling. The stacked cross-
correlation after outlier removal may still contain higher mode surface waves and body waves,
which can produce bias if strong enough compared with the fundamental mode surface wave.
Therefore it is important to isolate the wave packet of fundamental mode surface wave from all
the higher modes and body waves. This method doesn’t account for inter-station wave
focusing/defocusing due to heterogeneous velocity structure, the effect of which can be
estimated using numerical simulations. The method can separate attenuation from site
amplification assuming similar site conditions for all three stations.
134
References
Aki, K., 1957, Space and time spectra of stationary stochastic waves, with spectral reference to
microtremors: Bull. Earthq. Res. Inst., v. 35, p. 415–456.
Allam, A.A., and Ben-Zion, Y., 2012, Seismic velocity structures in the southern California
plate-boundary environment from double-difference tomography: Tomography of the San
Jacinto fault area: Geophysical Journal International, v. 190, no. 2, p. 1181–1196, doi:
10.1111/j.1365-246X.2012.05544.x.
Aster, R.C., and Shearer, P.M., 1991, High-frequency borehole seismograms recorded in the San
Jacinto Fault zone, Southern California Part 2. Attenuation and site effects: Bulletin of the
Seismological Society of America, v. 81, no. 4, p. 1081–1100.
Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro,
N.M., and Yang, Y., 2007, Processing seismic ambient noise data to obtain reliable broad-
band surface wave dispersion measurements: Geophysical Journal International, v. 169, no.
3, p. 1239–1260, doi: 10.1111/j.1365-246X.2007.03374.x.
Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J.G., Ozer, N., Michael, A.J., Baris,
S., and Aktar, M., 2003, A shallow fault-zone structure illuminated by trapped waves in the
Karadere-Duzce branch of the North Anatolian Fault, western Turkey: Geophysical Journal
International, v. 152, no. 3, p. 699–717, doi: 10.1046/j.1365-246X.2003.01870.x.
Ben-Zion, Y., Vernon, F.L., Ozakin, Y., Zigone, D., Ross, Z.E., Meng, H., White, M., Reyes, J.,
Hollis, D., and Barklage, M., 2015, Basic data features and results from a spatially dense
seismic array on the San Jacinto fault zone: Geophysical Journal International, v. 202, no.
1, p. 370–380, doi: 10.1093/gji/ggv142.
Boashash, B., 1992, Estimating and interpreting the instantaneous frequency of a signal. I.
Fundamentals: Proceedings of the IEEE, v. 80, no. 4, p. 520–538, doi: 10.1109/5.135376.
Boschi, L., Weemstra, C., Verbeke, J., Ekström, G., Zunino, A., and Giardini, D., 2013, On
measuring surface wave phase velocity from station–station cross-correlation of ambient
signal: Geophysical Journal International, v. 192, no. 1, p. 346–358, doi:
10.1093/gji/ggs023.
Boué, P., Roux, P., Campillo, M., and Briand, X., 2014, Phase velocity tomography of surface
waves using ambient noise cross correlation and array processing: Journal of Geophysical
Research: Solid Earth, v. 119, no. 1, p. 519–529, doi: 10.1002/2013JB010446.
135
Brenguier, F., Campillo, M., Hadziioannou, C., Shapiro, N.M., Nadeau, R.M., and Larose, E.,
2008, Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous
Seismological Observations: Science, v. 321, no. 5895, p. 1478–1481, doi:
10.1126/science.1160943.
Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N.M., Briand, X., Emoto, K., and
Miyake, H., 2014, Mapping pressurized volcanic fluids from induced crustal seismic
velocity drops: Science, v. 345, no. 6192, p. 80–82, doi: 10.1126/science.1254073.
Brenguier, F., Clarke, D., Aoki, Y., Shapiro, N.M., Campillo, M., and Ferrazzini, V., 2011,
Monitoring volcanoes using seismic noise correlations: ComptesRendus Geoscience, v.
343, no. 8-9, p. 633–638, doi: 10.1016/j.crte.2010.12.010.
Campillo, M., 2006, Phase and Correlation in `Random’ Seismic Fields and the Reconstruction
of the Green Function: Pure and Applied Geophysics, v. 163, no. 2-3, p. 475–502, doi:
10.1007/s00024-005-0032-8.
Campillo, M., Sato, H., Shapiro, N.M., and van der Hilst, R.D., 2011, New developments on
imaging and monitoring with seismic noise: ComptesRendus Geoscience, v. 343, no. 8-9,
p. 487–495, doi: 10.1016/j.crte.2011.07.007.
Claerbout, J.F., 1968, Synthesis of a layered medium from its acoustic transmission response:
Geophysics, v. 33, no. 2, p. 264–269, doi: 10.1190/1.1439927.
Cupillard, P., and Capdeville, Y., 2010, On the amplitude of surface waves obtained by noise
correlation and the capability to recover the attenuation: a numerical approach:
Geophysical Journal International, doi: 10.1111/j.1365-246X.2010.04586.x.
Dziewonski, A., Bloch, S., and Landisman, M., 1969, A technique for the analysis of transient
seismic signals: Bulletin of the seismological Society of America, v. 59, no. 1, p. 427–444.
Galetti, E., and Curtis, A., 2012, Generalised receiver functions and seismic interferometry:
Tectonophysics, v. 532–535, p. 1–26, doi: 10.1016/j.tecto.2011.12.004.
Gouédard, P., Stehly, L., Brenguier, F., Campillo, M., Colin de Verdière, Y., Larose, E.,
Margerin, L., Roux, P., Sánchez-Sesma, F.J., Shapiro, N.M., and Weaver, R.L., 2008,
Cross-correlation of random fields: mathematical approach and applications: Geophysical
Prospecting, v. 56, no. 3, p. 375–393, doi: 10.1111/j.1365-2478.2007.00684.x.
Hauksson, E., and Shearer, P.M., 2006, Attenuation models (QP and QS) in three dimensions of
the southern California crust: Inferred fluid saturation at seismogenic depths: Journal of
Geophysical Research: Solid Earth, v. 111, no. B5, p. n/a–n/a, doi: 10.1029/2005JB003947.
136
Hennino, R., Trégourès, N., Shapiro, N.M., Margerin, L., Campillo, M., van Tiggelen, B.A., and
Weaver, R.L., 2001, Observation of Equipartition of Seismic Waves: Physical Review
Letters, v. 86, no. 15, p. 3447–3450, doi: 10.1103/PhysRevLett.86.3447.
Herrmann, R.B., 1973, Some aspects of band-pass filtering of surface waves: Bulletin of the
Seismological Society of America, v. 63, no. 2, p. 663–671.
Hillers, G., Campillo, M., Ben-Zion, Y., and Roux, P., 2014, Seismic fault zone trapped noise:
Journal of Geophysical Research: Solid Earth, v. 119, no. 7, p. 5786–5799, doi:
10.1002/2014JB011217.
Hillers, G., and Ben-Zion, Y., 2011, Seasonal variations of observed noise amplitudes at 2-18 Hz
in southern California: Seasonal variations of high-f seismic noise: Geophysical Journal
International, v. 184, no. 2, p. 860–868, doi: 10.1111/j.1365-246X.2010.04886.x.
Hillers, G., Ben-Zion, Y., Campillo, M., and Zigone, D., 2015, Seasonal variations of seismic
velocities in the San Jacinto fault area observed with ambient seismic noise: Geophysical
Journal International, v. 202, no. 2, p. 920–932, doi: 10.1093/gji/ggv151.
Hillers, G., Ben-Zion, Y., Landès, M., and Campillo, M., 2013, Interaction of microseisms with
crustal heterogeneity: A case study from the San Jacinto fault zone area: Geochemistry,
Geophysics, Geosystems, v. 14, no. 7, p. 2182–2197, doi: 10.1002/ggge.20140.
Kawase, H., Matsushima, S., Satoh, T., and Sánchez-Sesma, F.J., 2015, Applicability of
Theoretical Horizontal-to-Vertical Ratio of Microtremors Based on the Diffuse Field
Concept to Previously Observed Data: Bulletin of the Seismological Society of America,
doi: 10.1785/0120150134.
Kimman, W.P., and Trampert, J., 2010, Approximations in seismic interferometry and their
effects on surface waves: Approximations in seismic interferometry: Geophysical Journal
International, p. no–no, doi: 10.1111/j.1365-246X.2010.04632.x.
Landès, M., Hubans, F., Shapiro, N.M., Paul, A., and Campillo, M., 2010, Origin of deep ocean
microseisms by using teleseismic body waves: Journal of Geophysical Research, v. 115,
no. B5, p. B05302, doi: 10.1029/2009JB006918.
Lawrence, J.F., and Prieto, G.A., 2011, Attenuation tomography of the western United States
from ambient seismic noise: Journal of Geophysical Research, v. 116, no. B6, p. B06302,
doi: 10.1029/2010JB007836.
137
Lewis, M.A., Peng, Z., Ben-Zion, Y., and Vernon, F.L., 2005, Shallow seismic trapping structure
in the San Jacinto fault zone near Anza, California: Geophysical Journal International, v.
162, no. 3, p. 867–881, doi: 10.1111/j.1365-246X.2005.02684.x.
Lin, F.-C., Moschetti, M.P., and Ritzwoller, M.H., 2008, Surface wave tomography of the
western United States from ambient seismic noise: Rayleigh and Love wave phase velocity
maps: Geophysical Journal International, v. 173, no. 1, p. 281–298, doi: 10.1111/j.1365-
246X.2008.03720.x.
Lin, F.-C., and Ritzwoller, M.H., 2011, Helmholtz surface wave tomography for isotropic and
azimuthally anisotropic structure: Helmholtz surface wave tomography: Geophysical
Journal International, v. 186, no. 3, p. 1104–1120, doi: 10.1111/j.1365-246X.2011.05070.x.
Lin, F.-C., and Tsai, V.C., 2013, Seismic interferometry with antipodal station pairs:
ANTIPODAL SEISMIC INTERFEROMETRY: Geophysical Research Letters, v. 40, no.
17, p. 4609–4613, doi: 10.1002/grl.50907.
Lin, F.-C., Tsai, V.C., and Ritzwoller, M.H., 2012, The local amplification of surface waves: A
new observable to constrain elastic velocities, density, and anelastic attenuation: Journal of
Geophysical Research, v. 117, no. B6, p. B06302, doi: 10.1029/2012JB009208.
Liu, X., and Ben-Zion, Y., 2013, Theoretical and numerical results on effects of attenuation on
correlation functions of ambient seismic noise: Geophysical Journal International, v. 194,
no. 3, p. 1966–1983, doi: 10.1093/gji/ggt215.
Liu, X., Ben-Zion, Y., and Zigone, D., 2015, Extracting seismic attenuation coefficients from
cross-correlations of ambient noise at linear triplets of stations: Geophysical Journal
International, v. 203, no. 2, p. 1149–1163, doi: 10.1093/gji/ggv357.
Liu, X., and Ben-Zion, Y., 2016, Estimating correlations of neighboring frequencies in ambient
seismic noise: Geophysical Journal International, p. in review.
Liu, X., Ben-Zion, Y., and Zigone, D., 2016, Frequency Domain Analysis of Errors in Cross-
Correlations of Ambient Seismic Noise: Geophysical Journal International, p. in review.
Lobkis, O.I., and Weaver, R.L., 2001, On the emergence of the Green’s function in the
correlations of a diffuse field: The Journal of the Acoustical Society of America, v. 110, no.
6, p. 3011, doi: 10.1121/1.1417528.
Margerin, L., Campillo, M., Van Tiggelen, B.A., and Hennino, R., 2009, Energy partition of
seismic coda waves in layered media: theory and application to Pinyon Flats Observatory:
138
Geophysical Journal International, v. 177, no. 2, p. 571–585, doi: 10.1111/j.1365-
246X.2008.04068.x.
Nolet, G., 1987, Seismic tomography: with applications in global seismology and exploration
geophysics: Springer Science & Business Media.
Pedersen, H.A., Mars, J.I., and Amblard, P.-O., 2003, Improving surface-wave group velocity
measurements by energy reassignment: Geophysics, v. 68, no. 2, p. 677–684.
Peng, Z., Ben-Zion, Y., Michael, A.J., and Zhu, L., 2003, Quantitative analysis of seismic fault
zone waves in the rupture zone of the 1992 Landers, California, earthquake: evidence for a
shallow trapping structure: Geophysical Journal International, v. 155, no. 3, p. 1021–1041,
doi: 10.1111/j.1365-246X.2003.02109.x.
Peterson, J.R., 1993, Observations and Modeling of Seismic Background Noise (No. 93-322):
Geological Survey (US).
Poli, P., Campillo, M., Pedersen, H., and LAPNET Working Group, 2012, Body-Wave Imaging
of Earth’s Mantle Discontinuities from Ambient Seismic Noise: Science, v. 338, no. 6110,
p. 1063–1065, doi: 10.1126/science.1228194.
Poli, P., Pedersen, H.A., Campillo;, M., and the POLENET/LAPNET Working Group, 2013,
Noise directivity and group velocity tomography in a region with small velocity contrasts:
the northern Baltic shield: Geophysical Journal International, v. 192, no. 1, p. 413–424,
doi: 10.1093/gji/ggs034.
Prieto, G.A., Denolle, M., Lawrence, J.F., and Beroza, G.C., 2011, On amplitude information
carried by the ambient seismic field: ComptesRendus Geoscience, v. 343, no. 8-9, p. 600–
614, doi: 10.1016/j.crte.2011.03.006.
Prieto, G.A., Lawrence, J.F., and Beroza, G.C., 2009, Anelastic Earth structure from the
coherency of the ambient seismic field: Journal of Geophysical Research, v. 114, no. B7, p.
B07303, doi: 10.1029/2008JB006067.
Qiu, H., Ben-Zion, Y., Ross, Z.E., Share, P., and Vernon, F.L., 2015, Internal structure of the
San Jacinto fault zone at Jackass Flat from earthquake data recorded by a dense linear
array: Abstract of the annual SSA meeting,.
Rickett, J.E., and Claerbout, J.F., 2000, Calculation of the sun’s acoustic impulse response by
multi-dimensional spectral factorization: Solar Physics, v. 192, no. 1, p. 203–210, doi:
10.1023/A:1005205406377.
139
Roux, P., Kuperman, W.A., Group, N., and others, 2004, Extracting coherent wave fronts from
acoustic ambient noise in the ocean: The Journal of the Acoustical Society of America, v.
116, no. 4, p. 1995–2003.
Roux, P., Sabra, K.G., Gerstoft, P., Kuperman, W.A., and Fehler, M.C., 2005, P-waves from
cross-correlation of seismic noise: P WAVES FROM NOISE CROSS-CORRELATIONS:
Geophysical Research Letters, v. 32, no. 19, p. n/a–n/a, doi: 10.1029/2005GL023803.
Sabra, K.G., 2005, Extracting time-domain Green’s function estimates from ambient seismic
noise: Geophysical Research Letters, v. 32, no. 3, doi: 10.1029/2004GL021862.
Sabra, K.G., Gerstoft, P., Fehler, M.C., Gerstoft, P., Roux, P., and Kuperman, W.A., 2005a,
Extracting time-domain Green’s function estimates from ambient seismic noise:
Geophysical Research Letters, v. 32, no. 3, doi: 10.1029/2004GL021862.
Sabra, K.G., Roux, P., and Kuperman, W.A., 2005b, Emergence rate of the time-domain Green’s
function from the ambient noise cross-correlation function: The Journal of the Acoustical
Society of America, v. 118, no. 6, p. 3524, doi: 10.1121/1.2109059.
Sanchez-Sesma, F.J., and Campillo, M., 2006, Retrieval of the Green’s Function from Cross
Correlation: The Canonical Elastic Problem: Bulletin of the Seismological Society of
America, v. 96, no. 3, p. 1182–1191, doi: 10.1785/0120050181.
Sánchez-Sesma, F.J., Pérez-Ruiz, J.A., Luzón, F., Campillo, M., and Rodríguez-Castellanos, A.,
2008, Diffuse fields in dynamic elasticity: Wave Motion, v. 45, no. 5, p. 641–654, doi:
10.1016/j.wavemoti.2007.07.005.
Schechinger, B., Goertz, A., Artman, B., Lambert, M.-A., Koerbe, M., Krajewski, P., and others,
2009, Extracting subsurface information from ambient seismic noise-an example from
Germany, in 2009 SEG Annual Meeting, Society of Exploration Geophysicists.
Schulte-Pelkum, V., Earle, P.S., and Vernon, F.L., 2004, Strong directivity of ocean-generated
seismic noise: OCEAN MICROSEISMIC NOISE: Geochemistry, Geophysics,
Geosystems, v. 5, no. 3, p. Q03004, doi: 10.1029/2003GC000520.
Sens-Schönfelder, C., Snieder, R., and Stähler, S.C., 2015, The lack of equipartitioning in global
body wave coda: Geophysical Research Letters, v. 42, no. 18, p. 7483–7489, doi:
10.1002/2015GL065108.
Sens-Schönfelder, C., and Wegler, U., 2006, Passive image interferometry and seasonal
variations of seismic velocities at Merapi Volcano, Indonesia: Geophysical Research
Letters, v. 33, no. 21, doi: 10.1029/2006GL027797.
140
Shapiro, N.M., and Campillo, M., 2004, Emergence of broadband Rayleigh waves from
correlations of the ambient seismic noise: Geophysical Research Letters, v. 31, no. 7, p.
L07614, doi: 10.1029/2004GL019491.
Shapiro, N.M., Campillo, M., Stehly, L., and Ritzwoller, M.H., 2005, High-resolution surface-
wave tomography from ambient seismic noise: Science, v. 307, no. 5715, p. 1615–1618,
doi: 10.1126/science.1108339.
Snieder, R., 2004, Extracting the Green’s function from the correlation of coda waves: A
derivation based on stationary phase: Physical Review E, v. 69, no. 4, p. 046610, doi:
10.1103/PhysRevE.69.046610.
Snieder, R., Sánchez-Sesma, F.J., and Wapenaar, K., 2009, Field Fluctuations, Imaging with
Backscattered Waves, a Generalized Energy Theorem, and the Optical Theorem: SIAM
Journal on Imaging Sciences, v. 2, no. 2, p. 763–776, doi: 10.1137/08072913X.
Stehly, L., Campillo, M., and Shapiro, N.M., 2006, A study of the seismic noise from its long-
range correlation properties: Journal of Geophysical Research, v. 111, no. B10, p. B10306,
doi: 10.1029/2005JB004237.
Van Tiggelen, B.A., 2003, Green Function Retrieval and Time Reversal in a Disordered World:
Physical Review Letters, v. 91, no. 24, doi: 10.1103/PhysRevLett.91.243904.
Traer, J., Gerstoft, P., Bromirski, P.D., and Shearer, P.M., 2012, Microseisms and hum from
ocean surface gravity waves: Journal of Geophysical Research, v. 117, no. B11, doi:
10.1029/2012JB009550.
Tsai, V.C., 2011, Understanding the amplitudes of noise correlation measurements: Journal of
Geophysical Research: Solid Earth, v. 116, no. B9, p. B09311, doi:
10.1029/2011JB008483.
Wapenaar, K., 2004, Retrieving the Elastodynamic Green’s Function of an Arbitrary
Inhomogeneous Medium by Cross Correlation: Physical Review Letters, v. 93, no. 25, doi:
10.1103/PhysRevLett.93.254301.
Weaver, R.L., 1982, On diffuse waves in solid media: The Journal of the Acoustical Society of
America, v. 71, no. 6, p. 1608, doi: 10.1121/1.387816.
Weaver, R.L., 2011, On the amplitudes of correlations and the inference of attenuations, specific
intensities and site factors from ambient noise: ComptesRendus Geoscience, v. 343, no. 8-
9, p. 615–622, doi: 10.1016/j.crte.2011.07.001.
141
Weaver, R.L., 2013, On the retrieval of attenuation and site amplifications from ambient noise
on linear arrays: further numerical simulations: Geophysical Journal International, v. 193,
no. 3, p. 1644–1657, doi: 10.1093/gji/ggt063.
Weaver, R., Froment, B., and Campillo, M., 2009, On the correlation of non-isotropically
distributed ballistic scalar diffuse waves: The Journal of the Acoustical Society of America,
v. 126, no. 4, p. 1817, doi: 10.1121/1.3203359.
Weaver, R.L., and Lobkis, O.I., 2004, Diffuse fields in open systems and the emergence of the
Green’s function (L): The Journal of the Acoustical Society of America, v. 116, no. 5, p.
2731, doi: 10.1121/1.1810232.
Weaver, R.L., and Lobkis, O.I., 2005, Fluctuations in diffuse field–field correlations and the
emergence of the Green’s function in open systems: The Journal of the Acoustical Society
of America, v. 117, no. 6, p. 3432, doi: 10.1121/1.1898683.
Webb, S.C., 2008, The Earth’s hum: the excitation of Earth normal modes by ocean waves:
Geophysical Journal International, v. 174, no. 2, p. 542–566, doi: 10.1111/j.1365-
246X.2008.03801.x.
Weemstra, C., Westra, W., Snieder, R., and Boschi, L., 2014, On estimating attenuation from the
amplitude of the spectrally whitened ambient seismic field: Geophysical Journal
International, v. 197, no. 3, p. 1770–1788, doi: 10.1093/gji/ggu088.
Yao, H., and Van Der Hilst, R.D., 2009, Analysis of ambient noise energy distribution and phase
velocity bias in ambient noise tomography, with application to SE Tibet: Geophysical
Journal International, v. 179, no. 2, p. 1113–1132, doi: 10.1111/j.1365-246X.2009.04329.x.
Yao, H., van Der Hilst, R.D., and De Hoop, M.V., 2006, Surface-wave array tomography in SE
Tibet from ambient seismic noise and two-station analysis–I. Phase velocity maps:
Geophysical Journal International, v. 166, no. 2, p. 732–744, doi: 10.1111/j.1365-
246X.2006.03028.x.
Zhan, Z., Tsai, V.C., and Clayton, R.W., 2013, Spurious velocity changes caused by temporal
variations in ambient noise frequency content: Geophysical Journal International, v. 194,
no. 3, p. 1574–1581, doi: 10.1093/gji/ggt170.
Zhan, Z., Tsai, V.C., Jackson, J.M., and Helmberger, D., 2014, Ambient noise correlation on the
Amery Ice Shelf, East Antarctica: Geophysical Journal International, v. 196, no. 3, p.
1796–1802, doi: 10.1093/gji/ggt488.
142
Zigone, D., Ben-Zion, Y., Campillo, M., and Roux, P., 2015, Seismic Tomography of the
Southern California Plate Boundary Region from Noise-Based Rayleigh and Love Waves:
Pure and Applied Geophysics, v. 172, no. 5, p. 1007–1032, doi: 10.1007/s00024-014-0872-
1.
Abstract (if available)
Abstract
I study statistical properties of ambient seismic noise wave field and possible errors and biases in the empirical Green’s function obtained from ambient seismic noise cross-correlations. The techniques I developed can be applied to understanding the spatio-temporal non-diffuse characteristics of ambient seismic noise field and estimating errors of phase and amplitude information of empirical Green’s functions (noise cross-correlations) computed from finite amount of noise data. In addition, I derive the theoretical form of noise cross-correlation functions for different inter-station/background attenuation Q values and develop a methodology for estimating inter-station attenuation Q values from linear array of three stations (triplet). ❧ A technique is developed for estimating the non-diffuse characteristics of ambient seismic noise by computing correlations of neighboring power spectral values at different frequencies. This technique is applied to broadband stations in southern California and show distinct features for primary and secondary ocean microseismic peaks. The observed features can be explained by adding cross-frequency random components to the fully diffuse noise through numerical simulation. ❧ A method is proposed for estimating the error of stacked noise cross-spectrum from evenly spaced noise windows and deriving phase and amplitude errors for tomography purposes. The method is verified through numerical simulations and then applied to southern California regional network stations across San Andreas fault and dense linear array stations across the San Jacinto fault zone with different frequency bands. Outliers of the ambient seismic noise are removed using a statistical approach based on probability distributions of noise cross-spectrum to improve the measurements of phase and amplitude from the stacked cross-spectrum. ❧ A theory on the effects of inter-station and background attenuation structure on the expected noise cross-correlation amplitude is constructed to improve the understanding on the anelastic decay of the noise amplitude cross-spectra. Based on the theory, a methodology is proposed to estimated inter-station attenuation from linear array of three stations, which is validated by numerical simulations and applied to array data in regional network and dense linear array in southern California. The attenuation Q⁻¹ values with confidence intervals are estimated through bootstrap resampling by days and the results for a dense linear array suggest the location of a fault damage zone.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Multi-scale imaging and monitoring of crustal and fault zone structures in southern California
PDF
Detailed properties of seismic waveforms and earthquake sources in Southern California
PDF
Elements of seismic structures near major faults from the surface to the Moho
PDF
High-resolution imaging and monitoring of fault zones and shallow structures: case studies in southern California and on Mars
PDF
Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure
PDF
Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults
Asset Metadata
Creator
Liu, Xin
(author),
刘昕
(author)
Core Title
Statistical analyses of ambient seismic noise spectra with applications to detecting imperfectly diffuse wave field and deriving attenuation and phase velocity information
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Geological Sciences
Publication Date
08/01/2016
Defense Date
05/06/2016
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
ambient seismic noise,amplitude information,attenuation,diffuse wavefield,empirical Green's function,error estimation,OAI-PMH Harvest,seismic interferometry
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Ben-Zion, Yehuda (
committee chair
), Haas, Stephan (
committee member
), Sammis, Charles (
committee member
)
Creator Email
liu409@usc.edu,liuxine@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c40-289262
Unique identifier
UC11280594
Identifier
etd-LiuXin-4692.pdf (filename),usctheses-c40-289262 (legacy record id)
Legacy Identifier
etd-LiuXin-4692.pdf
Dmrecord
289262
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Liu, Xin; 刘昕
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
ambient seismic noise
amplitude information
attenuation
diffuse wavefield
empirical Green's function
error estimation
seismic interferometry