Electronic Supplement to
A Waveform Detector That Targets Template-Decorrelated Signals and Achieves Its Predicted Performance, Part I: Demonstration with IMS Data

by Joshua D. Carmichael

This electronic supplement develops several theoretical and practical aspects of implementing the correlation and cone detector on geophysical waveform data. The Effective Degrees of Freedom section documents an estimator for the number of statistically independent samples N_E in a data stream and describes how it shapes the correlation statistic’s probability density function (PDF). In The Maximum-Likelihood Correlation Statistic section, we derive this PDF, which quantifies our correlation and cone detector’s ability to identify their respective target waveforms. The Cone Detector Statistic section then presents an alternative to the correlation detector hypothesis test and introduces additional, deterministic uncertainty into the target signal model. Therein, we derive the cone detector PDF, illustrate its geometric properties, and model its capability to detect template-decorrelated target waveforms. We support this derivation with Monte Carlo simulations (Fig. S1) and illustrate how deterministic template-target waveform decorrelation inflates thresholds relative to the correlation case (Fig. S2). The Relative Magnitude Parameterization of the Cone Detector Statistic section then relates the probability of detecting a template-decorrelated target waveform to the magnitude of the seismic source that produces this target waveform. Finally, we provide an analysis of three ostensible sources of error from our data processing implementation in the Error Analyses section and describe their influence on our threshold and capability analyses. We find that only data-variance pooling is likely to appreciably affect detector performance.

Effective Degrees of Freedom

This section outlines a method to estimate the effective degrees of freedom N_E of a data stream x as ${\widehat{N}}_E$. This estimate provides an ostensible alternative to a full covariance matrix $\mathbf{\Sigma }\ne {\sigma }^2\mathbf{I}$ for x that is generally required for temporally or spatially correlated data. Density functions for detection statistics that process x are then easily parameterized by the effective number of independent data-stream samples N_E. This scalar theoretically equates to twice the time (T) bandwidth (B) product (2TB) of the data stream over the temporal duration of the correlation template. Real data often show that $N_E\ll 2TB$, however. This occurs both naturally and through processing operations like band-pass filtering, which replace each sample with itself plus a weighted sum of its neighbors and thereby introduce intersample statistical dependence. To quantify the influence such correlation exerts on the shape of our detector’s PDF curves, we implement an empirical estimator for N_E, denoted ${\widehat{N}}_E$, to continuously update such PDF parameterizations (e.g., the correlation statistic PDF $f_R(r;{\mathcal{H}}_k)$). This estimator computes the sample correlation between the multichannel template waveform w and several hundred pseudorandom, commensurate data vectors drawn from nonintersecting segments of preprocessed, signal sparse data within x (see Weichecki-Vergara et al., 2001, their section 2.4). We compute the sample variance ${\widehat{\sigma }}_R^2$ of the resultant correlation time series using 99.9% of the data by excluding 0.01% of the extreme left and right tails of its histogram. This provides the needed statistic to estimate N_E as

We use ${\widehat{N}}_E$ to parameterize $f_R(r;{\mathcal{H}}_k)$ (equation S7), compute detector thresholds $\widehat{\eta }$, and quantify detector performance.

The Maximum-Likelihood Correlation Statistic

This section develops the PDF and detection performance of the sample correlation detection statistic r(x) (equation 3 of the main article). Because the form of the sample correlation PDF differs from that reported elsewhere, we derive the general form here.

The relative square error in approximating a data stream $\mathbf{x}=\mathbf{u}+\mathbf{n}$ (the ${\mathcal{H}}_1$ model of equation 15 of the main article), with a waveform template w, is the ratio of the least-squares error ${\Vert \mathbf{e}\Vert }^2$ to measured signal energy ${\Vert \mathbf{x}\Vert }^2$. This error may be rewritten as a ratio of quadratic forms

in which $\widehat{A}$ is the maximum-likelihood estimate (MLE) for template waveform amplitude. We now rewrite the right side of the second equality in equation (S2) as a ratio of subspace projections

in which W is the subspace span (w), $W^{\perp }$ is the orthogonal complement to W, P_W is the projector onto W, and $P_W^{\perp }$ is the projector onto $W^{\perp }$. The denominator follows from the Pythagorean identity for Hilbert spaces. We define two noncentrality parameters from these terms

in which the expected value and linear-projection operators commute. We combine the previous three equations to rewrite r²(x)

in which $\stackrel{d}{=}$ is distributional equality, ${\chi }_{N_E-1}^2({\lambda }^{\perp })$ is the noncentral chi-square distribution with N_E  −1 degrees of freedom and noncentrality parameter ${\lambda }^{\perp }$, and ${\chi }_1^2(\lambda )$ is the noncentral chi-square distribution with one degree of freedom and noncentrality parameter $\lambda$. From the definition of the beta distribution

(Kay, 1998), in which $\mathrm{B}(t,N_1,N_2,\alpha ,\beta )$ is the doubly noncentral beta distribution function. It is evaluated at $t$ (with the same domain as r²), has N₁ and N₂ degrees of freedom and noncentrality parameters $\alpha$ and $\beta$. The presence or absence of a target signal is indexed by the hypothesis ${\mathcal{H}}_k$ on the data. Hypothesis ${\mathcal{H}}_0$ symbolizes that the data consist of only noise, whereas ${\mathcal{H}}_1$ signifies that the data consist of a signal plus noise. The scalar N_E denotes the effective number of independent samples within x. When the data stream contains only noise, the hypothesis ${\mathcal{H}}_0$ is satisfied, and r² has a central beta distribution, in which ${\lambda }^{\perp }=\lambda =0$. In the presence of signal, the data stream x generally has nonzero projections $P_W(\mathbf{x})$ and $P_W^{\perp }(\mathbf{x})$ that are respectively onto and orthogonal to the noise-contaminated template data vector w. In this case, $\lambda$, ${\lambda }^{\perp }\ne 0$, and r² has doubly noncentral beta distribution. If the target signal is an amplitude-scaled copy of the template waveform, then $\mathbf{x}=A\mathbf{w}+\mathbf{n}$, ${\lambda }^{\perp }=0$, and r² has a noncentral beta distribution. The singly noncentral beta distribution therefore provides an absolute upper bound on the detection performance of a correlation detector consistent with assumptions underlying ${\mathcal{H}}_1$.

We derive the PDF for $r$ from the density of r² using a variable transformation; we additionally consider values r < 0, so that

The form of this distribution function differs from that derived in similar applications (Weichecki-Vergara et al., 2001). In the alternative case, the signal-present data stream was assumed to correlate sample-by-sample with the template waveform, and the test statistic had a Pearson moment-product distribution.

The Cone Detector Statistic

This section develops the cone detection statistic (equation 17 of the main article), its PDF, and illustrates its performance; it was originally introduced in the context of icequake detection (Carmichael, 2013). We first reference the two competing hypotheses of equation (15) of the main article

We further assume that preprocessing operations (like filtering) and naturally occuring temporal data correlation reduces the number of statistically independent samples in x to $N_E<N$. Given this parameterization, the PDF under ${\mathcal{H}}_0$ is denoted as $f_0(\mathbf{x};{\mathcal{H}}_0)$ and the PDF under ${\mathcal{H}}_1$ as $f_1(\mathbf{x};{\mathcal{H}}_1)$, in which

We estimate the unknown parameters under each model listed in equation (15) of the main article and select the applicable hypothesis using a log-generalized likelihood ratio test (log-GLRT). This test evaluates the PDFs for the competing hypotheses at the MLEs of their unknown parameter values and then compares their log ratio to a threshold value $\eta$. The test’s decision rule is a conditional, scalar inequality

The MLE of the variance under each hypothesis is

(Kay, 1993; Scharf and Friedlander, 1994), in which the subscripts on each sample variance estimate in equation (S11) indicate the applicable hypothesis, and N is the number of samples in x, not to be confused with N_E. To estimate the distribution mean under ${\mathcal{H}}_1$, we evaluate $f_1(\mathbf{x};{\mathcal{H}}_1)$ at the MLE for ${\sigma }_1$ and perform a constrained maximization of $\mathbf{u}\in \mathcal{C}$

The solution to this equation

is the MLE for u. That is, $\widehat{\mathbf{u}}$ is the vector that minimizes the distance between the observed data x and all points that constitute the multiplet set $\mathcal{C}$. This defines the projection of x onto $\mathcal{C}$ as the unique signal that is either interior to, or on the boundary $\partial \mathcal{C}$ of $\mathcal{C}$ (Stark and Yang, 1998). The sample variance and cone-element MLEs from equation (S11) reduce $\mathrm{\Lambda }(\mathbf{x})$ to

We obtained equation (S14) without specifying $\mathcal{C}$ aside from its convexity. This result therefore applies to any normally distributed data confined in a convex set. In the case that $\mathcal{C}$ is described by the correlation constraint of equation (14) of the main article, the projection energy has a conic decomposition

as given by Moreau’s decomposition theorem (Moreau, 1962), which is analogous to the orthogonal subspace decomposition from linear analysis (Stark and Yang, 1998). The log ratio is then expressible as

The right side of equation (S14) is of the form $-\ln (1-x^2)$, which is monotonically increasing for $0\le x\le 1$. Because ${\Vert P_{\mathcal{C}}(\mathbf{x})\Vert }^2/{\Vert \mathbf{x}\Vert }^2\le 1$, equation (S16) may be inverted for its argument to provide an equivalent test statistic s(x) for the decision rule introduced in equation (S10)

Equation (S17) demonstrates that s(x) compares the ratio of the projected signal energy to the original signal energy. The exact projection $P_{\mathcal{C}}(\mathbf{x})$ of x, from a general Hilbert space and onto $\mathcal{C}$, is derived in Stark and Yang (1998, pp. 111–113). Vector $P_{\mathcal{C}}(\mathbf{x})$ is a nonlinear projection that depends on the value of x and is either in $\mathcal{C}$, on its boundary $\partial \mathcal{C}$, or zero. We document an equivalent form of that projection here

The constants and vectors in equation (S18) are

The vector $\widehat{\mathbf{w}}$ in equation (S19) is the normalized multichannel template waveform defined by equation (14) of the main article, not to be confused with a parameter estimate. Using the definitions in equation (S19), the test statistic s(x) of equation (16) of the main article is then expressible as

The detection statistic s(x) must be equivalent to the correlation coefficient as ${\rho }_{\infty }$ approaches a limiting value of 1. To demonstrate this, we first note that any projected signal has a decreasing probability of lying inside the cone as ${\rho }_{\infty }$ decreases (equation S20). Similarly, any projected signal has a decreasing probability zero of lying on the cone vertex. It follows that only the projection onto the limiting cone boundary $\partial \mathcal{C}$ is nontrivial, so that

in which $\underset{{\rho }_{\infty }\to 1^{-}}{\lim }$ is the limit of ${\rho }_{\infty }$ approaching from values less than 1, whereby $c\to 0$.

We derive the requisite PDF for s(x) from its cumulative distribution function (CDF) $F_S(s;{\mathcal{H}}_k)$ using the law of total probability. In words, this law states that the probability that the cone statistic random variable S takes on a value as large as s(x) is the sum of three conditional probabilities: the probability that (1) $S<s(\mathbf{x})$, given $\mathrm{Pr}\{P_{\mathcal{C}}(\mathbf{x})=\mathbf{0}\}$, plus (2) the probability $S<s(\mathbf{x})$, given $\mathrm{Pr}\{P_{\mathcal{C}}(\mathbf{x})\in \partial \mathcal{C}\}$, plus (3) the probability that $S<s(\mathbf{x})$, given $\mathrm{Pr}\{P_{\mathcal{C}}(\mathbf{x})\in \mathcal{C}°\}$

in which $\mathrm{Pr}\{A|B\}$ is the conditional probability of A, given B is true and statement $P_{\mathcal{C}}(\mathbf{x})\in \mathcal{C}°$ is equivalent to $P_{\mathcal{C}}(\mathbf{x})=(\mathbf{x}\in \mathcal{C})$. To express $F_S(s;{\mathcal{H}}_k)$ in a computationally evaluable form, we write its three conditioning factors in terms of $r/\sqrt{1-r^2}$ through equation (S20). Two of these terms are trivial to evaluate from the definition for $s(\mathbf{x})=\Vert P_{\mathcal{C}}(\mathbf{x})\Vert /\Vert \mathbf{x}\Vert$

The other terms in equation (S22) are expressible using $r/\sqrt{1-r^2}$ through equation (S20)

We evaluate these probabilities by developing the CDF $F_Q(q;{\mathcal{H}}_k)$ for the ratio $q=r/\sqrt{1-r^2}$ ($k=\mathrm{0,1}$). To do so, we make a change of variables with $F_R(r;{\mathcal{H}}_k)$

in which $F_R(r;{\mathcal{H}}_k)$ is the CDF for correlation statistic $r$. We then use equation (S25) to write the identities from equation (S24) in terms of $F_Q(q;{\mathcal{H}}_k)$:

Last, we evaluate the derivative of $\mathrm{Pr}\{P_{\mathcal{C}}(\mathbf{x})\in \partial \mathcal{C}\}$ in equation (S22) through a variable change on $r$, in which $-c<r/\sqrt{1-r^2}<1/c$. To do so, we first note that $s(\mathbf{x})={\rho }_{\infty }(r-c\sqrt{1-r^2})$ is invertible over this domain and express r as a function of s

The PDF of s(x) over $-c<r/\sqrt{1-r^2}<1/c$ is therefore

in which $|P_{\mathcal{C}}(\mathbf{x})\in \partial \mathcal{C}$ indicates the restricted domain of $r$. To obtain the PDF for s(x) over $-1\le r\le 1$, we differentiate equation (S22) and substitute equations (S26) and (S28) into the results. This produces a density function that depends on only s, c, and $f_R(·;{\mathcal{H}}_k)$

because the top line of equation (S22) is constant and differentiates to zero.

We assessed the validity of our derivation using a Monte Carlo simulation, whereby we projected random noise and noise-contaminated signal vectors onto several convex cones of increasing aperture and computed the statistic s(x). This simulation demonstrates that equation (S29) agrees with our empirical histograms (Fig. S1).

We additionally compared our cone detector thresholds against constant false-alarm-on-noise constraints. We thereby inverted for cone detector thresholds using a fixed value for the effective degrees of freedom ($N_E=400$) over a grid of decorrelation parameters ($1-{\rho }_{\infty }$) that ranged from 0 to 0.25. We repeated this process for several false alarm rates (Fig. S2). The resulting detection thresholds increase most rapidly for small changes near ${\rho }_{\infty }=1$, at which the signal space of the cone geometrically collapses to the 1D subspace used by the associated correlation detector. The slope of the curves here may become undefined. We will explore this result more quantitatively in future work.

Relative Magnitude Parameterization of the Cone Detector Statistic

An unbiased estimator for the relative magnitude $\mathrm{\Delta }m$ between an explosion that generates a noisy waveform $\mathbf{x}=\mathbf{u}+\mathbf{n}$ and an explosion that generates a detector template w is

(Carmichael and Hartse, 2016, their equation B.8), in which ${\sigma }_1^2$ denotes an estimate for noise variance in x (subscript denotes hypothesis ${\mathcal{H}}_1$ is true) and N is the number of samples in x, not to be confused with N_E. To relate ${\Vert \mathbf{x}\Vert }^2$ to the noncentrality parameters $\lambda$ and ${\lambda }^{\perp }$ that shape the PDF $f_S(s;{\mathcal{H}}_0)$, we again use the standard Pythagorean identity for Hilbert spaces

in which the subspace projection terms are defined in equation (S3). Next, we use the definitions of u, $\lambda$, and ${\lambda }^{\perp }$ to rewrite ${\Vert \mathbf{x}\Vert }^2$ as

in which “noise term” is an inner product expression that includes n, and ${\sigma }_1^2$ is the true variance of the target data. The noncentrality parameters in equation (S32) are related by Carmichael and Hartse (2016, their equation A.11)

and the expected value of “noise term” in equation (S32) is

in which $\mathbb{E}\{⟨\mathbf{u},\mathbf{n}⟩\}=0$, because the noise is zero mean. Finally, we combine the preceding equations (this section), remove the noise-bias term ${\sigma }_1^{2}N$, and write $\lambda$ in terms of relative magnitude

Equation (S35) thereby expresses the noncentrality parameter for both the correlation and cone detection statistic in terms of the four fundamental scalars describing the wavefield: the deterministic correlation ${\rho }_{\infty }$ between the template and target waveforms, the signal energy ${\Vert \mathbf{w}\Vert }^2$ of the template waveform, the noise variance ${\sigma }_1^2$ of the target data, and the relative magnitude $\mathrm{\Delta }m$ between the template and target source.

For our purposes, equation (S35) parameterizes $\lambda$ by relative magnitude. In such cases, the term ${\rho }_{\infty }$ is estimable from equation (12) of the main article, ${\sigma }_1^2$ is estimable from the top term of equation (S11), and ${\Vert \mathbf{w}\Vert }^2$ is effectively deterministic because correlation detectors generally implement high-signal-to-noise ratio (SNR) templates.

Error Analyses

We identified three potentially significant sources (risks) of error in our detection capability assessments and empirical comparison. The first risk of error is attributable to certain details of template selection. Specifically, waveforms recorded on International Monitoring System (IMS) arrays with large differences in source–receiver separation distance show temporal gaps in start time of the high-SNR portions of the recorded signals and therefore require sample imputation. Zero padding such data to equalize length can lead to several biases, however. For example, supplementing data with a large number of zeros causes the empirical null PDF (histogram) for the correlation detection statistic to become more concentrated around its mean and thereby artificially lowers the detector’s threshold. The statistic histogram then fits the predicted PDF poorly and biases estimates of the degree of freedom parameter ${\widehat{N}}_E$, because the correlation variance is reduced by the presence of zeros (equation S1). To mitigate these problems and facilitate template scanning, we therefore abstained from zero-padding the intermediate, noisy portion of our waveform template. Instead, we shifted each seismogram channel by an amount equal to its start time, minus the earliest start time among all template channels (template and target data). We thereby maintained the signal-only length of our template, avoided padding, and kept all data aligned to the same origin time. In addition to preventing estimation bias, this shifting mitigated needless computation of padded data and prevented divide-by-zero errors. We therefore consider our template selection to present a low (direct) risk of error.

Estimation of ${\rho }_{\infty }$ presents a second potential source of error. This arises largely from ambiguous estimation schemes for waveform SNR, which influence the variability of ${\rho }_0$ (equation 11 of the main article), which, in turn, inversely scales ${\widehat{\rho }}_{\infty }$ (equation 12 of the main article). We mitigated ambiguity problems by carefully selecting an associated low-variance estimate for SNR, which normalizes ${\rho }_{\infty }$. One “common-sense” estimate for SNR is the ratio of an N-point sample variance estimate of the data proceeding the detected signal, divided by an N-point sample variance estimate of data preceding the detected signal (a renormalized short-term average/long-term average [STA/LTA]). It is obvious, however, that any such estimate will be biased by background signals contaminating the data steam, which reduce the resultant SNR estimate for u. It follows that ${\widehat{\rho }}_{\infty }$ will be a biased estimator and give lower-than-true deterministic correlation values. A better approach requires preprocessing target data with a power (STA/LTA) detector, removing samples that exceed a threshold for signal declaration, and then computing the noise variance ${\widehat{\sigma }}_0^2$ from this remaining data. The quotient

then gives a reduced-bias SNR estimate. We used this estimator to compute ${\rho }_0$ as:

We therefore consider our estimates of ${\rho }_0$ to be robust and unlikely to induce significant uncertainty.

Last, noise variance estimation presents an additional risk of error. As we noted before (see the Estimating Deterministic Decorrelation section in the main article), background seismicity adds signal to the time series and can thereby bias such estimates. We therefore processed our target data in parallel with an STA/LTA detector that operated at a ${10}^{-6}$ false-alarm-on-noise probability, removed data that exceeded the associated threshold, and used the remaining (almost) signal-free data to estimate noise variance on each IMS channel. This scheme thereby avoided bias from signal. A second form of bias was also present, however. This latter bias originated from the nonuniform number of channel samples processed by the detector. Specifically, we noted above that our template included large temporal gaps over portions of the waveform, and consequently, records from MJAR contribute less to each detection statistic than do longer duration records at MKAR. An alternative noise-variance estimate that accounts for such disparate channel contribution employs the pooled variance ${\widehat{\sigma }}_P^2$, given by

in which ${\widehat{\sigma }}_k^2$ is the noise-variance estimate from channel $k$. We performed a subsequent analysis using this estimator and found that it was often smaller than our naive estimation that used equally weighted data. This would have increased the effective size of the noncentrality parameter $\lambda$ that controls predicted detection power. It may explain why, at times, our observed detection capability exceeded the concurrent predicted detection capability. We therefore suggest that performance discrepancy may generally result from an inconsistency in noise variance estimation.

---

Figures

Figure S1. Null hypothesis PDFs for three cases of deterministic template-target waveform correlation uncertainty: ${\rho }_{\infty }=1$, ${\rho }_{\infty }=0.99$, ${\rho }_{\infty }=0.90$ (in which N_E = 400). The PDFs for r(x) and s(x) equate for identically shaped waveforms as shown by the purple and black curves. The shaded region shows a Monte Carlo simulation of the PDF for s(x) when ${\rho }_{\infty }=0.9$

Figure S2. Cone detector thresholds $\eta$ (equation 21 of the main article) for constant values of ${\mathrm{Pr}}_{\mathrm{FA}}$ compared against deterministic uncertainty in the template-target waveform cross correlation $1-{\rho }_{\infty }$, in which $N_E=400$. The leftmost tangent line to the ${10}^{-6}$ curve shows a rapid increase in $\eta$ for small uncertainties in deterministic template-to-target waveform correlation relative to near-linear increases in $\eta$ for larger uncertainties. ${\rho }_{\infty }=0$ corresponds to the correlation detector.

---

[ Back ]