# Frequency domain analysis of errors in cross-correlations of ambient seismic noise

Xin Liu, Yehuda Ben-Zion, Dimitri Zigone

► To cite this version:

Xin Liu, Yehuda Ben-Zion, Dimitri Zigone. Frequency domain analysis of errors in cross-correlations of ambient seismic noise. Geophysical Journal International, 2016, 207 (3), pp.1630 - 1652. 10.1093/gji/ggw361. hal-01614797

HAL Id: hal-01614797

https://hal.science/hal-01614797

Submitted on 9 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Frequency domain analysis of errors in cross-correlations of ambient seismic noise

Xin Liu, Yehuda Ben-Zion and Dimitri Zigone*

Department of Earth Sciences, University of Southern California, Los Angeles, CA 90089-0740, USA. E-mail: liu409@usc.edu

Accepted 2016 September 23. Received 2016 September 21; in original form 2016 March 4

## SUMMARY

We analyse random errors (variances) in cross-correlations of ambient seismic noise in the frequency domain, which differ from previous time domain methods. Extending previous theoretical results on ensemble averaged cross-spectrum, we estimate confidence interval of stacked cross-spectrum of finite amount of data at each frequency using non-overlapping windows with fixed length. The extended theory also connects amplitude and phase variances with the variance of each complex spectrum value. Analysis of synthetic stationary ambient noise is used to estimate the confidence interval of stacked cross-spectrum obtained with different length of noise data corresponding to different number of evenly spaced windows of the same duration. This method allows estimating Signal/Noise Ratio (SNR) of noise cross-correlation in the frequency domain, without specifying filter bandwidth or signal/noise windows that are needed for time domain SNR estimations. Based on synthetic ambient noise data, we also compare the probability distributions, causal part amplitude and SNR of stacked cross-spectrum function using one-bit normalization or pre-whitening with those obtained without these pre-processing steps. Natural continuous noise records contain both ambient noise and small earthquakes that are inseparable from the noise with the existing pre-processing steps. Using probability distributions of random cross-spectrum values based on the theoretical results provides an effective way to exclude such small earthquakes, and additional data segments (outliers) contaminated by signals of different statistics (e.g. rain, cultural noise), from continuous noise waveforms. This technique is applied to constrain values and uncertainties of amplitude and phase velocity of stacked noise cross-spectrum at different frequencies, using data from southern California at both regional scale (~35 km) and dense linear array (~20 m) across the plate-boundary faults. A block bootstrap resampling method is used to account for temporal correlation of noise cross-spectrum at low frequencies (0.05–0.2 Hz) near the ocean microseismic peaks.

**Key words:** Time-series analysis; Interferometry; Seismic tomography; Theoretical seismology; Wave scattering and diffraction; Wave propagation.

## 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 & 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 data at two different locations sampling 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 & Weaver 2001; Roux et al. 2005; Sánchez-Sesma & 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 & Ben-Zion 2016) due to imperfectly scattered ocean microseismic noise or other effects. A non-diffuse ambient noise field may produce biased empirical Green's function estimated from cross-correlations (Liu & Ben-Zion 2016). Furthermore, the distribution of ambient seismic noise sources

<footnote>* Now at: Institut de Physique du Globe de Strasbourg, Université de Strasbourg, EOST, CNRS, Strasbourg, France.</footnote>
can be non-isotropic (e.g. Campillo 2006; Weaver et al. 2009; Tsai 2011), which can also bias empirical Green's function estimates. Despite those complications, ambient noise 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 et al. 2009; Lin et al. 2012; Liu et al. 2015), and body wave phases retrieval (e.g. Poli et al. 2012; Lin & Tsai 2013). Noise correlations were also used to analyse temporal changes of seismic velocities (e.g. Sens-Schönfelder & 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 homogeneous 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 & 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) generalized the variance expression to heterogeneous media and coloured noise spectra assuming isotropic noise source distribution and stationary 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 narrow-band 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. Sabra et al. 2005b; Lin et al. 2008; 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. These methods use only information from the stacked/averaged cross-correlation functions, and they do not consider temporal correlation structure and non-stationarity in the noise recordings. In this study we use more information and propose a practical approach for computing averaged cross-spectra and estimating random errors (variances) with validations from block bootstrap resampling.

In Section 2, we develop a probabilistic description of stacked (averaged) cross-spectrum based on assumptions on Gaussian noise model, stationary noise random process and N independent, identically distributed (i.i.d) random cross-spectrum observations. These can be approximated from N evenly spaced time windows extracted from ambient noise data. The variance of stacked cross-spectrum is found to be equal to the product of power spectrum values of two stations at the same frequency divided by the number of observations N (Section 2.1), which is linked to the time domain variance expression mentioned above. The random errors in phase and causal (or anti-causal) amplitude based on standard errors of complex stacked cross-spectrum are derived in Section 2.2. Equations necessary for practical data processing based on this method are also given in Section 2.2. The remaining theoretical sections may be useful for future work. The SNR for narrowband Gaussian filtered cross-correlation is estimated in time domain and linked to the variance and mean of stacked cross-spectrum (Section 2.3). The mean and variance estimations for piecewise stationary ambient noise records are given in Section 2.4.

In Section 3, we simulate synthetic noise data at two stations and study the distribution and statistics of stacked cross-spectrum at each frequency computed from numerous non-overlapping windows. The modified random cross-spectrum distributions and causal amplitude spectral decay curves based on two widely used pre-processing 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. In Section 4, similar error analysis procedures are applied to three regional stations (~35 km spacing) of the Southern California regional network crossing the San Andreas Fault. The histogram of real cross-spectrum observations 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 stacked cross-spectrum compared with the results based on standard error of stacked cross-spectrum. The results suggest that temporal correlations of cross-spectrum observations among different windows increase the variance of stacked cross-spectrum (between 0.05 and 0.2 Hz). In Section 5, we analyse the errors in high frequency noise recorded by densely spaced (~20 m) arrays that cross the San Jacinto Fault Zone. Removing the outliers increases the SNR, which improves both the amplitude and phase velocity estimations. Implications of the results and prospects for future work are discussed in Section 6.

## 2 THEORY

We consider ambient noise wave field at two stations a and b. The random wave field $d_j(\omega)$ at each station j ($j = a$ or b) can be expressed as the sum of correlated noise $u_j(\omega)$, which stands out through cross-correlation, and uncorrelated noise $v_j(\omega)$, which cancels out after cross-correlation, $d_j(\omega) = u_j(\omega) + v_j(\omega)$. Assuming the correlated noise field at different sites have common noise sources, the correlated noise spectrum of station j at angular frequency $\omega$ is written as a sum of contributions from numerous noise sources $s_k(\omega)$, which include scatterers (discrete version of eq. 1 in Liu & Ben-Zion 2013),

$$u_j(\omega) = A_j(\omega) \sum_k \frac{1}{\sqrt{\Delta_{jk}}} \exp \left[ -\frac{i\omega\Delta_{jk}}{\tilde{c}(\omega)} \right] s_k(\omega), \quad (1)$$

where $\Delta_{jk}$ represents distance from noise source k to station j. The inverse complex phase velocity is defined as $1/\tilde{c}(\omega) = (1 - i \operatorname{sgn}(\omega)/2Q(\omega))/c(\omega)$ where Q is the spatially varying (heterogeneous) attenuation quality factor and c(\omega) is the real phase velocity. $A_j(\omega)$ represents the site amplification factor for station j.

The noise sources $s_k(\omega)$ are assumed to be circular Gaussian random variables that are mutually uncorrelated such that $E[s_k^*(\omega)s_l(\omega)] = B_k(\omega)\delta_{kl}$, where $B_k(\omega)$ is the spectral density of kth noise source and $\delta_{kl}$ is the Kronecker delta. Other factors in eq. (1) depending on frequency $\omega$ can be absorbed in the spectral density term $B_k(\omega)$ for simplicity. The noise sources distribution can be non-isotropic. The noise sources are
assumed stationary and wave components at different frequencies are not correlated (Liu & Ben-Zion 2016). As a result of the summation in eq. (1), the correlated noise spectrum $u_j(\omega)$ is a Gaussian random variable. The uncorrelated noise spectrum is also assumed circular Gaussian random variable and satisfies $E[u_k^*(\omega)v_l(\omega)] = 0$ and $E[v_k^*(\omega)v_l(\omega)] = V_l(\omega)\delta_{kl}$, where $V_l(\omega)$ is the spectral density of uncorrelated noise at station l.

## 2.1 Stacked (averaged) cross-spectrum and its mean and variance

In frequency domain, the cross-correlation operation can be expressed as multiplication of the spectrums at two stations. In realistic cases, there are only finite lengths of noise recordings available. Here we assume i.i.d spectrum observations (random variables) at the same frequency, and denote the nth observation of random spectrum on station j as $d_j^n(\omega)$. The cross-spectrum between stations a and b from the nth observation is $d_a^{n*}(\omega)d_b^n(\omega)$. Then the average of N observed i.i.d cross-spectrum values between stations a and b can be formulated as

$$R_{ab}^{N}(\omega) = \frac{1}{N} \sum_{n=1}^{N} d_{a}^{n*}(\omega) d_{b}^{n}(\omega). \qquad (2)$$

The expectation of the averaged cross-spectrum is

$$E[R_{ab}^{N}(\omega)] = \frac{1}{N} \sum_{n=1}^{N} E[d_{a}^{n*}(\omega) d_{b}^{n}(\omega)] = \frac{1}{N} \sum_{n=1}^{N} E[u_{a}^{n*}(\omega) u_{b}^{n}(\omega)] = C_{ab}(\omega) \qquad (3)$$

where $C_{ab}(\omega)$ is the expected cross-spectrum (ensemble average) of correlated noise recordings from two stations a and b (e.g. Weaver et al. 2009; Liu & Ben-Zion 2013; Liu et al. 2015). The uncorrelated noise terms $v_a(\omega)$ and $v_b(\omega)$ cancel out.

The stacked (averaged) cross-spectrum $R_{ab}^N(\omega)$ is a complex random variable that approaches a Gaussian distribution for a large number of observations N according to the Central Limit Theorem of probability. The variance of $R_{ab}^N(\omega)$ is (Appendix A),

$$\mathrm{Var}[R_{ab}^{N}(\omega)] = \frac{1}{N} E[d_{a}^{n*}(\omega) d_{a}(\omega)] E[d_{b}^{n*}(\omega) d_{b}(\omega)] = \frac{1}{N} A_{aa}(\omega) A_{bb}(\omega) \qquad (4)$$

where $A_{aa}(\omega)$ and $A_{bb}(\omega)$ are expected (ensemble-averaged) power spectrum functions for stations a and b, respectively. This expression is related to the time domain variance expression (e.g. Sabra et al. 2005a), which is also proportional to the product of two autocorrelation functions divided by length of recording T (T is related to number of observations N; explained in the numerical simulation Section 3 using N windows).

The pseudo-variance of the stacked cross-spectrum is (Appendix A),

$$\mathrm{pVar}[R_{ab}^{N}(\omega)] = \frac{1}{N} E[d_{a}^{n*}(\omega) d_{b}(\omega)] E[d_{a}^{n*}(\omega) d_{b}(\omega)] = \frac{1}{N} C_{ab}(\omega) C_{ab}(\omega) \qquad (5)$$

which is generally not zero. A non-zero pseudo-variance indicates that the complex random variable $R_{ab}^N(\omega)$ is not circular (Ollila 2008) and therefore its real and imaginary parts are correlated (Appendix B).

## 2.2 Errors of amplitude and phase measurements

For N i.i.d cross-spectrum observations (from finite amount of noise data), the standard errors for both real and imaginary parts of the averaged cross-spectrum can be measured directly from the distributions of these cross-spectrum observations. The standard errors describe the uncertainties in the averaged cross-spectrum. The averaged cross-spectrum can be represented by a sum of its expected value $C_{ab}$ and a zero-mean noise term $n_C(\omega)$ which is a complex random variable,

$$R_{ab}^{N}(\omega) = C_{ab}(\omega) + n_{C}(\omega) = C_{ab}(\omega) + n_{R}(\omega) + in_{I}(\omega) \qquad (6)$$

where $n_R(\omega)$ and $n_I(\omega)$ are real and imaginary parts of the complex residual noise $n_C(\omega)$, which has the same variance as the averaged cross-spectrum $R_{ab}^N(\omega)$, and is related to both the correlated noise $u_j$ and uncorrelated noise term $v_j$. The noise term $n_C(\omega)$ 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 ($t > 0$) and anti-causal ($t < 0$) parts propagating in opposite directions with the same speed (simplified from eq. 1 of Liu et al. 2015),

$$C_{ab}(\omega) = \alpha(\omega) \exp \left[ -i \frac{\omega x}{c(\omega)} + i \frac{\pi}{4} \right] + \beta(\omega) \exp \left[ i \frac{\omega x}{c(\omega)} - i \frac{\pi}{4} \right] \qquad (7)$$

where $\alpha(\omega)$ and $\beta(\omega)$ are two real amplitude terms including attenuation, site amplification and far-field noise source intensity for causal and anticausal parts, respectively. A phase shift of $\pi/4$ is equivalent to multiplying with $\sqrt{i}$ which is due to the approximation of Bessel functions when the inter-station distance is greater than the wavelength (Liu et al. 2015). Phase shift due to imperfectly isotropic noise sources is not considered to simplify derivations. In eq. (7), the wave phase speed $c(\omega)$ is assumed finite and greater than zero.
Substituting eq. (7) into eq. (6), we have the stacked cross-spectrum

$$R_{ab}^{N}(\omega) = [\alpha(\omega) + \beta(\omega)] \cos\left[\frac{\omega x}{c(\omega)} - \frac{\pi}{4}\right] + n_{R}(\omega) + in_{I}(\omega) + i[\beta(\omega) - \alpha(\omega)] \sin\left[\frac{\omega x}{c(\omega)} - \frac{\pi}{4}\right] \quad (8)$$

where the real and imaginary parts have $\beta(\omega) + \alpha(\omega)$ and $\beta(\omega) - \alpha(\omega)$ as their expected envelope functions, respectively. Hilbert transform can be used to separate the causal and the anti-causal parts because the analytic form of the averaged cross-spectra function is equivalent to isolating the causal part ($t > 0$) of the averaged cross-correlation. Applying frequency domain Hilbert transforms to the real and imaginary parts of the stacked cross-spectrum respectively, we get their analytic forms

$$\tilde{R}_{R}^{N}(\omega) = \tilde{C}_{ab,R}(\omega) + \tilde{n}_{R}(\omega) = [\alpha(\omega) + \beta(\omega)] \exp\left[i \frac{\omega x}{c(\omega)} - i \frac{\pi}{4}\right] + \tilde{n}_{R}(\omega) \\ \tilde{R}_{I}^{N}(\omega) = \tilde{C}_{ab,I}(\omega) + \tilde{n}_{I}(\omega) = [\beta(\omega) - \alpha(\omega)] \exp\left[i \frac{\omega x}{c(\omega)} - i \frac{\pi}{4} - i \frac{\pi}{2}\right] + \tilde{n}_{I}(\omega) \quad (9)$$

where $\tilde{C}_{ab,R}(\omega)$ and $\tilde{C}_{ab,I}(\omega)$ 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 $\beta(\omega) \pm \alpha(\omega)$ are zero for $|t| \ge \min[x/c(\omega)]$ 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 $\tilde{n}_{R}(\omega)$ and $\tilde{n}_{I}(\omega)$, respectively

$$\tilde{n}_{R}(\omega) = n_{R}(\omega) + iHn_{R}(\omega) \\ \tilde{n}_{I}(\omega) = n_{I}(\omega) + iHn_{I}(\omega) \quad (10)$$

where $Hn_{R}(\omega)$ and $Hn_{I}(\omega)$ are, respectively, the Hilbert transforms of $n_{R}(\omega)$ and $n_{I}(\omega)$. The Hilbert transforms $Hn_{R}(\omega)$ and $Hn_{I}(\omega)$ are uncorrelated with $n_{R}(\omega)$ and $n_{I}(\omega)$, respectively (Appendix C).

The wave packets for causal and anticausal parts of the cross-spectrum can be separated based on eq. (9),

$$2\hat{\alpha}(\omega) \exp[i\hat{\phi}] = \tilde{R}_{R}^{N}(\omega) - i\tilde{R}_{I}^{N}(\omega) \\ 2\hat{\beta}(\omega) \exp[i\hat{\phi}] = \tilde{R}_{R}^{N}(\omega) + i\tilde{R}_{I}^{N}(\omega) \quad (11)$$

From eq. (11) we can derive the mean values and uncertainties for the amplitude estimators $\hat{\alpha}(\omega)$ and $\hat{\beta}(\omega)$ as well as the phase angle estimator $\hat{\phi} = \omega x/c(\omega) - \pi/4$ based on the right hand side of eq. (11) consisting of Hilbert transforms of the real and imaginary parts of averaged cross-spectrum. The amplitude uncertainties for causal and anticausal parts are, respectively (Appendix B)

$$\mathrm{Var}[\hat{\alpha}(\omega)] = [\sigma_R^2(\omega) + \sigma_I^2(\omega)]/2 \\ \mathrm{Var}[\hat{\beta}(\omega)] = [\sigma_R^2(\omega) + \sigma_I^2(\omega)]/2 \quad (12)$$

where $\sigma_R^2(\omega) = \mathrm{Var}[n_R(\omega)]$ and $\sigma_I^2(\omega) = \mathrm{Var}[n_I(\omega)]$ are variances for the real and imaginary parts of stacked cross-spectrum, respectively. Eq. (12) suggests that the causal and anticausal parts have the same variance. The frequency domain SNR on causal or anti-causal part cross-spectrum can be conveniently defined as mean over square root of variance: $\hat{\alpha}(\omega)/\sqrt{\mathrm{Var}[\hat{\alpha}(\omega)]}$ or $\hat{\beta}(\omega)/\sqrt{\mathrm{Var}[\hat{\beta}(\omega)]}$.

The phase uncertainty can be derived based on error propagation in the nonlinear case (Appendix B),

$$\hat{\sigma}_{\psi|R}^2(\omega) = \frac{\sigma_R^2(\omega)}{[\hat{\alpha}(\omega) + \hat{\beta}(\omega)]^2} \\ \hat{\sigma}_{\psi|\alpha}^2(\omega) = \frac{[\sigma_R^2(\omega) + \sigma_I^2(\omega)]/4}{[\hat{\alpha}(\omega)]^2} \\ \hat{\sigma}_{\psi|\beta}^2(\omega) = \frac{[\sigma_R^2(\omega) + \sigma_I^2(\omega)]/4}{[\hat{\beta}(\omega)]^2} \quad (13)$$

where $\hat{\sigma}_{\psi|R}^2(\omega)$, $\hat{\sigma}_{\psi|\alpha}^2(\omega)$ and $\hat{\sigma}_{\psi|\beta}^2(\omega)$ 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.
## 2.3 SNR on narrow-band time domain cross-correlation

Assuming phase velocity $c(\omega)$ and amplitude terms $\alpha(\omega)$ and $\beta(\omega)$ 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,

$$ \tilde{R}_{ab}(t) = F^{-1}\{\tilde{A}(\omega - \omega_0) R_{ab}^{N}(\omega)\} = F^{-1}\{\tilde{A}(\omega - \omega_0)(C_{ab}(\omega) + n_C(\omega))\} \quad (14) $$

where the Gaussian filter is $\tilde{A}(\omega) = \exp(-a^2\omega^2)$ with $a$ the filter width parameter and $a > 0$. Applying convolution theorem and assuming the unfiltered stacked cross-spectrum satisfy Hermitian symmetry, the time domain narrow-band analytic cross-correlation becomes,

$$ \begin{aligned} \tilde{R}_{ab}(t) &= F^{-1}\{\tilde{A}(\omega - \omega_0)\}_{*t} F^{-1}\{C_{ab}(\omega) + n_C(\omega)\}_{\omega_0} \\ &= \exp(i\omega_0 t) A(t)_{*t} \left\{ \alpha(\omega_0) \delta\left(t - \frac{x}{c(\omega_0)}\right) \exp\left[i\frac{\pi}{4}\right] + \beta(\omega_0) \delta\left(t + \frac{x}{c(\omega_0)}\right) \exp\left[-i\frac{\pi}{4}\right] + n_C(t; \omega_0) \right\}, \end{aligned} \quad (15) $$

where $A(t) = \exp(-t^2/(4a^2))/(2a\sqrt{\pi})$ is the inverse Fourier transform of the filter $\tilde{A}(\omega)$, 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 $\omega_0$ is equal to $\alpha(\omega_0)/(2a\sqrt{\pi})$, and similarly for the anti-causal part the peak amplitude is $\beta(\omega_0)/(2a\sqrt{\pi})$.

The RMS estimation of noise in the real time domain cross-correlation $Re[\tilde{R}_{ab}(t)]$ is (Appendix D),

$$ \text{RMS}_{n,\text{TD}}(\omega_0) = \sqrt{\frac{1}{4a\sqrt{2\pi}}} E[n_C(\omega_0) n_C^*(\omega_0)] = \frac{1}{2} \sqrt{\frac{\sigma_R^2(\omega_0) + \sigma_I^2(\omega_0)}{a\sqrt{2\pi}}} \quad (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 eqs (15) and (16), the conventional SNR measurement in time domain on causal part can be related to (Appendix D),

$$ \text{SNR}_{\alpha}(\omega_0) = \sqrt{\frac{1}{a} \frac{\sqrt{2}}{\pi} \frac{\alpha(\omega_0)}{\sqrt{\sigma_R^2(\omega_0) + \sigma_I^2(\omega_0)}}} \quad (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.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 $\omega$ of jth day are,

$$ \begin{aligned} m_j &= E[X_j] \\ \sigma_j^2 &= \text{Var}[X_j]. \end{aligned} \quad (18) $$

They can be estimated from the sample mean and variance of observations 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 observations and the standard error from jth day is, respectively,

$$ \begin{aligned} \bar{X}_j &= \frac{1}{M} \sum_{k=1}^{M} X_j(t_k) = m_j \\ \sigma_{j,\text{stde}}^2 &= \frac{1}{M} \text{Var}[X_j] = \frac{1}{M} \sigma_j^2 \end{aligned} \quad (19) $$

where $t_k$ is the time of the kth time window in a day.

Let $X_S$ be the averaged stack of $N$ days of cross-spectrum,

$$ X_S = \frac{1}{N} \sum_{j=1}^{N} \bar{X}_j. \quad (20) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. (a) Geometry for numerical simulation of ambient noise at two stations. The interstation distance is 50 km and the interstation attenuation is 30. The background attenuation is 500. The far field noise source distribution is uniform in all directions and source-receiver distance is 8°. (b) Schematic plot for evenly spaced windows on noise recording. Window length is 100 s and gap is 20 s. (c) 3600 independent cross-spectrum curves computed from 3600 windows in 5 d. (d) Stacked cross-spectra with error bars from 3600 observations (windows). (e) Histograms of the real part and joint real & imaginary parts cross-spectrum from 3600 observations. (f) Stacked cross-correlation of 3600 windows from 5 d of data.</figcaption></figure>

Then the mean and variance of $X_S$ are, respectively

$$m_S = E[X_S] = \frac{1}{N} \sum_{j=1}^{N} m_j$$ $$\sigma_S^2 = \mathrm{Var}[X_S] = \frac{1}{N^2} \sum_{j=1}^{N} \sigma_{j,stde}^2 = \frac{1}{NM} \left( \frac{1}{N} \sum_{j=1}^{N} \sigma_j^2 \right) = \frac{1}{NM} \bar{\sigma}^2. \qquad (21)$$

where $\bar{\sigma}^2 = \sum_{j=1}^{N} \sigma_j^2/N$ is the average variance of N different days and NM is the total number of 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 data variance $\bar{\sigma}^2$ divided by total number of observations NM. There are two ways to reduce the variance of stacked cross-spectrum $X_S$: (1) increase the total number of observations NM or (2) remove or rescale the daily segments of very large or small variances in order to minimize $\bar{\sigma}^2$ relative to the absolute value of the mean $m_S$.

## 3 NUMERICAL SIMULATIONS

We perform numerical simulations of finite-length (25 d) stationary random noise records at two stations (Fig. 1a) and compute stacked (averaged) cross-spectrum of the two stations. The standard errors of the real and imaginary parts of each stacked cross-spectrum value are computed and the errors on amplitude and phase are estimated based on the stacked cross-spectrum variances.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** 1σ phase random error and errors for its derived quantities. (a) Phase error computed from the symmetric component (real part) of cross-spectra in Fig. 1(d). (b) Phase traveltime error. (c) Phase velocity dispersion curve with 1σ confidence interval. (d) Relative error in phase velocity.</figcaption></figure>

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\omega)$, which is exaggerated compared with the realistic case to illustrate how the uncorrelated noise 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° away from the centre 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 d 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 variance 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 observations 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 variances, we select 5-d long noise data with 3600 non-overlapping windows, which produce 3600 independent observations of random cross-spectra curves (Fig. 1c). The stacked cross-spectra of the 3600 independent cross-spectra curves are shown in Fig. 1(d) with error bars for real and imaginary parts. The standard errors for the real and imaginary parts at the same frequency are very similar but the amplitude of the real part is much higher than that on the imaginary part due to the isotopic noise source distribution. The real part and joint real-imaginary parts histograms of all observations 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 well constrained because it is 10 times more than its standard error 1.92e-6. The stacked cross-correlation function of the 5-d (3600 windows) data is shown in Fig. (1f).

The phase random error (Fig. 2a) can be derived from eq. (13) using standard errors from the stacked cross-spectrum variances. 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 $\pi/4$, otherwise it may cause phase unwrapping ambiguities for traveltime and phase velocity estimates. Assuming the phase error is much smaller than $\pi$, the error in traveltime can be conveniently derived from the linear relation $dt = \omega d\varphi$ (Fig. 2b) and can be used for tomography based on traveltime measurements on surface waves (e.g. Lin et al. 2009). In addition, assuming that the ray path is straight, phase velocity random errors (Figs 2c and d) can be derived from $dc = -c^2d\varphi/\omega 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 random errors can potentially be used for tomography based on straight ray path approximation and monitoring temporal velocity changes at different frequencies.

The effects of several commonly used pre-processing steps on distributions of cross-spectrum observations and amplitude of stacked cross-spectrum are evaluated using the 25-d long noise data with 18 000 windows. Figs 3(a)–(c) show results for no pre-processing, one-bit normalization and pre-whitening, respectively. The top panels contain the joint distributions of the real and imaginary parts of cross-spectrum
<figure><img src="image_5.png" /><figcaption>Figure 3. Comparison among different pre-processing techniques: no pre-processing (a), one-bit normalization (b) and pre-whitening (c). 18 000 windows (25 d) 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-spectrum amplitude curves with error bars.</figcaption></figure>

at 0.14 Hz. The bottom panels show the causal part of stacked cross-spectra amplitude decay curves computed from first line of eq. (11) with the different pre-processing steps previously mentioned. These results suggest that both one-bit and pre-whitening modify the distribution of cross-spectrum observations and produce larger random errors relative to their corresponding amplitude mean values. Moreover, the pre-whitening step also normalizes the absolute square of each observed cross-spectrum value to one (Fig. 3c)

$$R_{ab,\text{whiten}}^{N}(\omega) = \frac{1}{N} \sum_{n=1}^{N} \frac{d_{a}^{n*}(\omega) d_{b}^{n}(\omega)}{|d_{a}^{n*}(\omega)| |d_{b}^{n}(\omega)|} \quad (22)$$

Therefore, the joint distribution of real and imaginary parts with pre-whitening is non-zero only on the unit circle in the complex plane. The pre-whitening step significantly modifies the amplitude decay curve (Fig. 3f) and does not allow for accurate amplitude estimations.

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 coda noise. The narrow-band Gaussian filter width is $a = 28.29$ and the noise window is selected between 50 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 observations (windows) are used for each frequency. For a fixed number of observations, 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 relative to correlated noise power and consequently slower convergence rate according to eq. (4). The time domain SNR is unstable and fluctuates with increasing trend as the number of observations increases while the frequency domain SNR increases steadily with more observations (Figs 4a and b). 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 d (18 000 windows) are shown in Figs 4(c) and (d) 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 observations, however, the maximum SNR values (14.5 for one bit and 15 for pre-whitening) are much lower than in the case without pre-processing (Fig. 4b, maximum SNR is 18) because of the nonlinear nature of the one-bit or pre-whitening pre-processing steps.

## 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 d in the year 2014 (days 50–300) are divided into 4-hour long segments. Segments dominated by significant earthquakes are removed
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Causal part SNR convergence images for 25 d/18 000 windows noise data. (a) and (b) are SNR convergence images of causal part cross-correlation amplitude estimated in time and frequency domains, respectively. SNR is defined as amplitude divided by corresponding square-root of variance. (c) and (d) are stacked cross-correlations from 18 000 windows (observations) for one-bit normalization and pre-whitening, respectively. (e) Frequency domain SNR convergence image of causal part cross-correlation amplitude with one-bit normalization. (f) Frequency domain SNR convergence image of causal part cross-correlation amplitude with pre-whitening.</figcaption></figure>

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 pre-processing step to preserve the original amplitude information. A band pass filter is applied between 0.05 and 0.6 Hz. The cross-correlation functions of raw 250-d stack for the three station pairs are shown in Fig. 5(b).

The analysis of random errors primarily focuses on the station pair CHF-SBB2. As in the numerical simulations, the non-overlapping window length is 100 s with gap of 20 s between windows (see sketch in Fig. 1b). The 250-d data are equivalent to 162 250 windows (observations) after earthquake segment removal. The histogram of the real part cross-spectrum at 0.24 Hz is shown in Fig. 6(a) and suggests the existence of many outlier signals with different statistics creating the long tails of the histogram. The observations outside of three Median Absolute Deviation (MAD) are assumed to be outliers and removed. After outlier removal, the new distribution of 154 137 remaining observations 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 from stacked cross-spectrum. 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).
<figure><img src="image_7.png" /><figcaption>**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 colours show low and high elevation respectively. (b) Cross-correlations of three station pairs computed from 250 d in 2014.</figcaption></figure>

Outliers for each frequency can be identified based on MAD in a similar way as in Figs 6(a) and (b). The identified outliers for each frequency can then be associated with their respective noise window. For each window, the percentage of outliers over all frequencies in the frequency passband of interest is an indicator of the amount of abnormal signal within the window. We select windows with maximum 5 per cent outliers over all frequencies between 0.05 and 0.6 Hz for stacking. The 115 661 selected windows produce the same number of observations of random cross-spectra curves (Fig. 6d) while the 46 589 unselected (outlier) windows give random cross-spectra curves of different statistics (Fig. 6e). A simple stack of the selected windows produces cross-spectra with associated standard errors computed from the distribution of cross-spectrum at each frequency (Fig. 6f). Similarly, Fig. 6(g) shows the stacked cross-spectra for outlier windows, which presents 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 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.

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/windowing 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 stack (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 (Figs 7a and b, respectively) are computed in the same way as in the synthetic test. The colour 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 similar peak frequencies of SNR values, suggesting multiple noise sources at different frequencies other than the main ocean microseism peak. However, the relative amplitudes of the peaks are different for time and frequency domain SNR figures.

By slicing both SNR images at three different frequencies (0.15, 0.20 and 0.32 Hz), the corresponding SNR values are plotted as functions of number of days in Figs 7(c) and (d) for time and frequency domain SNR, respectively. At 0.15 and 0.32 Hz, the SNR values increase with number of days because they correspond to the noise peaks in Figs 7(a) and (b). At 0.20 Hz, however, the SNR stops increasing beyond 35 d probably because of a gap between two noise sources in which the noise random process is highly non-stationary. By slicing the SNR images in Figs 7(a) and (b) at three different day counts, we plot in Figs 7(e) and (f) SNR functions of centre frequency for time and frequency domain SNR measurements, respectively. For 5 d 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 d is only 40). For 125 and 250 d 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.3, 0.52 Hz) the SNR values keep increasing up to 250 d. For frequency domain SNR (Fig. 7f), all three SNR curves for 5, 125 and 250 d 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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Raw distribution of 162 250 real part cross-spectrum observations from same number of windows at 0.24 Hz for station pair CHF-SBB2. (b) Distribution of real part cross-spectrum (154 137 observations) 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 115 661 selected windows. (e) First 2000 cross-spectra curves from 46 589 outlier (unselected) windows. (f) Stacked cross-spectra of CHF-SBB2 for selected windows (115 661 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 (46 589 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 unit variance.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) and (b) are SNR convergence images of causal part cross-correlation amplitude for 250 d estimated in time and frequency domains, respectively. Only selected windows (defined in Fig. 6f) are used. SNR is defined as amplitude divided by corresponding square-root of variance. (c) and (d) are SNR convergence curves sliced at different 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 d) for time and frequency domain measurements, respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Block bootstrap (hourly block) estimation of empirical confidence intervals of 250-d stacked cross-spectrum for CHF-SBB2. (a) Temporal correlation of noise revealed by autocorrelation of real 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) Histograms of bootstrap averaged real part cross-spectrum at 0.16 Hz and 0.24 Hz, respectively. The bootstrap standard deviation provides estimation for the uncertainty of stacked cross-spectra. (c) Ratio of real/imaginary part block bootstrap standard deviation (more realistic) over standard error (assuming i.i.d observations from different windows). Bootstrap method suggests larger variance between 0.05 and 0.2 Hz due to temporal correlation. The ratio approaches to one for frequency between 0.2 and 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 and 0.2 Hz are lower than the frequency domain SNR based on standard error (Fig. 7b).</figcaption></figure>

from frequency domain method probably because the frequency domain method uses variance of cross-spectrum in addition to mean (the stacked) cross-spectrum.

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. 8(a), 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 (49 windows offset with 100 s window length and 20 s gap). 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 and 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 (115 661) for pair CHF-SBB2 as in the previous SNR example. The purpose of hourly block bootstrap is to capture some temporal correlation among the neighbouring windows if they exist. For each one-hour block, cross-spectra curves from all selected windows within the block are averaged. Bootstrap resampling is applied by randomly pick the same number of blocks from available 1 hr blocks with replacement. At each frequency, a bootstrap mean cross-spectrum observation is computed by averaging randomly drawn block cross-spectrum values. There are 4000 bootstrap averaged cross-spectrum observations (Fig. 8b) and the standard deviation of these observations approximates the uncertainty of the stacked cross-spectrum. The bootstrap histograms (Fig. 8b) approach to Gaussian distribution that confirm with the statistical property of stacked cross-spectrum in theory.

The ratio of 250-d block bootstrap standard deviation over standard error (which implicitly assumes i.i.d windows) for the real and imaginary parts of cross-spectra are plotted in Fig. 8(d). 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 *N* i.i.d cross-spectrum observations at the same frequency from different windows offers reliable estimation of the data variance between 0.2 and 0.6 Hz. The ratio for the real and imaginary parts between 0.05 and 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).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Phase and amplitude random errors ($1\sigma$) of cross-spectra based on block bootstrap resampling of 11 5661 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 traveltime error. (c) Phase velocity dispersion curve with $1\sigma$ 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).</figcaption></figure>

Applying block bootstrap to different number of days (1–250), the mean and standard deviation of bootstrap distribution for stacked cross-spectrum value can be easily computed with increasing number of days 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 observations at each frequency from different windows. For frequency between 0.05 and 0.2 Hz, however, the SNR values derived from bootstrap are evidently lower (~20 per cent 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 better capture the temporal correlation at this frequency and therefore further increase the variance and reduce the corresponding SNR value.

The phase error curve (Fig. 9a) is derived directly from the bootstrap variance of the real part of stacked cross-spectrum according to eq. (13). The maximum phase error ($1\sigma$) is 0.11 radians, which is much less than phase angle period $2\pi$. This indicates that the phase angle can be reliably unwrapped without $2\pi$ ambiguities. Following similar procedure as in numerical simulation (Section 3), the phase traveltime uncertainty can be estimated directly from the phase error. Assuming straight ray path, the phase angle errors can be mapped into phase velocity uncertainties (Figs 9c and d) and the largest uncertainty is 0.4 per cent velocity variation within $1\sigma$ interval. The causal part amplitude cross-spectrum (Figs 9e and 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.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Random error analysis using 1-d 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 second 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 observations from same number of windows in 1 d for station pair JF00-JFS1. (e) Distribution of real part cross-spectrum (7185 observations) for JF00-JFS1 after outlier removal. (f) First 2000 of 6220 selected cross-spectrum curves from 6220 windows. (g) Raw stack of cross-spectrum curves computed from all 7332 windows in the second day of 2012. (h) Stacked cross-spectra of 6220 selected cross-spectra curves (panel e) from same number of windows after removing outliers.</figcaption></figure>

## 5 APPLICATION TO DENSELY SPACED STATIONS WITH HIGH FREQUENCY DATA

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 separated by about 20 m (Fig. 10a). The cross-correlation functions of first 10 d in the year 2012 are plotted in Fig. 10(b). 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 1 d (the second 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 observations 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 observations (windows), the new distribution (Fig. 10d) of real part cross-spectrum with 7185 observations has a higher ratio of mean over standard error of 13.14. A selected window should contain no more than 5 per cent 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. 10g) have more fluctuations and larger error bars relative to spectral
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Amplitude SNR and phase random error based on 1-d 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 traveltime error derived based phase error. (e). Phase velocity curves. (f) Phase velocity relative error.</figcaption></figure>

amplitude than the stacked cross-spectra after outlier removal (Fig. 10h). There is an anomalous peak at 20 Hz in Fig. (10h), which might be related to some monochromatic local noise source at the same frequency.

Figs 11(a) and (b) compare the amplitude SNR values based on raw stacked (with outliers) cross-spectra and outlier-removed cross-spectra obtained from causal part (and anticausal part) of 1-d stacked cross-spectra for pair JF00-JFS1. After outlier removal, the SNR of causal part amplitude is 4–10 times greater than that before outlier removal between 17 and 45 Hz (Fig. 11a). For anticausal part amplitude, the outlier-removed SNR is 4–10 times better than outlier-included SNR between 10 and 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.

The outliers also affect phase velocity retrieved from 1-d stack of cross-spectra for the pair JF00-JFS1. The phase random 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 ambiguities. Based on phase error information, the phase traveltime error (Fig. 11d) and phase velocity error (Fig. 11f) can be derived following the same way as in numerical simulation (Section 3). The phase velocity dispersion curves based on stacked cross-spectra with and without outliers show 1–4 per cent discrepancy between 11 and 20 Hz (Fig. 11e), which is probably due to systematic errors brought by outliers across different frequencies.

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 12a and b) compared to the phase velocity errors from the second 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 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 and d) and JFS1-JFS2 (Figs 12e and f) show improved phase error estimations and less spikes after outlier removal.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Phase velocity and relative phase velocity random error based on 1-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) and (b) are phase velocity and corresponding relative error on JF00-JFS1. (c) and (d) are phase velocity and corresponding relative error on JF00-JFS2. (e) and (f) are phase velocity and corresponding relative error on JFS1-JFS2.</figcaption></figure>

# 6 DISCUSSION

Analysis of random errors in ambient noise cross-correlation can provide important information on the deviation of the stacked noise cross-correlation of finite amount of data from the expected cross-correlation, which is linked through time derivative to the empirical Green's function. We develop a formulation to estimate the variances 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 traveltime and phase velocity random 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 stacked noise cross-correlation, while the frequency domain noise variance can be derived directly from the standard error of numerous i.i.d cross-spectrum observations at the same frequency in different windows, which uses more information from the observations. The variance of stacked cross-spectrum at each frequency is easier to derive and can provide more information than the time domain variance, as it is a statistical quantity based on the distribution of i.i.d cross-spectrum observations and can be validated using power spectrum functions in the eq. (4).
At the same frequency, the cross-spectrum observations 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.

Numerical simulation with 25-d stationary noise produced from Gaussian noise model are used to study statistical properties of stacked cross-spectrum under ideal conditions. As can be inferred from the histograms of real and imaginary parts cross-spectrum observations 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 observations, the distribution remains the same (due to the stationary property) but the standard error decreases and the distribution of stacked cross-spectrum approaches to Gaussian.

Effects of two pre-processing techniques, one-bit normalization and pre-whitening, are compared to the case without pre-processing normalization in terms of distribution of random cross-spectrum observations and causal part amplitude decay. The results show that both pre-processing methods modify the distribution of cross-spectrum values and increase the standard errors by ~15 per cent (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 attenuation estimations based on amplitude information. Together, these results highlight the strong influence of the pre-processing procedures on the estimation of any amplitude related quantities. The proposed methodology allows assessing the discrepancy of amplitude estimations between different pre-processing methods with confidence intervals. Such information can potentially 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 numerical simulation, and the time domain SNR is non-unique depending on the narrowband filter width and coda noise window definition (e.g. Sabra et al. 2005a,b; Zigone et al. 2015).

The distribution of cross-spectrum of real data contains outliers that could 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 various activities leading to outliers that can bias the estimations of empirical Green's function from noise cross-correlation.

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 that the SNR should be relatively lower than predicted by the frequency domain SNR image because of the temporal correlation in ambient seismic noise data 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 and 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 and 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 to capture this temporal correlation for more accurate variance estimation at those frequencies. For cross-spectrum between 0.2 and 0.6 Hz, the temporal correlation is close to zero (<0.01) and the bootstrap standard deviation converges to the standard error that assumes i.i.d cross-spectrum observations from different windows. The frequency band of temporal correlation coincides with correlated neighbouring frequency band (0.07–0.2 Hz) near the secondary ocean microseismic peak (Liu & Ben-Zion 2016), which contains non-diffuse noise components. Another limitation of this paper is that it only discusses random errors of statistical quantity (e.g. the variance of stacked cross-spectrum) with respect to its expectation. We do not estimate any systematic bias of the expectation of a statistical estimator from the true physical quantity (e.g. the analysis does not consider small phase shift in the stacked cross-spectrum due to strong directional far-field noise sources).

Based on the bootstrap cross-spectrum variances, the phase velocity uncertainty for 250-d stack on selected windows is less than 0.23 per cent for frequency between 0.15 and 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 uncertainty to claim temporal change of phase velocity at certain frequency. The phase uncertainty 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 the 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 1-d stack suggest that after outlier removal, the amplitude SNRs of stacked cross-spectra are improved by a factor of 4–10 and phase velocity uncertainties 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 data 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 confidence intervals before outlier removal are ~2–6 times greater than those after removing outliers.

The presented error analysis method is promising for analysing the variances of stacked noise cross-correlation data at various scales and frequencies with applications to amplitude (attenuation/site amplification) and phase (velocity/traveltime) results. 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 issue and additional related work are left for a future study.

ACKNOWLEDGEMENTS

We thank Robert A. Scholtz and Zachary Ross for discussions on time-frequency duality, non-stationary noise and temporal correlations. We also thank Toshiro Tanimoto, and anonymous referee and editor Michael Ritzwoller for useful comments that helped us to improve the paper. The study was supported by the Southern California Earthquake Center (based on National Science Foundation Cooperative Agreement EAR-0529922 and U.S. Geological Survey Cooperative Agreement 07HQAC0008) and the Department of Energy (award DE-SC0016520).

REFERENCES

Aki, K., 1957. Space and time spectra of stationary stochastic waves, with special reference to microtremors, *Bull. Earthq. Res. Inst.*, **35**, 415–457.

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys. J. Int.*, **169**(3), 1239–1260.

Ben-Zion, Y. et al., 2015. Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone, *Geophys. J. Int.*, **202**(1), 370–380.

Boashash, B., 1992. Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals, *Proc. IEEE*, **80**(4), 520–538.

Brenguier, F., Campillo, M., Hadzioinnou, C., Shapiro, N.M., Nadeau, R.M. & Larose, E., 2008. Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations, *Science*, **321**(5895), 1478–1481.

Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N.M., Briand, X., Emoto, K. & Miyake, H., 2014. Mapping pressurized volcanic fluids from induced crustal seismic velocity drops, *Science*, **345**(6192), 80–82.

Campillo, M., 2006. Phase and correlation in ‘random’ seismic fields and the reconstruction of the green function, *Pure appl. Geophys.*, **163**(2–3), 475–502.

Gouédard, P. et al., 2008. Cross-correlation of random fields: mathematical approach and applications, *Geophys. Prospect.*, **56**(3), 375–393.

Hillers, G., Campillo, M., Ben-Zion, Y. & Roux, P., 2014. Seismic fault zone trapped noise, *J. geophys. Res.*, **119**, 5786–5799.

Hillers, G., Ben-Zion, Y., Campillo, M. & Zigone, D., 2015. Seasonal variations of seismic velocities in the San Jacinto Fault area observed with ambient seismic noise, *Geophys. J. Int.*, **202**, 920–932.

Lin, F.-C. & Tsai, V.C., 2013. Seismic interferometry with antipodal station pairs: antipodal seismic interferometry, *Geophys. Res. Lett.*, **40**(17), 4609–4613.

Lin, F.-C., Moschetti, M.P. & Ritzwoller, M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, *Geophys. J. Int.*, **173**(1), 281–298.

Lin, F.-C., Ritzwoller, M.H. & Snieder, R., 2009. Eikonal tomography: surface wave tomography by phase front tracking across a regional broadband seismic array, *Geophys. J. Int.*, **177**(3), 1091–1110.

Lin, F.-C., Tsai, V.C. & Ritzwoller, M.H., 2012. The local amplification of surface waves: a new observable to constrain elastic velocities, density, and anelastic attenuation, *J. geophys. Res.*, **117**(B6), doi:10.1029/2012JB009208.

Liu, X. & Ben-Zion, Y., 2013. Theoretical and numerical results on effects of attenuation on correlation functions of ambient seismic noise, *Geophys. J. Int.*, **194**(3), 1966–1983.

Liu, X. & Ben-Zion, Y., 2016. Estimating correlations of neighbouring frequencies in ambient seismic noise, *Geophys. J. Int.*, **206**(2), 1065–1075.

Liu, X., Ben-Zion, Y. & Zigone, D., 2015. Extracting seismic attenuation coefficients from cross-correlations of ambient noise at linear triplets of stations, *Geophys. J. Int.*, **203**(2), 1149–1163.

Lobkis, O.I. & Weaver, R.L., 2001. On the emergence of the Green's function in the correlations of a diffuse field, *J. acoust. Soc. Am.*, **110**(6), 3011, doi:10.1121/1.1417528.

Ollila, E., 2008. On the circularity of a complex random variable, *IEEE Signal Process. Lett.*, **15**, 841–844.

Poli, P., Campillo, M. & Pedersen, H. & LAPNET Working Group, 2012. Body-wave imaging of Earth's mantle discontinuities from ambient seismic noise, *Science*, **338**(6110), 1063–1065.

Poli, P., Pedersen, H.A. & Campillo, M. & the POLENET/LAPNET Working Group, 2013. Noise directivity and group velocity tomography in a region with small velocity contrasts: the northern Baltic shield, *Geophys. J. Int.*, **192**(1), 413–424.

Prieto, G.A., Lawrence, J.F. & Beroza, G.C., 2009. Anelastic Earth structure from the coherency of the ambient seismic field, *J. geophys. Res.*, **114**(B7), doi:10.1029/2008JB006067.

Roux, P., Sabra, K.G., Gerstoft, P., Kuperman, W.A. & Fehler, M.C., 2005. P-waves from cross-correlation of seismic noise, *Geophys. Res. Lett.*, **32**(19), L19303, doi:10.1029/2005GL023803.

Sabra, K.G., Roux, P. & Kuperman, W.A., 2005a. Emergence rate of the time-domain Green's function from the ambient noise cross-correlation function, *J. acoust. Soc. Am.*, **118**(6), 3524, doi:10.1121/1.2109059.

Sabra, K.G., Gerstoft, P., Fehler, M.C., Roux, P. & Kuperman, W.A., 2005b. Extracting time-domain Green's function estimates from ambient seismic noise, *Geophys. Res. Lett.*, **32**(3), doi:10.1029/2004GL021862.

Sánchez-Sesma, F.J. & Campillo, M., 2006. Retrieval of the Green's function from cross correlation: the canonical elastic problem, *Bull. seism. Soc. Am.*, **96**(3), 1182–1191.

Seats, K.J., Lawrence, J.F. & Prieto, G.A., 2012. Improved ambient noise correlation functions using Welch's method, *Geophys. J. Int.*, **188**(2), 513–523.

Sens-Schönfelder, C. & Wegler, U., 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia, *Geophys. Res. Lett.*, **33**(21), doi:10.1029/2006GL027797.

Sens-Schönfelder, C., Snieder, R. & Stähler, S.C., 2015. The lack of equipartitioning in global body wave coda, *Geophys. Res. Lett.*, **42**(18), 7483–7489.

Shapiro, N.M. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.*, **31**(7), L07614, doi:10.1029/2004GL019491.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**(5715), 1615–1618.

Snieder, R., 2004. Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase. *Phys. Rev. E*, **69**(4), 046610, doi:10.1103/PhysRevE.69.046610.
Tsai, V.C., 2011. Understanding the amplitudes of noise correlation measurements, *J. geophys. Res.*, **116**(B9), B09311, doi:10.1029/2011JB008483.

Weaver, R., Froment, B. & Campillo, M., 2009. On the correlation of non-isotropically distributed ballistic scalar diffuse waves, *J. acoust. Soc. Am.*, **126**, 1817–1826.

Weaver, R.L. & Lobkis, O.I., 2005. Fluctuations in diffuse field–field correlations and the emergence of the Green’s function in open systems, *J. acoust. Soc. Am.*, **117**(6), 3432, doi:10.1121/1.1898683.

Zigone, D., Ben-Zion, Y., Campillo, M. & Roux, P., 2015. Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and Love Waves, *Pure appl. Geophys.*, **172**, 1007–1032.

## APPENDIX A: DERIVATION OF VARIANCE AND PSEUDO-VARIANCE OF STACKED CROSS-SPECTRUM

The variance of the stacked cross-spectrum $R_{ab}^N(\omega)$ is by definition

$$ \mathrm{Var}[R_{ab}^N(\omega)] = E[R_{ab}^N(\omega) R_{ab}^{N*}(\omega)] - E[R_{ab}^N(\omega)] E[R_{ab}^{N*}(\omega)]. \quad (A1) $$

Combined with eq. (2), the correlation of stacked cross-spectrum is

$$ \begin{aligned} E[R_{ab}^N(\omega) R_{ab}^{N*}(\omega)] &= E\left[\frac{1}{N} \sum_{n=1}^N d_a^{n*}(\omega) d_b^n(\omega) \frac{1}{N} \sum_{p=1}^N d_a^p(\omega) d_b^{p*}(\omega)\right] \\ &= \frac{1}{N^2} \sum_{n=1}^N E[d_a^{n*}(\omega) d_b^n(\omega)] \sum_{p=1}^N E[d_a^p(\omega) d_b^{p*}(\omega)] \\ &\quad + \frac{1}{N^2} \sum_{n=1}^N E[d_a^{n*}(\omega) d_a^n(\omega)] E[d_b^n(\omega) d_b^{n*}(\omega)]. \end{aligned} \quad (A2) $$

Normally, the expectation of the product of four Gaussian random variables in eq. (A2) is equal to the sum of three terms of second-order moments of Gaussian processes. One of the three terms regarding the convolution of noise recordings on stations a and b is equal to zero.

Substituting eq. (A2) back into eq. (A1), we get the variance of the stacked cross-spectrum:

$$ \mathrm{Var}[R_{ab}^N(\omega)] = \frac{N^2}{N^2} C_{ab}(\omega) C_{ab}^*(\omega) + \frac{1}{N} A_{aa} A_{bb} - C_{ab}(\omega) C_{ab}^*(\omega) = \frac{1}{N} A_{aa} A_{bb}, \quad (A3) $$

where the expected power spectrum on station a is $A_{aa} = E[d_a^*(\omega)d_a(\omega)]$ and it takes a similar form on station b.

The pseudo-variance of the stacked cross-spectrum is defined as

$$ \mathrm{pVar}[R_{ab}^N(\omega)] = E[R_{ab}^N(\omega) R_{ab}^N(\omega)] - E[R_{ab}^N(\omega)]^2 \quad (A4) $$

and similarly, when combined with eq. (2), the pseudo-correlation is

$$ \begin{aligned} E[R_{ab}^N(\omega) R_{ab}^N(\omega)] &= E\left[\frac{1}{N} \sum_{n=1}^N d_a^{n*}(\omega) d_b^n(\omega) \frac{1}{N} \sum_{p=1}^N d_a^{p*}(\omega) d_b^p(\omega)\right] \\ &= \frac{1}{N^2} \sum_{n=1}^N E[d_a^{n*}(\omega) d_b^n(\omega)] \sum_{p=1}^N E[d_a^{p*}(\omega) d_b^p(\omega)] \\ &\quad + \frac{1}{N^2} \sum_{n=1}^N E[d_a^{n*}(\omega) d_b^n(\omega)] E[d_a^{n*}(\omega) d_b^n(\omega)]. \end{aligned} \quad (A5) $$

Substituting eq. (A5) back into eq. (A4), we get the pseudo-variance of the stacked cross-spectrum:

$$ \mathrm{pVar}[R_{ab}^N(\omega)] = \frac{N^2}{N^2} C_{ab}(\omega) C_{ab}(\omega) + \frac{1}{N} C_{ab}(\omega) C_{ab}(\omega) - C_{ab}(\omega) C_{ab}(\omega) = \frac{1}{N} C_{ab}(\omega) C_{ab}(\omega). \quad (A6) $$

## APPENDIX B: DERIVATIONS FOR ERRORS OF AMPLITUDE AND PHASE FROM THE SECOND MOMENT STATISTICS OF THE STACKED CROSS-SPECTRUM

Based on eq. (6), variance and pseudo-variance of stacked cross-spectrum can be related to the standard errors of its real and imaginary parts,

$$ \mathrm{Var}[R_{ab}^N(\omega)] = E[n_C(\omega) n_C^*(\omega)] = \sigma_R^2(\omega) + \sigma_I^2(\omega) $$

$$ \mathrm{pVar}[R_{ab}^N(\omega)] = E[n_C(\omega) n_C(\omega)] = \sigma_R^2(\omega) - \sigma_I^2(\omega) + i2E[n_R(\omega) n_I(\omega)] \quad (B1) $$
where $\sigma_R^2(\omega)$ and $\sigma_I^2(\omega)$ are variances for the real and imaginary parts of stacked cross-spectrum, respectively. Combining eq. (B1) with eq. (5), we have

$$2E[n_R(\omega)n_I(\omega)] = \mathrm{Im}\left[\frac{1}{N}C_{ab}(\omega)C_{ab}(\omega)\right] \qquad (B2)$$

which means that $n_R(\omega)$ and $n_I(\omega)$ are weakly correlated.

Next the variance of causal part estimator is

$$\begin{aligned}
4\mathrm{Var}[\hat{\alpha}(\omega)] &= E[|\tilde{n}_R(\omega) - i\tilde{n}_I(\omega)|^2] \\
&= E[|n_R(\omega) + iHn_R(\omega) - i[n_I(\omega) + iHn_I(\omega)]|^2] \\
&= E[n_R^2(\omega) + Hn_I^2(\omega) + 2n_R(\omega)Hn_I(\omega) + Hn_R^2(\omega) + n_I^2(\omega) - 2Hn_R(\omega)n_I(\omega)] \\
&= 2E[n_R^2(\omega)] + 2E[n_I^2(\omega)] = 2[\sigma_R^2(\omega) + \sigma_I^2(\omega)]
\end{aligned}
\qquad (B3)$$

where the cross terms have no linear correlations due to the Hilbert transform (see Appendix C): $E[2n_R(\omega)Hn_I(\omega)] = 0$ and $E[2Hn_R(\omega)n_I(\omega)] = 0$. The Hilbert transform also preserves the distribution of the zero-mean random variables $n_R(\omega)$ and $n_I(\omega)$ 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(\omega)$ and $n_I(\omega)$ are only approximately the same as those of $n_R(\omega)$ and $n_I(\omega)$. The variance for anticausal part estimator can be derived in a similar way,

$$\begin{aligned}
4\mathrm{Var}[\hat{\beta}(\omega)] &= E[|\tilde{n}_R(\omega) + i\tilde{n}_I(\omega)|^2] \\
&= E[|n_R(\omega) + iHn_R(\omega) + i[n_I(\omega) + iHn_I(\omega)]|^2] \\
&= E[n_R^2(\omega) + Hn_I^2(\omega) - 2n_R(\omega)Hn_I(\omega) + Hn_R^2(\omega) + n_I^2(\omega) + 2Hn_R(\omega)n_I(\omega)] \\
&= 2E[n_R^2(\omega)] + 2E[n_I^2(\omega)] = 2[\sigma_R^2(\omega) + \sigma_I^2(\omega)].
\end{aligned}
\qquad (B4)$$

The phase angle information can be derived from the analytic form of the real part of stacked cross-spectrum $\tilde{R}_R^N(\omega)$ in eq. (9),

$$\tan\varphi(\omega) = \frac{\mathrm{Im}\tilde{R}_R^N(\omega)}{\mathrm{Re}\tilde{R}_R^N(\omega)} = \frac{\mathrm{Im}\tilde{C}_{ab,R}(\omega) + Hn_R(\omega)}{\mathrm{Re}\tilde{C}_{ab,R}(\omega) + n_R(\omega)}. \qquad (B5)$$

The random variables on the right-hand side are $n_R(\omega)$ and its Hilbert transform $Hn_R(\omega)$, which are not correlated (according to Appendix C).

The phase variance can be estimated with error propagation equation,

$$\sigma_{\psi}^{2}(\omega) = \left(\frac{\partial\varphi}{\partial n_{R}}\right)^{2}\sigma_{n_{R}}^{2}(\omega) + \left(\frac{\partial\varphi}{\partial Hn_{R}}\right)^{2}\sigma_{Hn_{R}}^{2}(\omega). \qquad (B6)$$

Evaluating partial derivatives in eq. (B6) and knowing $n_R(\omega)$ has the same variance with its Hilbert transform $Hn_R(\omega)$, the phase variance becomes,

$$\sigma_{\psi|R}^2(\omega) = \frac{\sigma_{n_R}^2(\omega)}{[\hat{\alpha}(\omega) + \hat{\beta}(\omega)]^2} \qquad (B7)$$

where $\hat{\alpha}(\omega) + \hat{\beta}(\omega) = |\tilde{R}_R^N(\omega)|$ is simply the envelope of the real part of stacked cross-spectrum. Eq. (B7) can be viewed intuitively as inverse square of SNR of the real part stacked cross-spectrum.

Following similar approaches, the phase variances for causal and anticausal parts are, respectively,

$$\begin{aligned}
\hat{\sigma}_{\psi|\alpha}^2(\omega) &= \frac{[\sigma_R^2(\omega) + \sigma_I^2(\omega)]/4}{[\hat{\alpha}(\omega)]^2} \\
\hat{\sigma}_{\psi|\beta}^2(\omega) &= \frac{[\sigma_R^2(\omega) + \sigma_I^2(\omega)]/4}{[\hat{\beta}(\omega)]^2}.
\end{aligned}
\qquad (B8)$$

## APPENDIX C: LINEAR 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,

$$C_{XY}(f_1, f_1 + \Delta f) = D_{XY}(f_1) \delta_D(\Delta f), \qquad (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

$$ 
\begin{aligned}
C_{XY}(f_1, f_1 + \Delta f) &= E[X^*(f_1) Y(f_1 + \Delta f)] \\
&= E\left[\int_{-\infty}^{\infty} x^*(t) \exp(i2\pi f_1 t) dt \int_{-\infty}^{\infty} y(t') \exp[-i2\pi (f_1 + \Delta f) t'] dt'\right] \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} E[x^*(t)y(t')] \exp[i2\pi f_1 (t-t')] \exp(-i2\pi \Delta f t') dt dt' \\
&= \int_{-\infty}^{\infty} \left\{ \int_{-\infty}^{\infty} E[x^*(t'+w)y(t')] \exp[i2\pi f_1 w] dw \right\} \exp(-i2\pi \Delta f t') dt',
\end{aligned}
\quad (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 eqs (C1) and (C2), we have

$$ D_{XY}(f_1) = \int_{-\infty}^{\infty} E[x^*(t' + w)y(t')] \exp[i2\pi f_1 w] dw, \quad (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[-i \operatorname{sgn}(t)y(t)]. \quad (C4) $$

Therefore the covariance function of $X(f)$ and $HY(f)$ is,

$$ 
\begin{aligned}
C_{XHY}(f_1, f_1 + \Delta f) &= E[X^*(f_1) HY(f_1 + \Delta f)] \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} E[x^*(t)(-i) \operatorname{sgn}(t') y(t')] \exp[i2\pi f_1 (t-t')] \exp(-i2\pi \Delta f t') dt dt' \\
&= \int_{-\infty}^{\infty} \left\{ \int_{-\infty}^{\infty} E[x^*(t'+w) y(t')] \exp[i2\pi f_1 w] dw \right\} (-i) \operatorname{sgn}(t') \exp(-i2\pi \Delta f t') dt' \\
&= D_{XY}(f_1) \int_{-\infty}^{\infty} (-i) \operatorname{sgn}(t') \exp(-i2\pi \Delta f t') dt' = 0.
\end{aligned}
\quad (C5) $$

Therefore $X(f)$ and $HY(f)$ are not linearly correlated, which still holds if the noise covariance function $C_{XY}(f_1, f_1 + \Delta f)$ is a real symmetric function of $\Delta f$. As a special case, $X(f)$ and $Y(f)$ can be the same function. Therefore, $X(f)$ and $HY(f)$ are not linearly correlated when condition in eq. (C1) holds.

## APPENDIX D: 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,

$$ R_{n,\text{ft}}(\tau; \omega_0) = \frac{1}{2\pi} \int_{-\infty}^{\infty} S_{n,\text{ft}}(\omega; \omega_0) \exp[i\omega\tau] d\omega \quad (D1) $$

where $S_{n,\text{ft}}(\omega; \omega_0)$ is the power spectral density of the noise narrow-band filtered at centre frequency $\omega_0$.

The real part of the narrow-band filtered noise term on time domain cross-correlation according to eq. (15) is

$$ n_{C,\text{ft}}(t; \omega_0) = \operatorname{Re}[\exp(i\omega_0 t)] A(t) *_C n_C(t; \omega_0). \quad (D2) $$

Applying inverse Fourier transform to eq. (D2), the narrow-band filtered noise spectra are

$$ n_{C,\text{ft}}(\omega; \omega_0) = \frac{1}{2} [\tilde{A}(\omega - \omega_0) + \tilde{A}(\omega + \omega_0)] n_C(\omega_0). \quad (D3) $$
Its power spectral density is

$$S_{n,\text{ft}}(\omega; \omega_0) = \frac{1}{4} [\tilde{A}^2(\omega - \omega_0) + \tilde{A}^2(\omega + \omega_0)] E[|n_C(\omega_0)|^2]. \qquad (D4)$$

Substituting eq. (D4) back into eq. (D1), the correlation function can be evaluated as

$$\begin{aligned} R_{n,\text{ft}}(\tau; \omega_0) &= \frac{1}{4\sqrt{2}} [\exp(i\omega_0\tau) A(\tau) + \exp(-i\omega_0\tau) A(\tau)] E[|n_C(\omega_0)|^2] \\ &= \frac{1}{2\sqrt{2}} \cos(\omega_0\tau) A(\tau) [\sigma_R^2(\omega_0) + \sigma_I^2(\omega_0)]. \end{aligned} \qquad (D5)$$

Assuming the noise on real time domain cross-correlation in eq. (D2) is an ergodic random process, the Mean-Square of it is

$$\mathrm{MS}_{n,\mathrm{TD}}(\omega_0) = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} |n_{C,\mathrm{ft}}(t; \omega_0)|^2 dt = E\left[\left[n_{C,\mathrm{ft}}(t; \omega_0)\right]^2\right] = R_{n,\mathrm{ft}}(\tau = 0; \omega_0). \qquad (D6)$$

Based on eq. (D6), the Root Mean-Square (RMS) of the time domain noise is

$$\mathrm{RMS}_{n,\mathrm{TD}}(\omega_0) = \frac{1}{2} \sqrt{\frac{\sigma_R^2(\omega_0) + \sigma_I^2(\omega_0)}{a\sqrt{2\pi}}}. \qquad (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,

$$\mathrm{SNR}_{\alpha}(\omega_0) = \sqrt{\frac{\frac{1}{a}\sqrt{\frac{2}{\pi}}\alpha(\omega_0)}{\sqrt{\sigma_R^2(\omega_0) + \sigma_I^2(\omega_0)}}}. \qquad (D8)$$

Similarly the SNR of the anticausal part cross-correlation in time domain is

$$\mathrm{SNR}_{\beta}(\omega_0) = \sqrt{\frac{\frac{1}{a}\sqrt{\frac{2}{\pi}}\beta(\omega_0)}{\sqrt{\sigma_R^2(\omega_0) + \sigma_I^2(\omega_0)}}}. \qquad (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).