Hypothesis Tests on Rayleigh Wave Radiation Pattern Shapes: A Theoretical Assessment of Idealized Source Screening

The most general symmetric moment tensor for a single point source is a superposition of an explosion with isotropic moment $M_{I}$ that effectively colocates with a non-VDS faulting source of moment $M_{0}$ and a damage source, which is equivalent to a compensated linear vector dipole (CLVD) with moment $M_{\text{CLVD}}$. Moment $M_{\text{CLVD}}$ is positive if explosions create dilatation along the vertical axis of damage, and negative if explosions drive contraction along that axis. We assume that the fault is parameterized by strike ($\phi_{s}$), rake ($\lambda$) and dip ($\delta$) angles. The radiation pattern coefficients from eq. 2 are then (algebra omitted):

where we have used a moment decomposition (Patton and Taylor, , 2011, eq. 15). In eq. 4, scalar $SS$ quantifies the strength of the strike-slip faulting component, whereas scalar $DS$ quantifies the strength of the dip-slip component. These terms are (Ekström and Richards, , 1994, eq. 21 and eq. 22):

where the equivalent expression for $SS$ elsewhere (Patton and Taylor, , 2011, eq. 16) is missing a factor of two in the dip angle. We express the full Rayleigh wave displacement by combining eq. 2, eq. 4, and eq. 5. The superposition of the explosive, damage and non-VDS faulting sources thereby creates a Rayleigh wave displacement field that is:

where $\psi$ $=$ $\phi$ $-$ $\phi_{s}$ is the difference between the azimuthal and strike angles; because we use $\psi$ hereon, we refer to it as the “azimuthal” angle, despite abuse in terminology. Regardless, the function that scales $\boldsymbol{g}$ is the radiation pattern $\mathcal{R}$, and includes the same functional form as the scalar factor that appears in eq. 2:

Eq. 4 with eq. 7 demonstrates that $\bar{\mathcal{R}}$ is azimuthally constant, and that azimuthal variability arises from non-VDS faults where $\delta$ $\neq$ $\pm\pi$/2, $\lambda$ $\neq$ $\pm\pi$/2. Eq. 7 is also the starting point for estimating the optimal receiver sampling of the Rayleigh wave radiation pattern of shallow sources (Section 4.3).

We first assume a network of $N\geq 8$ sensors records a Rayleigh wave radiation pattern with known process variance ($\sigma_{\mathcal{R}}^{2}$ in eq. 16), unknown faulting parameters ($\boldsymbol{\theta}$ in eq. 11), and unknown stochastic variance ($\sigma_{M}^{2}$ in eq. 16). This case idealizes explosion monitoring scenarios in which an underground testing complex colocates with a shallow fault of unknown geometry. This case also includes the scenario that all variability is stochastic in origin ($\sigma_{\mathcal{R}}^{2}$ $=$ $0$) as an example. Under this general model, a sensor network then records a low magnitude seismic event that originates near a test site, and an observer must estimate the probability that this event is explosive or tectonically sourced by a non-VDS fault (hypothesis $\mathcal{H}_{0}$ or $\mathcal{H}_{1}$).

To determine which hypothesis Rayleigh wave records likely describe, we use a generalized log-likelihood ratio test (log-GLRT). This test compares the logarithmic ratio of the alternative hypothesis PDF to the null hypothesis PDF, where each density is evaluated at maximum likelihood estimates (MLE) of its unknown parameters to form a test statistic. The log-GLRT is equivalent to a conditional decision rule on this scalar test statistic $L_{\boldsymbol{\theta}}^{(I)}\left(\boldsymbol{\mathcal{R}}\right)$:

We compute the MLEs of the parameters in Appendix B using projector matrices and reduce the scalar statistic into a ratio of two scaled, statistically independent terms. This provides a test statistic $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ equivalent to eq. 17 for screening circular from non-circular radiation patterns:

The rank-1 projector matrix $\boldsymbol{P}_{\boldsymbol{X}}$ in the numerator of eq. 18 is:

The rank-3 projector matrix $\boldsymbol{P}_{\boldsymbol{H}}$ $=$ $\boldsymbol{H}\left(\boldsymbol{H}^{\text{T}}\boldsymbol{H}\right)^{-1}\boldsymbol{H}^{\text{T}}$ projects vectors onto the subspace spanned by $\boldsymbol{H}$, $\text{span}\{\boldsymbol{H}\}$. Projector matrix $\boldsymbol{P}_{\boldsymbol{H}}^{\perp}$ $=$ $\boldsymbol{I}$ - $\boldsymbol{P}_{\boldsymbol{H}}$ in the denominator of eq. 18 then projects vectors onto the space orthogonal to $\text{span}\{\boldsymbol{H}\}$. We thereby interpret the screening statistic $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ in eq. 18 as a scaled ratio of the radiation pattern energy due to non-VDS faulting (note structure of $\boldsymbol{A}$), divided by the residual radiation pattern energy not attributable to the Rayleigh wave model (the system $\boldsymbol{H}$). Probability theory predicts that this quotient has a noncentral-$F$ distribution with one, and $N-3$ degrees of freedom. The density for $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ is then shaped by a scalar noncentrality parameter $\Lambda$ if $\mathcal{R}_{k}$ includes Gaussian noise ($1\leq k\leq N$) and the radiation pattern is non-circular. We write its distributional dependence as $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ $\sim$ $\mathcal{F}_{1,N-3}\left(L,\Lambda>0\right)$ and the cumulative, noncentral $\mathcal{F}$ distribution for $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ as:

These distributions have a well-defined variance when their second degree of freedom exceeds four, so that the screening statistic $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ is only well-defined when $N\geq 8$. This sampling implies that at least two sensors sample each quadrant of the Rayleigh radiation pattern (given full azimuthal coverage by sensors), which mitigates spatial aliasing of four-lobed radiation patterns. This statistic is identical in form to the subspace detector (Scharf and Friedlander, , 1994) that is often applied in seismic monitoring, and has completely analogous interpretations for waveform detection. Consequently, many standard detection theory texts document the properties and analytical form for $F_{1,N-3}\left(L,\Lambda\right)$, and computational packages like MATLAB include routines to calculate random variables and parameters of the noncentral-$F$ distribution.

Crucially, this theory shows that the scalar $\Lambda$ in eq. 20 completely quantifies the screening capability of $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ to test between $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ at fixed $N$; eq. B.5 shows its general analytical form. Here, we use those results to write $\Lambda$ as the product of one factor that depends on deployment azimuth $\psi$ and a second factor that depends on source and noise parameters:

This factorization shows that $\Lambda$ is the product of a scalar term $\left(\boldsymbol{A}\left[\boldsymbol{H}^{\text{T}}\boldsymbol{H}\right]^{-1}\boldsymbol{A}^{\text{T}}\right)^{-1}$ that stores the $N$ sensor deployment locations, and a term $\left(DS+SS\right)^{2}/(\sigma_{M}^{2}+\sigma_{\mathcal{R}}^{2})$ that quantifies the non-VDS faulting signal energy, relative to noise and random process variance. It is analogous to the signal-to-noise plus interference ratio (SNIR) of waveform processing. Explosive or VDS sources that trigger no Rayleigh waves include parameters $\sigma_{\mathcal{R}}^{2}$ $=$ $DS$ $=$ $SS$ $=$ $0$, and $\Lambda$ $=$ $0$ (as expected). Strike-slip faults maximize $\Lambda$ among other faulting sources, given equal sensor placement, because $\left(DS+SS\right)^{2}/\left(\sigma_{M}^{2}+\sigma_{\mathcal{R}}^{2}\right)$ $\leq$ $M_{0}^{2}/\left(\sigma_{M}^{2}+\sigma_{\mathcal{R}}^{2}\right)$ (see eq. 5). Importantly, eq. 18 and eq. 21 demonstrate that discriminating between a circular and non-circular Rayleigh wave radiation pattern under the Case I assumptions is entirely equivalent to a signal detection problem. We further note the random process variance $\sigma_{\mathcal{R}}^{2}$ acts to simply inflate the effective noise variance and reduce any observed, non-VDS faulting signal SNR under $\mathcal{H}_{1}$. This asymmetry in variance between the null and alternative hypotheses means that a high SNR faulting signal with high random signal variance is indistinguishable from a lower SNR faulting signal with zero random signal variance. We will assume $\sigma_{\mathcal{R}}^{2}=0$ hereon, so that $\Lambda$ and $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ are each simplified from their more general forms. This simplification preserves our physical insight into the radiation pattern testing problem without adding cumbersome notation. The faulting SNIR then reduces to just SNR.

Under these simplifying (but still rather general) assumptions, the probability $\text{Pr}_{D}$ of detecting any general faulting signal, and then deciding that $\boldsymbol{\mathcal{R}}$ records a non-circular radiation pattern is the probability $\text{Pr}_{D}$ of measuring a value of $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ that exceeds a threshold $\eta$ (eq. 18):

We select a threshold $\eta$ for deciding a radiation pattern is non-circular using the Neyman-Pearson criteria. This criteria uses a constant false-attribution probability $\text{Pr}_{FA}$ to invert for $\eta$:

where we set $\text{Pr}_{FA}$ to an acceptable value (e.g., $10^{-3}$) and $F_{1,N-3}^{-1}\left(\bullet\right)$ is the inverse of the central $F$ distribution with one and $N-3$ degrees of freedom. The screening probability $\text{Pr}_{D}$ in eq. 22 monotonically increases with $\Lambda$ at fixed $N$ and $\eta$. The locus of points $\left(\Lambda,\text{Pr}_{D}\right)$ thereby defines performance curves that measure the screening capability of the decision rule in eq. 18 versus SNR of the faulting signal at fixed $\text{Pr}_{FA}$ (Fig. 3). Deployment constraints and tectonic settings can each limit the permissible parameter space of $\Lambda$ to a subset $\mathcal{S}\left(\Lambda\right)$. These limits include the azimuthal deployment range of sensor numbers and locations that quantify the leftmost factor in eq. 21, and variability in the source focal mechanism, which quantifies the rightmost factor in eq. 21. In such cases, we estimate the expected screening capability $\bar{\text{Pr}}_{D}$ of eq. 18 as the average of the detection rate $\text{Pr}_{D}$ over such parametric constraints. We use eq. 22 to write this average as:

in which $\lvert\mathcal{S}\left(\Lambda\right)\rvert$ is the size of the set (e.g., the number of realizations in a Monte Carlo simulation). The sum in eq. 24 equates to an integration when the subset of permissible values of $\Lambda$ is a continuum (e.g., $\lvert\psi\rvert$ $\leq$ $3\pi/4$) and gives a marginal distribution for $\bar{\text{Pr}}_{D}$ (e.g., $\bar{\text{Pr}}_{D}$ $=$ $4/\left(3\pi\right)\int_{-3\pi/4}^{+3\pi/4}d\psi\left(1-F_{1,N-3}\left(\eta,\Lambda\right)\right)$.

Fig. 4, Fig. 5 and Fig. 6 each show qualitatively similar, mean screening curves (eq. 24) that each index distinct deployment and source parameters. Fig. 4 shows the probability that eq. 24 will screen strike-slip earthquakes from explosions with $N$ $=$ $12$ sensors deployed over three distinct azimuthal ranges against faulting signal (pictured at right). Each of these three curves show the dependence of screening probability on the faulting signal SNR $\left(DS+SS\right)^{2}/\sigma_{M}^{2}$, averaged over 100 individual screening curves. Scalar $\Lambda$ uniformly samples receiver deployments confined to these azimuthal ranges where source parameters are $\delta$ $=$ $\frac{\pi}{2}$, $\lambda$ $=$ $0$, $(DS+SS)^{2}$ $=$ $M_{0}$, and faulting SNR $=$ $M_{0}^{2}/\sigma_{M}^{2}$. The three curves demonstrate that screening probability markedly decreases as azimuthal gap increases with a fixed number of stations, even as station density over deployment azimuth (sensor number per arc length) increases. Shading about the top blue curve associated with full azimuthal sampling marks the one-standard deviation ($\pm\sigma$) variability.

Fig. 5 illustrates that earthquake signals triggered by other focal mechanisms are more difficult to screen from explosions relative to strike-slip events ($\delta$ $=$ $\pm\pi/2$, $\lambda$ $=$ $0,\pi$). Specifically, eq. 18 illustrates predicted curves that quantify the capability of statistic $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ to screen distinct earthquakes from explosions using the same deployment strategy with $N$ $=$ $12$ sensors. In this case, each curve is an average of $100$ screening curves that randomly sample a circular instrument deployment around three different sources (focal mechanisms at right). These curves demonstrate that the radiation pattern test screens non-strike-slip earthquakes from explosions less effectively than strike-slip earthquakes, even with identical station coverage. Seismic events with a Rayleigh wave radiation pattern like that of the 2003 Lop Nor tectonic event that is sampled by an unrestricted deployment of $N$ $=$ $12$ sensors is particularly difficult to screen from an explosion. Shading about the curve that associates with full azimuthal sampling marks the $\pm\sigma$ variability in screening sources like Event 1 from the Rock Valley earthquake.

An increase in sensor deployments ($N$ $>$ $12$) correspondingly increases screening probability between non-VDS faulting and explosion sources. Fig. 6 shows three distinct averages of $100$ screening curves that each associate with three distinct sensor deployments ($N$ $=$ $8,16,24$). Each deployment randomly samples a circle centered about a shallow source like that of the 2003 Lop Nor tectonic event. In this case, a screening statistic that uses $N$ $=$ $24$ sensors can screen non-VDS faulting from explosion sources at roughly the same mean rate as a statistic that uses $12$ receivers to screen sources like Event 1 from the Rock Valley earthquake sequence (Fig. 5). While a $24$ sensor deployment screens shallow sources like that show in Fig. 6 with relatively low variance, the test statistic obviously requires twice the number of observations.

Our interpretation of screening curves that compare probabilities like $\bar{\text{Pr}}_{D}$ against faulting signal SNIR require comment (see Fig. 4, Fig. 5 and Fig. 6). The binary decision rules necessarily decide between a continuum of source types that are not identically explosions or faults. For example, the decision rule in eq. 18 will screen data that records Rayleigh motion sourced by an explosion that superimposes with non-zero contributions from a non-VDS faulting signal with a probability that much greater than the false-attribution probability $\text{Pr}_{FA}=10^{-3}$. A proper false attribution probability constraint that considers an acceptable amount of faulting signal with an explosion would more appropriately marginalize out some faulting effects with a prior on $\Lambda$ from this continuum (Berger and Delampady, , 1987). We address such source-type errors and ambiguities in Section 5.4 of our Discussion. The binary hypothesis testing theory that we use here is standard, however, and we proceed with the caveat that neither the $\mathcal{H}_{0}$ nor $\mathcal{H}_{1}$ in our models fully capture the continuum of source type hypotheses.

We next explore if more optimal sensor placement can improve screening rates of target sources that are particularly challenging to screen from cylindrically symmetric sources.

We suppose an existing $N\geq 8$ receiver network with limited azimuthal coverage requires an auxiliary receiver to maximize its capability to screen shallow non-VDS sources from explosions, through application of $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ as a discriminant. To develop an objective function for this discriminant, we add a row $\left[1,\cos(2\psi)-c,\sin(2\psi)\right]$ to $\boldsymbol{H}$ and make a new matrix $\boldsymbol{H}_{+1}$. We then maximize $\text{Pr}_{D}$ over additional sensor deployment azimuths $\psi_{+1}$ that we estimate as $\hat{\psi}_{+1}$. Such an optimal auxiliary sensor location estimate supplements the existing network to form an updated $N+1$ network that maximizes the probability of detecting a non-VDS signal in the radiation pattern data $\boldsymbol{\mathcal{R}}$. Specifically, this new observation point updates the discriminant in eq. 18, while eq. 23 updates the threshold $\eta$ to maintain the same false attribution probability $\text{Pr}_{FA}$. The optimal solution $\hat{\psi}_{+1}$ depends on both $\eta$ and $\Lambda$, and maximizes $\text{Pr}_{D}$ where the total differential $d\text{Pr}_{D}$ peaks:

The decision threshold decreases by an increment $\Delta\eta$:

whereas the faulting signal changes by increment $\Delta\Lambda$:

The $\arg\!\max$ argument on the right-hand-side eq. 25 is a directional derivative. It is maximal at parameter increments $\left[\Delta\eta,\,\Delta\Lambda\right]$ in the direction of gradient $\nabla\text{Pr}_{D}$ $=$ $\left[\frac{\partial\text{Pr}_{D}}{\partial\eta},\,\frac{\partial\text{Pr}_{D}}{\partial\Lambda}\right]$, where $\text{sign}\left\{\frac{\partial\text{Pr}_{D}}{\partial\eta}\right\}$ $=$ $-1$ and $\text{sign}\left\{\frac{\partial\text{Pr}_{D}}{\partial\Lambda}\right\}$ $=$ $+1$, and only increment $\Delta\Lambda$ depends on $\psi_{+1}$. Moreover, the quadratic forms in 27 are both strictly positive, since $\left[\boldsymbol{H}^{\text{T}}\boldsymbol{H}\right]^{-1}$ is positive definite. This collective system structure implies that the bracketed difference between the geometric factors is always non-negative (e.g., adding more sensors cannot reduce the source screening probability), so that:

The last term of eq. 28 defines the objective function for Rayleigh wave sampling. It omits explicit dependence on data $\boldsymbol{\mathcal{R}}$. Therefore, the optimal strike-relative deployment azimuth estimate $\hat{\psi}_{+1}$ for a receiver when $\sigma_{M}^{2}$ is unknown depends only on current receiver placement and the relative seismic wave structure ($\alpha$ and $\beta$). Other values of sensor placement where $\boldsymbol{A}\left[\boldsymbol{H}_{+1}^{\text{T}}\boldsymbol{H}_{+1}\right]^{-1}\boldsymbol{A}^{\text{T}}$ $<$ $1$ increase $\Lambda$ from its prior value (eq. 21). Auxiliary sensor azimuths where $\boldsymbol{A}\left[\boldsymbol{H}_{+1}^{\text{T}}\boldsymbol{H}_{+1}\right]^{-1}\boldsymbol{A}^{\text{T}}$ $>$ $1$ decrease $\Lambda$ from its prior value. To quantify the change in screening power, we reapply eq. 23 to estimate the updated threshold $\eta$ and reapply eq. 22 to compute source screening rates with $N+1$ sensors; we then update the decision rule (eq. 18).

Fig. 7a shows $\pi$-periodic optimal and sub-optimal solutions for extra receiver locations that supplement $16$ existing receivers uniformly deployed along a $75^{\circ}$ azimuthal arc (gray circles, Fig. 7b). The objective function shows local and global extrema of $\boldsymbol{A}\left[\boldsymbol{H}_{+1}^{\text{T}}\boldsymbol{H}_{+1}\right]^{-1}\boldsymbol{A}^{\text{T}}$. The maxima occur at deployment azimuths that minimize $\Lambda$ and therefore the probability of screening binary source types; these sensor locations provide the poorest improvement in source-type screening. Local minima (triangles) provide a locally optimal sensor placements near these worst-case sensor locations. These local minima solutions may apply to scenarios in which deployment outside restricted azimuthal ranges is impossible (e.g., inaccessible areas). Our global solutions (circles) confirm expectation that optimal receiver deployment should include points that include no current receiver coverage. These optimal deployment locations are not orthogonal to the mid-point of the azimuthal arc, but depend on the relative values of the body wave speeds (parameter $c$). Fig. 7b shows screening curves for the decision rule 18 when a receiver network contains the backbone of existing sensors (inset gray circles), supplemented with poorly placed sensors (squared one and two), versus optimally placed extra sensors (circled 1 and 2). The optimally placed sensor that supplements this $16$ sensor backbone network forms a $17$ sensor network the improves screening rates of the faulting source by $>6\times$ that of the backbone network, for sources with the greatest faulting signal shown (circled 1). An additional, optimally placed sensor (circled 2) provides a greater gain in screening performance, relative to an $18$ sensor network that includes two poorly placed auxiliary sensors (squared 2). The screening performance our test achieves with this azimuthally limited, 18 sensor arrangement outperforms average test that uses the $24$ sensor deployment in Fig. 6.

Lastly, we emphasize that the update to $\eta$ maintains a consistent false-attribution probability for the discriminant before and after adding any additional receivers. Failure to update this threshold when $\Delta\Lambda$ $<$ $1$ could lead an observer to erroneously conclude that adding receivers at particular values of $\psi_{+1}$ (shaded domain of Fig. 7(a)) can reduce the screening power of the same test.

We now quantify the largest faulting source signal (as defined by $\left(DS+SS\right)^{2}/\sigma_{M}^{2}$) that could be mistakenly attributed to an explosion, at a fixed false attribution probability $\text{Pr}_{FA}$. This application is particular important to underground explosion monitoring applications because it provides fundamental physical limits on the the discrimination power of eq. 16 tests to screen explosive from non-VDS faulting sources.

To compute this signal, we first note that the probability of mistaking a non-VDS faulting source for an explosive source is equivalent to the event that $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ $<$ $\eta$ when $\mathcal{H}_{1}$ is true (eq. 18). Here, threshold $\eta$ relates to the CDF under hypothesis $\mathcal{H}_{0}$ through eq. 23. The probability of miss-attributing a non-VDS faulting source to an explosion then equates to $1-\text{Pr}_{D}$ and relates to the faulting signal through a nonlinear equation that requires inversion of the noncentral-$F$ CDF for $\Lambda$:

The term $1-\text{Pr}_{D}$ is equivalent to a “miss” in waveform detection theory. The faulting signal SNR that solves eq. 29 relates to noncentrality parameter estimate $\hat{\Lambda}$ through:

Our estimate for the source term of $\hat{\Lambda}$, as measured by a given sensor deployment locations $\psi_{m}$ ($m=1,2,\cdots,N$) that are stored in $\boldsymbol{H}$, is:

Fig. 8 illustrates three solution curves to eq. 31 for a strike-slip faulting signal when $1-\text{Pr}_{D}$ $=$ $0.1$ (for a $0.9$ correct attribution probability), which depicts averages over $100$ random, azimuthally constrained sensor deployments. These solutions provide a lower bound on threshold source size, since $\Lambda$ is largest for strike-slip events ($\left(DS+SS\right)^{2}/\sigma_{M}^{2}$ $=$ $M_{0}^{2}/\sigma_{M}^{2}$). Vertical differences between distinct curves suggest that false attribution rates of a non-VDS faulting signal to an explosion source (via eq. 18) increases with azimuthal gap.

Fig. 9(a) better quantifies this gap effect on screening explosive from strike-slip sources. Each family of curves show the maximum, mean and minimum solutions to eq. 31 per azimuthal gap value. The vertical range of SNR values between each curve measures the variability that results from the $100$ randomly placed sensors. In particular, these solution curves indicate that strike-slip source size (SNR $=$ $M_{0}^{2}/\sigma_{M}^{2}$) increases exponentially as deployment gap increases beyond $90^{\circ}$ (equivalent to a $270^{\circ}$ azimuthal coverage). Moreover, the variability in such threshold estimates is asymmetric: the log-scale maximum threshold estimates (top curve) depart from the mean (middle curve) significantly more than do log-scale minimum threshold estimates (bottom curve). Noisy radiation patterns that associate with these solutions in Fig. 9(b) show that the decision defined by eq. 18 can mistake (with probability $1-\text{Pr}_{D}$ $=$ $0.1$) the subtle four-lobed shape for a noisy circular radiation pattern, even with full azimuthal sensor coverage.

We progress from Case I and now assume that a sensor network records a Rayleigh wave radiation pattern with unknown random process variability ($\sigma_{\mathcal{R}}^{2}$ in eq. 16), unknown faulting parameters ($\boldsymbol{\theta}$ in eq. 11), and a stochastic variance $\sigma_{M}^{2}$ that is either known, or has a significantly lower estimator variance than than that of $\hat{\sigma}_{\mathcal{R}}^{2}$. This case applies to idealized monitoring scenarios in which an underground testing complex colocates with a shallow fault of unknown geometry and radiation variability, but with well-characterized stochastic variability $\sigma_{M}^{2}$ that is perhaps quantified from historical records of proximal explosions. Sensor networks that record a low magnitude seismic event that originates from this test site faces the same challenge as that in Case I, but the screening statistic has a distinct form. We treat Case II in less detail than Case I because it appears less practically applicable than Case I while remaining more complicated. Appendix C derives several results that we state here. For Case II, the log-GLRT equivalent to eq. 17 is:

where eq. 14 and eq. 15 define $f_{\boldsymbol{\mathcal{R}}}\left(\boldsymbol{\mathcal{R}};\mathcal{H}_{0}\right)$ and eq. 13 defines $f_{\boldsymbol{\mathcal{R}}}\left(\boldsymbol{\mathcal{R}};\mathcal{H}_{1}\right)$. The MLE for $\sigma_{\mathcal{R}}^{2}$ is (algebra omitted):

in which the text that follows eq. 19 defines the projector matrix $\boldsymbol{P}_{\boldsymbol{H}}^{\perp}$. Our eq. B.2 defines the MLEs $\hat{\boldsymbol{\theta}}_{k}$ for parameter vectors $\boldsymbol{\theta}_{k}$ ($k=0,1$). We input these parameters to eq. 33 and reuse the symbol $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ to write the resultant, equivalent GLRT statistic (algebra omitted):

We use the MLEs for $\hat{\boldsymbol{\theta}}_{k}$ to express $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ as a sum of two statistically independent random variables ($X$ and $Y$):

We interpret the first quadratic form ($\|\boldsymbol{P}_{\boldsymbol{H}}^{\perp}\boldsymbol{\mathcal{R}}\|^{2}//\sigma_{M}^{2}$) in eq. 35 as the residual radiation pattern energy not attributable to the Rayleigh wave model, relative to random noise variance. The second term ($\|\boldsymbol{P}_{\boldsymbol{X}}\boldsymbol{\mathcal{R}}\|^{2}/\sigma_{M}^{2}$) is energy due to non-VDS faulting, relative to the same noise variance. The decision rule that tests the function of these two terms in eq. 35 is:

The threshold $\hat{\eta}$ in eq. 36 is strictly an estimate that maintains a false-attribution probability $\text{Pr}_{FA}$, consistent with Case I. We describe our estimation strategy for this threshold in Appendix C (see eq. C.5). To compute the screening performance for Case II (SNIR versus $\text{Pr}_{D}$), we note that the terms $X$ and $Y$ are independent, piecewise invertible functions. In particular, the term $X$ is piecewise invertible in argument $s$ $=$ $\|\boldsymbol{P}_{\boldsymbol{H}}^{\perp}\boldsymbol{\mathcal{R}}\|^{2}/\sigma_{M}^{2}$ since function $x(s)$ $=$ $s$ $-$ $N\ln(s)$ monotonically decreases over $s$ $<$ $N$ and monotonically increases when $s$ $>$ $N$. The ratio $s$ $\sim$ $\chi_{N-3}^{2}(0)$, where $\chi_{N-3}^{2}(0)$ is chi-squared distribution function with $N-3$ degrees of freedom and a zero noncentrality parameter; the distributional form of $s$ is identical for both circular and non-circular radiation patterns. Unlike the Case I test statistic, these PDFs are well defined if $N\geq 4$ sensors rather than eight sensors. The ratio $Y$ $\sim$ $\chi_{1}^{2}(\Lambda)$ is a chi-squared distribution function with $1$ degrees of freedom and a nonzero noncentrality ($\Lambda$ $>$ $0$) under $\mathcal{H}_{1}$. We use these results to compute the PDF of the sum $Z$ $=$ $X$ $+$ $Y$ through a variable transformation on the original, statistically independent random variables. Given that $X$ and $Y$ respectively have PDFs $f_{X}\left(x\right)$ and $f_{Y}\left(y\right)$, the PDF $f_{Z}\left(z;\mathcal{H}_{i}\right)$ for $Z$ is their convolution $f_{X}\left(x\right)\ast f_{Y}\left(y\right)$ under hypothesis $\mathcal{H}_{i}$ ($i=0,1$). Omitting the explicit algebra, we symbolize the PDF $f_{Z}\left(z;\mathcal{H}_{i}\right)$ as:

We compute eq. 37 in the frequency domain with the characteristic function method (eq. C.4). The PDF $f_{\chi_{1}^{2}}(\bullet)$ includes a finite noncentrality parameter under the alternative hypothesis ($i=1$). The theory we detail in Appendix C demonstrates that this parameter equates to the same scalar $\Lambda$ in eq. 21 that we derived from Case I. This means that $\Lambda$ completely quantifies the screening capability of $L_{\boldsymbol{\theta}}\left(\boldsymbol{\mathcal{R}}\right)$ to test between $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$. The PDF $f_{X}\left(x\right)$ is identical under $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$. We interpret its effect on $f_{Z}\left(z;\mathcal{H}_{i}\right)$ as a smoothing window in the convolution operation in eq. 37 that is independent of the Rayleigh wave radiation pattern model.

Fig. 10 shows the screening power of the Case II test statistic compared against that of the Case I test statistic, for three source types and a fixed deployment constraint (an azimuthal gap of $90^{\circ}$) for $N=12$ sensors. The graph format is identical to that shown by Fig. 4, Fig. 5 and Fig. 6, in which each curve represents an average (eq. 24) over 100 deployment realizations. Each such realization effectively recomputes the deployment factor $\left(\boldsymbol{A}\left[\boldsymbol{H}^{\text{T}}\boldsymbol{H}\right]^{-1}\boldsymbol{A}^{\text{T}}\right)^{-1}$ of $\Lambda$ and re-convolves the PDFs $f_{Y}\left(\bullet\right)$ and $f_{X}\left(\bullet\right)$ in the frequency domain as a spectral product. We document the additional numerical details that we implemented to construct the screening curves shown by Fig. 10 in text following eq. C.5.

Importantly, these curves demonstrate that the Case II decision rule provides a superior source-screening capability, at least for the same faulting signal SNIR and thresholds required to maintain a false attribution probability of $\text{Pr}_{FA}$ $=$ $10^{-3}$. We expect the Case II test statistic to provide such an improved, average performance because the Case I statistic additionally requires an MLE of the noise variance ($\hat{\sigma}_{M}^{2}$) under $\mathcal{H}_{0}$; in general, test statistics with unknown parameters show poorer screening performance than competing test statistics with fewer unknown parameters (Kay, , 1998, pg. 195). The variability in screening performance with sensor placement over permissible azimuths (note the $90^{\circ}$ gap), however, is comparable between both Case I and Case II (note shading). Both tests are also similar in the relative order of screening rates with focal mechanism. In particular, the Case II statistic shows a qualitatively similar ordering in screening probability when compared to the Case I statistic. The Case II decision rule (eq. 36) that tests radiation pattern data that resembles the Rock Valley earthquake (for example) shows a much higher probability of correctly screening such events from isotropic sources than, say, data that records an event that resembles the 2003-03-13 Lop Nor earthquake.