Application of subspace methods to detect and characterize coal mine related seismicity in the Western United States

An approach for subspace detection and magnitude estimation of small seismic events is proposed. The process is used to identify mining related seismicity from a surface coal mine and an underground coal mining district, both located in the Western U.S. Using a blasting log and a locally derived seismic catalog as ground truth, the detector performance is assessed in terms of verified detections, false positives, and failed detections. Over 95% of the surface coal mine blasts and about 33% of the events from the underground mining district are correctly identified. The number of potential false positives are kept relatively low by requiring detections to simultaneously occur on two stations. Many of the potential false detections for the underground coal district are genuine events missed by the local seismic network, demonstrating the usefulness of regional subspace detectors in augmenting local catalogs. A trade-off in detection performance between stations at smaller source-receiver distances, which have increased signal to noise ratios, and stations at larger distances, which have greater waveform similarity, is observed. The increased detection capabilities of a single higher dimension subspace detector, compared to multiple lower dimension detectors, are explored in identifying events that can be described as linear combinations of training events. In this data set, such an advantage can be significant, justifying the use of a subspace detection scheme over conventional correlation methods.

1 APPLICATION OF SUBSPACE METHODS TO DETECT AND CHARACTERIZE COAL MINE RELATED SEISMICITY IN THE WESTERN UNITED STATES by Derrick James Allen Chambers A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science Department of Mining Engineering The University of Utah August 2015 2 Copyright © Derrick James Allen Chambers 2015 All Rights Reserved vi The University of Utah Graduate School STATEMENT OF THESIS APPROVAL The thesis of Derrick James Allen Chambers has been approved by the following supervisory committee members: Michael K. McCarter , Chair 06/05/15 Date Approved Keith D. Koper , Co-chair 06/05/15 Date Approved Michael G. Nelson , Member 06/05/15 Date Approved and by Michael G. Nelson , Chair/Dean of the Department/College/School of Mining Engineering and by David B. Kieda, Dean of The Graduate School. vi ABSTRACT An approach for subspace detection and magnitude estimation of small seismic events is proposed. The process is used to identify mining related seismicity from a surface coal mine and an underground coal mining district, both located in the Western U.S. Using a blasting log and a locally derived seismic catalog as ground truth, the detector performance is assessed in terms of verified detections, false positives, and failed detections. Over 95% of the surface coal mine blasts and about 33% of the events from the underground mining district are correctly identified. The number of potential false positives are kept relatively low by requiring detections to simultaneously occur on two stations. Many of the potential false detections for the underground coal district are genuine events missed by the local seismic network, demonstrating the usefulness of regional subspace detectors in augmenting local catalogs. A trade-off in detection performance between stations at smaller source-receiver distances, which have increased signal to noise ratios, and stations at larger distances, which have greater waveform similarity, is observed. The increased detection capabilities of a single higher dimension subspace detector, compared to multiple lower dimension detectors, are explored in identifying events that can be described as linear combinations of training events. In this data set, such an advantage can be significant, justifying the use of a subspace detection scheme over conventional correlation methods. vi This work is dedicated to Amy Jo, who encourages me to strive for excellence vi TABLE OF CONTENTS ABSTRACT ................................................................................................................... iii LIST OF TABLES ........................................................................................................ vii LIST OF FIGURES ...................................................................................................... viii ACKNOWLEDGEMENTS ............................................................................................ ix Chapters 1. INTRODUCTION ..................................................................................................... 1 2. PRODUCTION BLASTS FROM AN OPEN PIT COAL MINE ............................... 5 2.1 Seismic Date from the Earthscope Transportable Array .................................. 5 2.2 ANF Seismic Catalog ..................................................................................... 6 2.3 Log of Mine Blasts ......................................................................................... 6 2.4 Subspace Construction ................................................................................... 7 2.5 Setting a Detection Threshold ........................................................................ 8 2.6 Detection and Event Association .................................................................. 10 2.7 Magnitude Estimation of Detected Events .................................................... 11 2.8 Results of Subspace Detection ..................................................................... 13 3. INDUCED SEISMICITY FROM AN UNDERGROUND COAL MINING DISTRICT ............................................................................................................ 22 3.1 Seismic Data from the TA ............................................................................ 22 3.2 ANF Seismic Catalog ................................................................................... 22 3.3 OMSHR Seismicity Catalog ......................................................................... 23 3.4 Subspace Construction, Threshold Determination, and Association .............. 23 3.5 Magnitude Estimation of Detected Events .................................................... 24 3.6 Results of Subspace Detection ..................................................................... 24 4. DISCUSSION ......................................................................................................... 31 4.1 Assessment of Detector Performance ........................................................... 31vi 4.2 Comparison of 2D Subspace Detector with two 1D Detectors ...................... 34 5. CONCLUSIONS ..................................................................................................... 41 APPENDIX: THEORY REVIEW ................................................................................. 43 REFERENCES.............................................................................................................. 47 vi LIST OF TABLES 2.1 Surface mine subspace construction summary ....................................................... 15 2.2 Surface mine subspace detection results ................................................................. 16 3.1 Underground mining district subspace construction summary ................................ 26 3.2 Underground mining district subspace detection results ......................................... 27 vi LIST OF FIGURES 1.1 Location of mines, TA stations, and ANF events ..................................................... 4 2.1 Subspace creation and detection threshold determination ....................................... 17 2.2 Detection association ............................................................................................. 18 2.3 Magnitude estimations for surface mine events ...................................................... 19 2.4 Subspace detector performance for the surface mine .............................................. 20 2.5 Example of an undetected event ............................................................................. 21 3.1 Magnitude estimations for the underground coal district events ............................. 28 3.2 Subspace performance for the underground coal district ........................................ 29 3.3 Spatial relations of subspace detections on Q21A .................................................. 30 4.1 Conceptualization of subspace detection and waveform correlation ....................... 39 4.2 Pit locations for surface coal mine ......................................................................... 40 vi ACKNOWLEDGEMENTS I would like to thank Dr. Kim McCarter for giving me the opportunity to study mining related seismicity, his patience as I explored various thesis topics, and for his genuine interest in his students. I am grateful for the guidance provided by Keith Koper and Kris Pankow, both of whom I greatly respect and admire. It has been a privilege to work with my very talented peers Tex Kubacki, Shawn Boltz, Jessica Wempen, and Jared Stein. I would also like to acknowledge the support of NIOSH for not only funding my research but providing the excellent mentorship of Pete Swanson and Jeff Whyatt. The conclusions expressed in this paper, however, are those of the author and does not represent the opinions or policies of NIOSH. 1 CHAPTER 1 INTRODUCTION Matched filters have a wide application in the analysis of seismic events with nearby locations because such events tend to produce highly similar waveforms at receivers that are located much farther away than the interevent distance (Geller and Mueller 1980; Schaff and Richards 2004; Dodge and Walter 2015). Waveform correlation, a common type of matched filter, has been used for many purposes in seismology, including improved detection of foreshocks and aftershocks (Peng and Zhao 2009; Schaff 2010), high-precision relocation of seismicity (Shearer 2002; Richards et al. 2006), monitoring of nuclear weapon test sites (Gibbons and Ringdal 2006; Zhang and Wen 2015), measuring changes in Earth structure (Poupinet et al. 1984; Zhang et al. 2005), identifying nonvolcanic tremor (Shelly et al. 2007), and waveform similarity clustering (Houser et al. 2008). Subspace detection-a multidimensional generalization of waveform correlation (Harris 2006)-has proven useful for many of the same purposes as waveform correlation (Maceira et al. 2010; Harris and Dodge 2011; Harris et al. 2012; Barrett and Beroza 2014; Song et al. 2014). Although more complicated to implement, in many cases subspace detection has advantages over waveform correlation including reduced computational expense, improved ability to detect events described by linear 2 combinations of training events, and a greater tolerance for source-location variations, especially for temporally migrating sources (e.g., most aftershock sequences, Wesson 1987). Matched filters have been applied to mining induced seismicity (MIS, see Gibowicz (2009) for a review) for a variety of objectives. Several studies have focused on validating matched filter techniques with mining seismicity because mining events often have good ground truth (e.g., Joswig and Schulte-Theis 1993; Gibbons and Ringdal 2005; Slinkard et al. 2014). Another category of work seeks to use seismic cross correlation to understand geological conditions and stress states associated with mining. For example, Spottiswoode and Milev (1998) used cross correlation to locate and classify MIS into clusters in order to identify geological structures that could adversely affect mining activity. Grêt et al. (2006) showed in an experimental mine that correlating the coda of repeating sources could be used to infer stress states of mine pillars. Kubacki et al. (2014) used cross correlation to identify and locate previously undetected foreshocks and aftershocks from a large coal mine collapse, the results of which shed light on local structure and collapse conditions. In a similar fashion, Pankow et al. (2014) employed waveform correlation to find several previously unidentified earthquakes induced by a large surface mine slope failure. For other disciplines, mining seismicity has proven problematic by contaminating data sets of interest. In particular, earthquake monitoring organizations often find it difficult and time-consuming to discriminate surface blasts and MIS from tectonic earthquakes. Wüster (1993) and Stein et al. (2015) found waveform similarity clustering to be one effective measure for successfully discriminating mining related seismicity 3 from tectonic earthquakes. Likewise, discrimination between large mining blasts and small potential nuclear explosions can be problematic in nuclear test ban treaty monitoring, but waveform characterization of blasting events, and by extension matched filters, can aid in discrimination efforts (Smith 1993). Moreover, Ferretti et al. (2005) demonstrated that clustering based on waveform similarity is useful in classifying many families of seismic events. In this study, a subspace procedure (Harris 2006; Harris and Paik 2006) is implemented to detect and characterize coal mine related seismicity in the western U.S. The subspace detection method is applied to data from two distinct mining environments, an open pit coal mine in southwestern Wyoming in which surface blasts are routinely carried out, and an underground coal mining district in western Colorado in which MIS is prevalent (Figure 1.1). Importantly, in both of these cases independent information about the local seismicity is available for use as ground truth in evaluating the performance of the detectors. The overall performance of both single station and two station detectors is assessed, the effect of source-receiver distance on detector effectiveness is quantified, and the cause of failed detections of known events is examined. Additionally, the effectiveness of a two-dimensional subspace detector relative to two one-dimensional subspace detectors (i.e., conventional waveform correlators) is evaluated. 4 Figure 1.1 Location of mines, transportable array (TA) stations, and ANF events. Map includes mining district locations (red and green borders), TA stations used in this study (white inverted triangles), mine permit boundaries (blue polygons), and seismic events in the Array Network Facilities (ANF) catalog (dots colored according to local magnitude). 5 CHAPTER 2 PRODUCTION BLASTS FROM AN OPEN PIT COAL MINE 2.1 Seismic Data from the Earthscope Transportable Array Station coverage near the open pit coal mine in southwestern Wyoming was temporarily increased in the years of 2007 to 2009 when stations from the Earthscope Transportable Array (TA) experiment were deployed in the region. The two nearest TA stations were M17A (Evanston, Wyoming), located approximately 20 km south of the mine, and M18A (Lyman, Wyoming), located about 60 km southeast of the mine (Figure 1.1). Both M17A and M18A recorded three-component, broadband data at 40 samples per second. For use in the subspace detector, hour-long segments of continuous data from these two stations for the entire year of 2008 are downloaded using the Incorporated Research Institutions for Seismology (IRIS) web services (www.ds.iris.edu). Most of the continuous data files are complete and without issues, but a few files, less than 5%, have gaps, clipped waveforms, or missing channels. The hour-long files with such problems are discarded. After assessing data quality, the instrument response is removed and the data are filtered between 1 and 10 Hz. 6 2.2. ANF Seismic Catalog Although the mine is located within the authoritative region monitored by the University of Utah Seismograph Stations (UUSS, www.seis.utah.edu), the most detailed seismicity catalog for the mine in 2008 is that produced by the Array Network Facility (ANF, www.anf.ucsd.edu), which manages the TA data. This is because the UUSS seismicity catalog is intended to be an earthquake catalog and suspected mine blasts are purposely excluded, while the ANF catalog includes all seismic events. The ANF processing and location procedures are described in Astiz et al. (2014). The ANF catalog is searched for events located in proximity to mine permit boundaries that are also recorded on both M17A and M18A (Figure 1.1). Eighty-one events ranging in magnitude from 1.74 to 2.88 ML (local magnitude) and occurring in late 2007 to early 2009 meet these criteria. The magnitude of completeness of the ANF catalog is approximately 2.2 ML for the local daylight hours during which blasting is permitted (see Figure 13 of Astiz et al. 2014). For each station and event, five minutes of three-component seismic data are download starting one minute before the reported origin time in the ANF catalog. The instrument responses are removed and the data filtered between 1 and 10 Hz. These data are used to construct subspace detectors. 2.3. Log of Mine Blasts In many mining operations, including the surface mine in this study, the use of chemical blasts is necessary to facilitate extraction and comminution of the resource (Darling 2011). A log of all the blasts initiated by the mine in 2008 was provided to the authors by a blasting engineer. The log indicates in which of four pits each shot 7 occurred, the approximate detonation time, the amount and type of explosive used, and a few other parameters. For the year 2008, there are 381 blasts recorded in the log, 50 of which also appear as seismic events in the ANF catalog. In order to verify the accuracy of the blast log, it is first confirmed that the smallest blasts in terms of charge-weight (pounds of explosive detonated in an 8 ms period) from all four pits are in fact visible at the nearest station, M18A. Missing or misrecorded blasts are accounted for by examining the seismic records of M18A and removing the blasts from the log that produce no visible signal or that occur in times in which the continuous data are discarded due to quality issues. Eleven such cases are discarded from the log. Next, all but the largest blasts that were reportedly initiated at the same time are discarded resulting in the removal of 13 more log entries. This second redaction is performed because the subspace algorithm used in this study only declares a detection every few seconds in order to avoid double counting events. The remaining 357 log entries are then used to verify detections and assess the number of failed detections and false positives returned by the detector. 2.4. Subspace Construction A simplified workflow, similar to the one outlined in Harris (2006), is followed in order to construct subspace detectors from the 81 ANF events-which are assumed to be mine blasts-that locate close to the surface mine. The five-minute data from all three components are multiplexed to form a waveform that includes all observable phases. Similarity matrices are then created independently for each station by cross correlating all the multiplexed waveforms (Figure 2.1A). The events are grouped together based on 8 dissimilarity using a single link, hierarchical, clustering algorithm, and dendrograms are created to visualize the groupings (Figure 2.1B). The required correlation coefficient to define the clusters is allowed to vary between stations in order to avoid forming an excessive number of small clusters, and to make the grouping of events about the same for both stations. For each cluster, the lowest correlation coefficient of the linking pairs is used as a measure of cluster "tightness." The higher this quantity is, the higher the overall similarity of the cluster is; for a cluster of identical waveforms the lowest correlation coefficient of the linking pairs would be 1.0. For each of the clusters, the following processing is carried out. The waveforms are aligned and trimmed to a length of 30 seconds (Figure 2.1C) starting from a user-defined pick time, which ideally corresponds to the first arrival. Singular value decomposition (SVD) is performed on the aligned waveforms in order to calculate an orthonormal basis (Figure 2.1D) that spans the training events. The minimum numbers of vectors from the SVD bases required to represent the groups of waveforms adequately are estimated by calculating the fractional energy captured as a function of dimension of representation (i.e., the dimension of the subspace; Figure 2.1E). The fewest singular vectors able to capture an average of 90% of the waveform energy are selected as the subspace (Table 2.1; M18A subspace 2 corresponds to Figure 2.1). 2.5. Setting a Detection Threshold The detection statistic is the ratio of vector energy projected into a subspace over the original vector energy, and therefore can vary from 0.0 to 1.0 (Appendix A). Like the correlation coefficient, the detection statistic is used to test the hypothesis that a vector 9 (i.e., a segment of continuous data) contains signal (similar to the training events of the subspace) and noise, or noise only. Just as in cross correlation, the value of the detection statistic deemed significant depends strongly on the product of the frequency band and the duration of the detection window (Schaff 2008). The dimension of representation also has a significant effect; higher dimension detectors generally produce higher detection statistic values. The following procedure is used to estimate detection thresholds. An empirical estimation of the detection statistic resulting from seismic noise (i.e., recorded seismic data with no apparent transient signals) is made by calculating the detection statistic at each time step, for each subspace-station pair, using 100 hours of continuous data. Only hours of continuous data with no value above eight times the root-mean-square amplitude of the entire hour are used in an attempt to avoid transients. A beta distribution is then fit to the resulting detection statistic values using the maximum likelihood method. A value of the detection statistic that sets the upper tail probability of the beta distribution to 10-12 is selected as the threshold for declaring detections (Figure 2.1F) for each subspace-station pair (Table 2.1). If the distribution describes the data perfectly, any detection statistic value above the prescribed threshold would have a 10-12 likelihood of being a false detection due to incoherent random noise. However, the misfit in the tail of the distributions (Figure 2.1F and Kubacki et al. (2014)) is a significant and persistent problem with using statistical models to set detection thresholds. The misfit is generally thought to be due to transient signals that are more similar to the detector than noncoherent noise, but that do not originate from the same area as the training events and are therefore considered false 10 detections (Slinkard et al. 2014). Later on, strategies for dealing with this issue are discussed. A sliding window scheme is applied to the continuous data to calculate the detection statistic for every time step for all available continuous data of each station-subspace pair. The absolute times of detections, detection statistic value, station, subspace, estimated magnitude, signal to noise ratio (SNR), and predicted origin times are then recorded. 2.6 Detection and Event Associations Using the difference between a user-defined onset time for all of the events in the subspace and the ANF event origin times, a maximum and minimum offset time, essentially travel time bounds, are calculated for every subspace (Figure 2.2). The extents of maximum and minimum offset times are used to estimate a range of possible origin times for each detection. Overlaps in the range of predicted origin times are used to associate detections from different stations together. Additionally, the association window is widened by 1.0 second in order to mitigate missed detections of events that may have differential arrival times slightly outside those of the subspace training events, and to account for potential origin time error in the ANF catalog. The blast initiation times in the log were recorded only to about 5-minute accuracy so any detection occurring either 5 minutes before or 5 minutes after the recorded blasting time is counted as verified. If more than one detection on the same station meets this requirement for a blasting event, then only the detection with the largest detection statistic is associated with the blast. Likewise, if a detection could be 11 associated with two blasts, then it is only assigned to the blast with the greatest charge-weight. 2.7. Magnitude Estimation of Detected Events A relative magnitude estimation method using SVD output has been demonstrated by Rubinstein and Ellsworth (2010). This method allows for very little variation in source-location, so it is not as useful for subspace detection applications where the objective is to find as many events as possible, especially those with significant variations in source-location parameters. Schaff and Richards (2014) show that simple waveform scaling, when used to estimate relative magnitudes between two events, often introduces significant bias from variation in the waveforms. In order to calculate a more accurate scaling factor, Schaff and Richards (2014) suggest using the log of an L2 norm ratio of the two waveforms, which is a nearly bias-free indicator when signal to noise ratios are not too low. In the context of subspace detection, this method provides a much better estimate of magnitude than waveform scaling for events that have significant signal energy that is not captured by the subspace, 𝑈. More precisely, the method accounts for energy in an unknown orthogonal signal subspace, 𝑈⊥, in magnitude estimates. This discrimination is important because the energy in 𝑈⊥ can be significant, at times exceeding the energy in 𝑈. However, the L2 norm method can suffer from a potential SNR bias because it makes no differentiation between energy in 𝑈⊥ and energy of seismic noise. As a result, magnitudes are overestimated when the SNR of the continuous data is too low; in the data used in this study, the SNR bias is negligible. Two relatively simple methods of estimating the magnitudes of the detected 12 events are explored. The first is a subspace adaptation of waveform scaling. The magnitude (𝑀𝑖) of detection 𝑖 is estimated by first calculating a scaling factor using the ratio of energy projected into the current subspace (𝑈) of the continuous multiplexed data that triggered the detector (𝐶𝑖) and each of the training events used to form the subspace (𝑇𝑗). The scaling factors is used to estimate a magnitude based on the corresponding reference magnitude reported in the ANF catalog (𝑚𝑎𝑔𝑗). The result is then weighted by the square of the normalized correlation coefficient (𝑐𝑐𝑖𝑗) between 𝑇𝑗 and 𝐶𝑖: 𝑀𝑖=Σ(log10(√(𝑈𝑇𝐶𝑖)𝑇𝑈𝑇𝐶𝑖(𝑈𝑇𝑇𝑗)𝑇𝑈𝑇𝑇𝑗⁄)+𝑚𝑎𝑔𝑗)𝑐𝑐𝑖𝑗2𝐽𝑗=0 Σ𝑐𝑐𝑖𝑗2𝐽𝑗=0 (2.1) The second method, similar to the L2 norm method, is a weighted standard deviation (std) ratio scheme of the form: 𝑀𝑖=Σ(𝑙𝑜𝑔10(𝑠𝑡𝑑(𝐶𝑖)𝑠𝑡𝑑(𝑇𝑗)⁄)+ 𝑚𝑎𝑔𝑗)𝑐𝑐𝑖𝑗2𝐽𝑗=0 Σ𝑐𝑐𝑖𝑗2𝐽𝑗=0 (2.2) Because the most significant variables in predicting peak particle velocity caused by a blast is the distance from the blast to the point of interest and the charge-weight (ISEE 2011), the log base ten of the charge-weight and the estimated local magnitude should be approximately linearly related. The standard error of linear regression-the root mean 13 square of the residuals-between the log of charge-weight and estimated local magnitude is used to evaluate the two estimation schemes. The standard error for Equation 2.1 is 0.047 and for Equation 2.2 is 0.038 (Figure 2.3), indicating the standard deviation method provides a better estimate. Further evidence for this conclusion is presented in Section 3.4. For both data sets, the squared weighting scheme of the correlation coefficient (𝑐𝑐𝑖𝑗) outperforms other basic weighting schemes in minimizing the standard deviation of the misfit function, or the standard error of the regression line, for both Equations 2.1 and 2.2. The other weighting schemes evaluated include: only using the best correlated template (i.e., the 𝑖-𝑗 pair with the highest 𝑐𝑐𝑖𝑗), linear weighting, and cubic weighting of 𝑐𝑐𝑖𝑗. It should be noted that both Equations 2.1 and 2.2 require waveforms to be fairly well windowed, and that the source-receiver distance not vary greatly between events 𝑖 and 𝑗. In the context of subspace detection, however, good event windowing is already required, and only events relatively near the training events will trigger the detector at sufficiently conservative detection statistic thresholds or if source-locations are constrained by use of a multistation detector. Also, both Equation 2.1 and 2.2 depend on good reference magnitudes as any error in training event magnitudes propagates to the magnitude estimate of the detected event. 2.8. Results of Subspace Detection Nearly all the blasts are identified on both stations. Using the two-station detector, a total of 345 blasting events are detected out of 357 (~97%) recorded in the 14 revised log. There are 11 detections that occur during local business hours and do not correspond to blasts in the log (Figure 2.4). These are classified as "unknown" (UK), because they could be blasts that are misreported or unrecorded. By combining the two stations, the number of sure false detections, or events occurring outside local business hours, is reduced to zero. When using only a single station, the number of false detections is as high as 329 (Table 2.2). Of the 11 shots not detected by both stations, nine are successfully detected by M18A (the more distant station) but not by M17A. For these nine events, spikes in the detection statistic vector are visible on M17A but they do not reach the detection threshold (Figure 2.5). Using lower acceptable probability of false detection in the subspace construction does allow more events to be detected, with the maximum of 356 of the 357 events being detected on M18A (Table 2.2), but the potential false detections become considerably higher. 15 Table 2.1 Surface mine subspace construction summary M17A Subspace # of Events Dimension of Rep. Lowest CC* DS Threshold˟ 1 13 8 0.40 0.141 2 29 20 0.41 0.256 3 28 15 0.44 0.215 4 8 4 0.48 0.079 M18A Subspace # of Events Dimension of Rep. Lowest CC DS Threshold 1 28 11 0.58 0.088 2 8 4 0.59 0.054 3 28 11 0.60 0.086 4 11 5 0.61 0.061 * CC = Correlation Coefficient ˟ DS = Detection Statistic 16 Table 2.2 Surface mine subspace detection results Station False Detections Unknowns Verified Detections Missed Events Probability of False Detection (Pf) = 10-8 M17A 2709 985 352 5 M18A 1435 2374 356 1 Both 7 39 352 5 Probability of False Detection (Pf) = 10-10 M17A 890 306 348 9 M18A 633 1155 355 2 Both 1 18 348 9 Probability of False Detection (Pf) = 10-12 M17A 329 124 345 12 M18A 322 600 355 2 Both 0 11 345 12 Probability of False Detection (Pf) = 10-14 M17A 62 28 340 17 M18A 104 211 351 6 Both 0 10 337 20 17 Figure 2.1 Subspace creation and detection threshold determination. Panel A shows the correlation coefficient of each event pair organized by similarity grouping using the 81 ANF events recorded at M18A. Diagonal elements are autocorrelations. Panel B shows the hierarchical clustering of the waveforms based on dissimilarity (four clusters, each colored differently), with 6 events left as unclustered singletons. Panel C shows the 30-second time window of three component multiplexed data, after being aligned and trimmed to 30 seconds (waveforms between green vertical lines), belonging to the red cluster (subspace two) in B. Panel D contains the orthonormal basis vectors that are selected (blue) and unused (grey). Panel E shows the criteria for selecting dimension of representation (how many basis vectors to use, green vertical line) requiring the average (red line) of the fractional energy captured of each event (grey dotted lines) to be at least 0.9. Panel F depicts the selection of a detection threshold (green vertical line) based on a beta distribution (black line) fitted to detection statistic values of the subspace with continuous data with no high amplitude signals (blue histogram). 18 Figure 2.2 Detection association. Aligned multiplexed waveforms from subspace 2 (red group Figure 2.1B) of station M18A. The maximum and minimum difference between onset time of the aligned waveforms (green line) and corresponding origin time of each event (red stars; event labeled 1 defines maximum offset time and event labeled 3 defines minimum offset time) are used to predict an origin time range (red lines) for each detection. Detections are then associated together across a network based on overlap between minimum and maximum estimations of origin times. 19 Figure 2.3 Magnitude estimates for surface mine events. Magnitude estimates (y axis) from Equation 2.1 (top) and 2.2 (bottom) are compared against charge-weight (x axis) as reported in the blasting log. Standard error (S) values of linear regression indicate a better fit for Equation 2.1. Estimated magnitudes are based on scaling the corresponding magnitudes in the ANF catalog and are therefore local magnitudes (ML). 20 Figure 2.4 Subspace detector performance for the surface mine. The number of events from each pit detected by each subspace for the combined station results are shown. Events labeled UK are unknown and FD are false detections. 21 Figure 2.5 Example of an undetected event. The two-station detector fails to identify this known event. The panels show the waveform of the continuous data in black, and the detection statistic vector of the four subspaces in blue for station M17A (left) and M18A (right), respectively. The dotted black line is the threshold of each subspace. 22 CHAPTER 3 INDUCED SEISMICITY FROM AN UNDERGROUND COAL MINING DISTRICT 3.1. Seismic Data from the TA Stations Q20A (Grand Junction, Colorado) and Q21A (Paonia, Colorado) of the TA network are used to study seismicity originating from the underground, longwall coal mining district. Q20A and Q21A were located approximately 50 km and 10 km from mining activity, respectively (Figure 1.1). Again, three-component broadband data sampled at 40 Hz are downloaded in hour segments and filtered between 1 and 10 Hz. Files with any gaps, clipped waveforms, or missing channels are discarded. Less than 5% of the data are either missing or rejected. 3.2. ANF Seismic Catalog There are 92 events in the ANF catalog that occurred from late 2007 to early 2009 that were recorded on Q20A and Q21A. All 92 are used in the subspace construction (although only 80 of the events occur in 2008). The local magnitudes (ML) range from 1.2 to 2.75. This portion of Colorado is also routinely monitored by the USGS National Earthquake Information Center (NEIC) in Golden, Colorado; however, the mission of NEIC is focused on moderate-to-large earthquakes. 23 3.3. OMSHR Seismicity Catalog The Office of Mine Safety and Health Research (OMSHR) provided the 2008 catalog of seismic events created using data from an agency-run local array deployed above the underground coal mining district. During 2008, the array consisted of 13 stations each with an Episensor, 3-component accelerometer and a vertical component 1 Hz L4-C seismometer. Digital waveforms were recorded at 100 Hz. The resulting catalog (NF catalog) was produced using the automatic detection and location algorithms of the Earthworm system (Johnson et al. 1996). The software was configured to perform well on a mine-wide scale (Swanson et al. 2008), but the NF catalog has not undergone quality review by analysts. Helicorder images of the continuous data recorded by each station were also provided, which were used to quantify how many stations produced healthy-looking continuous data for every 6-hour period in 2008 in order to assess network health. For verification, only NF events and subspace detections occurring during times when at least 11 of the 13 stations appear operational are considered. There are 9,010 NF events that meet this condition, ranging in duration magnitude (Md) from 0.0 to 3.1. 3.4. Subspace Construction, Determination of Threshold, and Association The method outlined in Section 2.4 is followed to form subspace detectors and the procedure outlined in Section 2.5 is followed to set detection thresholds. The results of subspace construction are similar for both Q20A and Q21A (Table 3.1). The detections are associated across stations in the same way described in Section 2.6. For detections 24 association with verified events, however, the procedure is slightly different. A detection is only associated with an NF event if the reported NF origin time falls within a second of the predicted origin time range for the detection. 3.5. Magnitude Estimation of Detected Events Magnitude estimates for each detection from January and February of 2008 are made using Equations 2.1 and 2.2. from Section 2.7. Although the magnitudes reported in the NF catalog are duration magnitudes (Md), they are similar to the local magnitudes (ML) reported in the ANF catalog for the training events found in both catalogs (Figure 3.1A and 3.1B). Because of this similarity, comparing the NF duration magnitudes and the predicted magnitudes (which use the ANF local magnitudes) is reasonable. For the events in the NF catalog that were identified by the subspace detector, the standard deviation method (Eq. 2.2; Figure 3.1C and 3.1D) clearly outperforms the projected energy ratios method (Eq. 2.1; Figure 3.1E and 3.1F). For this data set, Equation 2.2 yields a zero-mean distribution of misfit to two decimal places, but this lack of bias will not always be the case, especially when detections have low SNR. 3.6. Subspace Results for the Underground Coal Mining District For station Q21A, 6,567 of 9,010 (~73%) of the events in the NF catalog were identified with 42,578 potential false detections. Using the more distant station, Q20A, 2,954 (~33%) NF events were detected with 8,340 false detections. When both stations are used 2,923 (~32%) of the NF events are identified with 1,947 potential false detections (Table 3.2). For the two-station detector, if the detection threshold were 25 increased, then both the number of detected events and potential false detections would decrease at approximately the same rate (Figure 3.2A). Interestingly, some of the potential false detections have very high detection statistics. The 18 potential false detections with the highest estimated magnitudes (Md > 1.25) are all visually identifiable upon inspecting the NF helicorder records, indicating that the examined potential false detections are genuine events missed by the NF array. Additionally, since the potential false detection curve and the verified detection curve are nearly identical in shape (Figure 3.2A), most of the potential false detections are likely genuine events not recorded in the NF catalog. Assuming equal data quality for Q20A and Q21A, the source-receiver distance plays a large role in the magnitude of events the detectors are able to identify; the closer station (Q21A) was much more effective in finding events with magnitudes 0.5 to 1.25 (Figure 3.2B) than the more distant station (Q20A). The subspaces generally detect events that are closely related spatially (Figure 3.3). Note that Figure 3.3 only shows the results for the periods in which 11 NF stations were recording data, and therefore not all detections or ANF events are shown in the figure. 26 Table 3.1 Underground mining district subspace construction summary Q20A Subspace # of Events Dimension of Rep. Lowest CC* DS Threshold˟ 1 56 26 0.39 0.069 2 3 2 0.43 0.023 3 4 3 0.43 0.024 4 19 6 0.42 0.036 Q21A Subspace # of Events Dimension of Rep. Lowest CC DS Threshold 1 60 20 0.48 0.076 2 4 3 0.52 0.022 3 18 4 0.51 0.025 4 3 2 0.51 0.018 * CC = Correlation Coefficient ˟ DS = Detection Statistic 27 Table 3.2 Underground mining district subspace detection results Station Potential False Detections Verified Detections Missed Events Probability of False Detection (Pf) = 10-8 Q20A 23609 3571 5439 Q21A 82183 7073 1937 Both 3292 3571 5439 Probability of False Detection (Pf) = 10-10 Q20A 13108 3208 5802 Q21A 55514 6839 2171 Both 2458 3197 5813 Probability of False Detection (Pf) = 10-12 Q20A 8340 2954 6056 Q21A 42578 6567 2443 Both 1947 2923 6087 Probability of False Detection (Pf) = 10-14 Q20A 4523 2648 6362 Q21A 29609 5972 3038 Both 1461 2523 6487 28 Figure 3.1 Magnitude estimations for the underground coal district events. Panels A and B show a high similarity between the NF duration magnitude and the ANF local magnitudes. C and D show the magnitude estimation using Equation 2 (blue) and E and F show the estimates using Equation 1 (green). The reported magnitudes in the North Fork catalog (left panels, x axis) and the predicted magnitudes or ANF magnitudes with the two methods (left panels, y axis) are compared against a perfect fit (red line). The right panels show histograms of the misfit with the residual on the x axis and the occurrence on the y axis. 29 Figure 3.2 Subspace performance for the underground coal district. Panel A shows the number of verified detections (blue line), potential false detections (green line), and events in the NF catalog (black solid line) as a function of detection statistic for the two station detector with a probability of false detection set at 10-12. Panel B shows the percentage of NF events detected by magnitude for each station binned in increments of 0.25. The difference between the red line and the green line, and the blue line and the green line, show the events missed by the two station detector but detected by Q21A and Q20A, respectively. 30 Figure 3.3 Spatial relations of subspace detections on Q21A. Verified detections are plotted by color according to subspace using the NF locations.31 CHAPTER 4 DISCUSSION 4.1 Assessment of Detector Performance The two-station detector at regional distances performs very well in finding surface blasts (~97% of events identified, Figure 2.4), and moderately well for finding MIS produced by the coal mines (~32% of events identified, Figure 3.2). In both settings, when only a single station is used, and the detection thresholds are set via the method outline in Section 2.5, the number of false positives becomes very large. The large number of false positives is often apparent in the variation of the observed data from the tail of the fitted theoretical distribution (Figure 2.1F). The misfit is generally attributed to rogue seismic events-transient signals not originating from the same source/location as the training events-that share some small degree of similarity with the training events. This explanation is reasonable because events originating in different regions can sometimes be moderately well correlated (Dodge and Walter 2015). The effects of the undesired signals on the detection results can be mitigated or eliminated by only accepting detections that occur on more than one station (Slinkard et al. 2014). Multistation detectors produce far fewer false detections for two reasons: First, differential arrival time requirements can be imposed to constrain the source location to the same area the training events originated. Second, adding more complexity 32 to the detector greatly reduces the chance that noncoherent random noise will exceed the threshold. In this study, requiring detections to occur on two stations is successful in greatly reducing the number of false detections. The disadvantage in this case, however, is that the number of potential detections, and the detectable magnitude range, is limited to that of the poorest performing station (Figure 3.2B). Nine of the 11 missed events from the surface mine were not detected by M17A but were detected on the more distant station M18A. The closer station showed no advantage over the more distant station because the magnitude of most of the events (estimated using Equation 2.2; Figure 2.3) is above 1.25, large enough to be easily detected by both stations. Conversely, for the underground mining district, the closest station (Q21A) performed much better than the more distant station (Q20A) but there were events detected by Q20A and not by Q21A (Figure 3.2B; difference between blue line and green line). A station located farther away from a source region can be a better matched filter detector than a station located closer to the source region because the observed waveforms tend to be more similar at the more distant station. The increased similarity is due to several reasons, including a higher percentage of shared ray paths, closer focal sphere sampling, increased phase separation, and greater attenuation of higher frequency energy that tends to be more sensitive to geological variation. If the station is too far away, however, the SNR of the unidentified signals can become too low to be detected, particularly for low magnitude events (e.g., Q21A is better at detecting lower magnitude events than Q20A, Figure 3.2B). Therefore, there is tradeoff between stations with larger 33 source-receiver distance performing better by recording more similar waveforms and closer stations performing better due to higher SNR. Moreover, a detector's failure to identify an event could be caused by variation in source mechanisms from those characterized by the training events. For example, studies have found that surface blasts occurring in nearly the same location can produce very different waveforms (Bonner et al. 2003, McLaughlin et al. 2004) possibly caused by variation in the free-face orientation, hole delays, explosive agent, or a variety of other blasting practices (Stump et al., 2001). In a synthetic data set, Baisch et al. (2008) also found that varying source mechanisms, even when the events are collocated, can cause the resulting waveforms to decorrelate very quickly. The subspace detectors are able to detect several events with moderate magnitudes (Md > 1.25) that are not in the NF catalog, but were easily identifiable in the helicorder records, and so are not false detections. For the combined results of Q20A and Q21A, the number of verified event detections vs. detection threshold curve matched closely in shape to the number of possible false detections curve (Figure 3.2A). This implies that many, perhaps most, of the possible false detections occurring on the two station detector are genuine events from the mining district that were missed by the NF network. In this case, regional subspace detection is useful in augmenting the conventional seismic catalog by finding unrecorded events. In addition, the events detected by a subspace that are located far from most other events detected by the same subspace (Figure 3.3) are possibly mislocated. Such events could be flagged and subjected to extra scrutiny in order to improve the quality of locations in the catalog. 34 The spatial correlation between events detected by a given subspace provides a way to estimate epicentral locations, as well as to broadly classify the events into some pre-determined group, in this case to a particular pit or section of a mine, with data from as little as one station. As other studies have found, it is apparent in this study that event classification through subspace detection and waveform similarity can be a powerful tool in discriminating one class of event from another (e.g., MIS from earthquakes or chemical blasts) because waveform similarity is sensitive to both source type and location. Although the subspace methods are not able to identify all of the events in the blasting log or NF catalog, they are able to greatly increase the detection capabilities of the regional TA stations relative to conventional methods. The same subspace methods might be applied to local array data in order to further increase catalog completeness. If sufficient stations are used in the detection process, double difference locations (e.g., Waldhauser and Ellsworth 2000) of the newly detected events could be carried out. 4.2. Comparison of a 2D Subspace Detector with Two 1D Correlation Detectors In theory, a two-dimensional subspace detector is more effective than two one-dimensional subspace detectors (i.e., waveform correlators) in detecting seismic events best described by linear combinations of shared training events. To illustrate the concept, consider two training events that are approximately orthogonal, or, more realistically, are the resulting orthogonal bases from a SVD of the two training events. A two-dimensional subspace detector and two one-dimensional subspace detectors, with 35 window lengths of N samples, could be constructed from the training events. As the detectors slide across continuous data, each observation is a vector of N elements that is to be tested for signal. Each of the orthogonal training events can be visualized as a horizontal axis (event X and event Y for the X and Y axes in Figure 4.1) and the vertical (Z) axis can represent all orthogonal bases. In this case, each possible observation from the sliding window operation on the continuous data must plot somewhere on the unit sphere. Using only the absolute values of the X, Y, and Z coordinates every point will plot in the first octant, which is more convenient for visualization. The majority of observations would likely only contain noise with no coherent signal (Figure 4.1A, grey spheres) and, consequently, would cluster near the top of the sphere. Another population of observations would be the events that contain coherent signal as well as noise. Some of these may be similar to the master events, and others will not be, thus allowing these events to appear more uniformly on the sphere (Figure 4.1B, yellow spheres). The detection statistic between the one-dimensional subspace composed of event X (subspace X) and an observation is simply the square of the X value of the observation (the correlation coefficient would be the X value). The events detectable by subspace X can be represented by the observations in contact with a circle projected onto the unit sphere, with the origin at the X axis. The diameter of the circle would grow as the threshold is lowered. The detection statistic of the two-dimensional subspace (subspace X-Y) would be one minus the square of the Z value of the observation. The events detected by subspace X-Y are the observations in contact with a vertical cylinder projected onto the sphere, the height of which would grow as the detection threshold is lowered. If there exists an event that is well described by a linear combination of the two 36 horizontal bases, it will plot on the sphere between both horizontal axes (Figure 4.1C red sphere). In order for subspace X, or Y, to detect the event, the threshold would have to be lowered sufficiently to also return many false detections (Figure 4.1D, blue circular shell). Subspace X-Y, however, can successfully detect the event without including nearly as many false detection (Figure 4.1D, blue horizontal shell). As a result, the two-dimensional subspace has greater detection capabilities for this class of events than two one-dimensional detectors. In contrast, if the event of interest were located primarily on the X axis (Figure 4.1E), then subspace X would return less false detections than subspace X-Y (Figure 4.1F). Therefore, augmenting a subspace with more dimensions is only advantageous if the additional dimensions can meaningfully contribute to the description of the events of interest. Otherwise, the additional dimensions reduce the effectiveness of the detector. The same concept can be seen in Figure 9 of Harris (2009) and Figure 11 of Song et al. (2014) when the addition of dimensions to a subspace detector decreases the probability of successful detection. Note, however, that subspace X-Y is able to detect the event in Figure 4.1E with far fewer false detections than subspace X can detect the event in Figure 4.1C. This difference in false detections returned would decrease as more dimensions are added to the detector, but as long as the degree of representation of the subspace is much less than the effective dimension of the embedding space, the higher dimensional subspace would still have an advantage. Consequently, if the linear combination events exist, and their detection is important, a higher dimensional detector is a better choice than a lower one, as long as each dimension is significant in capturing at least some of the detections. This generalization holds true even if subspace X and subspace Y are not 37 orthogonal; there could always be some detections on the unit sphere that can be better identified by subspace X-Y. In order to test if the linear combination effect can be significant in real data, two events are chosen as training events, one each from pits A and B (Figure 4.2), and 12 target events located in pit D, which is situated between pits A and B. Subspace detectors are then created and run from the two training events, with a very liberal detection threshold, over the 12 days in which the pit D blasts occurred and were recorded by station M17A. Any detections with a detection statistic below the minimum value required to identify all 12 of the pit D events are then discarded. A similar procedure is conducted using the training events as templates in waveform correlation, and likewise, any detections below the correlation coefficient required to identify all 12 of the pit D events with at least one of the training events is discarded. The two-dimensional subspace detector returns 4,312 detections that do not correspond to any events in the blasting log, and the waveform correlators return 11,784 such detections (only counting detections made by both waveforms once). Both detection methods successfully identified 10 additional events in the log (from pit A and B) apart from the 12 blasts in pit D. The experiment is then repeated using the two largest events from pit D and forcing the detectors to identify the remaining 10. In this case, the subspace detector returns 1,740 unverified detections and the waveform correlators return 2,662. Therefore, in these experiments, subspace detection does have a clear advantage over waveform correlation. If maximizing detections is the priority, and false detections can be effectively handled, a higher dimension subspace detector with a lower detection threshold may be 38 preferable to a several lower dimension detectors (including waveform correlators) given two conditions are met: (1) the dimension of the higher order detector is much less than effective dimension of the embedding space, and (2) each added dimension captures some meaningful class of signal and not just noise. 39 Figure 4.1 Conceptualization of subspace detection and waveform correlation. See the text in Section 4.2. 40 Figure 4.2 Pit locations of the surface coal mine.41 CHAPTER 5 CONCLUSIONS A subspace detection technique is described and applied to continuous seismic data recorded by stations of the Transportable Array in order to identify seismicity from two distinct mining environments in the Western U.S. Independent ground truth information about the seismicity allows the assessment of the performance of the detector. In both cases, the subspace technique greatly increases the detection capabilities of the TA array. The detections and training events of each subspace are generally associated with specific areas of the mining operations, clustering fairly well in space with a few exceptions. The most effective method for mitigating false detections is to use a multiple station detector. Requiring overlaps in a detection's estimated origin time range is an effective way to associate detections together from multiple stations. The ratio of standard deviations performs well in estimating relative magnitudes (Equation 2.2). The standard deviation scheme outperforms the subspace equivalent of waveform scaling by accounting for signal energy not captured by the subspace. It is also shown that, in some cases, higher dimensional subspace detectors outperform lower dimension detectors (including waveform correlators) using the same training events. A tradeoff in detector performance between more distant stations, which tend to observe more similar waveforms, and closer stations, which tend to have higher SNR is apparent. 42 This tradeoff is an important design parameter that should be considered when establishing an array that will employ matched filter techniques to study seismicity from a specific area. The regional subspace detector identifies many events missed by the local array. Therefore, regional subspace detection may be useful in evaluating and augmenting catalogs created with local arrays and conventional network location procedures. 43 APPENDIX THEORY REVIEW A few concepts merit review for those not familiar with waveform similarity detection methods. Although not exhaustive or mathematically rigorous, the following brief introduction should be sufficient to understand the underlying principles of subspace detection. A.1 Seismic Waveforms The waveforms that are observed on a seismogram are a measure of ground motion (displacement, velocity, or acceleration) as a function of time. The data are usually digitally recorded with a fixed sampling interval (Stein and Wysession 2009). After removing the filtering effect of the seismic instrument, a given waveform (𝑠) is composed of two elements (Harris, 2006): 𝑠(𝑡)=Σ∫ℎ𝑖(𝜏)𝑔𝑖(𝑡−𝜏,𝜉,𝜒)𝑑𝜏𝑖 A.1 where ℎ𝑖 is the source time history and 𝑔𝑖 is the Green's function for a given source 𝑖. The source time history accounts for the physical phenomenon that releases seismic44 energy (slip on a fault, an explosion, etc.) and its progression through time. The Green's function accounts for the effects of propagation from the source location (𝜉) to the receiver (𝜒) including scattering, reflection, refraction, etc. If two waveforms are similar, the source time history and Green's function are also similar, implying that the two events occurred in similar locations and were caused by similar sources. Sources that vary only in energy amplitude but that occur in the same location tend to produce similar waveforms that differ only by a scaling factor over some frequency band (Richards et al., 2006 and references therein), which allows the use of waveform matching techniques as an effective detection method (Anstey, 1996). A.2 Cross Correlation In order to quantify amplitude-independent similarity between two digitized waveforms, vectors X and Y of length m and n where m ≤ n, the normalized cross correlation operator can be used, which is defined as: 𝐶[𝑖]=(𝑋[0:𝑚]−𝑋̅[0:𝑚])∗(𝑌[𝑖:𝑚+𝑖]−𝑌̅[𝑖:𝑚+𝑖])𝑚∗𝜎𝑋[0:𝑚]∗𝜎𝑌[𝑖:𝑚+𝑖] A.2 where brackets indicate indicial ranges and the cross correlation index, i, varies from 0 to n-m +1; sigma indicates the standard deviation; 𝑋̅ and 𝑌̅ indicate averages; and the star operator, ∗, represents multiplication if applied to scalars, or the dot product if applied to vectors. The value of the normalized cross correlation (𝐶) ranges from -1.0 to 1.0. A value of 1.0 indicates that the vectors are identical if one is multiplied by a positive 45 scalar, a value of 0 indicates the vectors are orthogonal, and a value of -1 indicates the vectors are identical if one is multiplied by a negative scalar. In addition to seismic event detection, cross correlation can be used to quantify similarity between a number of known seismic waveforms to form groups, or clusters, based on some similarity requirement (Houser et al., 2008), as well as for more exotic applications such as extracting empirical Green's functions from ambient noise (Bensen et al., 2007). A.3 Singular Value Decomposition Singular value decomposition (SVD) is a factorization of a matrix, similar to eigenvalue decomposition. Following the notation of Harris (2006), a matrix (𝑋)-the columns of which, in this application, are aligned waveforms-can be written in terms of an orthonormal basis (𝑊), a diagonal matrix of singular values (Σ), and a unitary matrix (𝑉): 𝑋= 𝑊Σ𝑉𝑇 A.3 where the superscript T denotes the transpose operation. Often the original matrix (𝑋) can be approximated by: 𝑋 ≈ 𝑊𝑑Σ𝑑𝑉𝑑𝑇 A.4 46 where the subscript d denotes truncation: only the d left-most columns of the matrices are retained. The value d, therefore, indicates how many orthonormal vectors are in the basis 𝑊𝑑, which is formally defined as the dimension of representation. A.4 Vector Projections A vector (V) of length n can be projected on an n-by-m basis (U): 𝑉𝑝 =𝑈𝑈𝑇𝑉 A.5 A useful measure of how well the basis U is able to reproduce V, or how well V "fits" into U, is a ratio of the projected vector energy over the original vector energy: 𝐶=𝑉𝑝𝑇𝑉𝑝𝑉𝑇𝑉 A.6 where C is referred to by Harris (2006) as the fractional energy capture if U was created using Equation 4 (𝑈 = 𝑊𝑑) and V is a column vector of 𝑋. The fraction of energy captured is one gauge of the validity of Equation 4 for a given 𝑋 and d. Alternatively, C may be used to test the hypothesis that a signal composed of a linear combination of the column vectors of 𝑈 is present in V. In the latter case, C is referred to as the detection statistic. 47 REFERENCES Anstey, N.A. (1966). Correlation techniques-A review, Can. J. Expl. Geo-phys., 2, 55-82. Astiz, L., Eakins, J. A., Martynov V. G., Cox, T. A., Tytell, J., Reyes, J. C., Newman, R. L., Gulsum, H. K., Mulder, T., White, M., Davis, G. A., Busby R. W., Hafner, K., Meyer J. C., and Vernon, F. L., 2014. The array network facility seismic bulletin: Products and an unbiased view of United States seismicity, Seismol. Res. Lett., 85, 576-593. Baisch, S., Ceranna, L., and Harjes, H. P., 2008. Earthquake cluster: What can we learn from waveform similarity?, Bull. Seism. Soc. Am., 98, 2806-2814. Barrett, S. A., and Beroza, G. C., 2014. An empirical approach to subspace detection, Seism. Res. Lett., 85, 594-600. Bensen, G. D., Ritzwoller, M. H., Barmin, M. P., Levshin, A. L., Lin, F., Moschetti, M. P., Shapiro, N. M., and Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int., 169, 1239-1260. Bonner, J. L., Pearson, D. C., and Blomberg, W. S., 2003. Azimuthal variation of short-period Rayleigh waves from cast blasts in Northern Arizona, Bull. Seism. Soc. Am., 93, 724-736. Darling, P. (Ed.), 2011. SME mining engineering handbook (Vol. 1). Society for Mining, Metallurgy and Exploration. Denver Colorado. Dodge, D., and Walter, W. R., 2015. Initial global seismic cross-correlation results: Implications for empirical signal detectors, Bull. Seism. Soc. Am., 105, 240-256. Ferretti, G., Massa, M., and Solarino, S., 2005. An improved method for the recognition of seismic families: Application to the Garfagnana-Lunigiana area, Italy, Bull. Seism. Soc. Am., 95, 1903-1915. Geller, R. J., and Mueller, C. S., 1980. Four similar earthquakes in central California, Geophys. Res. Lett., 7, 821-824.48 Gibbons, S. J., and Ringdal, F., 2005. The detection of rockbursts at the Barentsburg coal mine, Spitsbergen, using waveform correlation on SPITS array data, NORSAR Scientific Report: Semiannual Technical Summary No. 1 - 2005, NORSAR, Kjeller, Norway, pp. 35-48. Gibbons, S. J., and Ringdal, F., 2006. The detection of low magnitude seismic events using array-based waveform correlation, Geophys. J. Int., 165, 149-166. Gibbons, S. J., and Ringdal, F., 2012. Seismic monitoring of the North Korea nuclear test site using a multichannel correlation detector, Geoscience and Remote Sensing, IEEE Transactions on, 50, 1897-1909. Gibowicz, S. J., 2009. Seismicity induced by mining: Recent research, Adv. Geophys., 51, 1-53. Grêt, A., Snieder, R., and Özbay, U., 2006. Monitoring in situ stress changes in a mining environment with coda wave interferometry, Geophys. J. Int., 167, 504-508. Harris, D., 2006. Subspace Detectors: Theory, Lawrence Livermore Natl. Lab. Rep. UCRL-TR-222758, 46 pp., Lawrence Livermore National Laboratory, Livermore, California. Harris, D., and Dodge, D. A., 2011. An autonomous system for grouping events in a developing aftershock sequence, Bull. Seism. Soc. Am., 101, 763-774. Harris, D. B., Gibbons, S. J., Rodgers, A. J., and Pasyanos, M. E., 2012. Nuclear test ban treaty verification: Improving test ban monitoring with empirical and model-based signal processing, Signal Processing Magazine, IEEE, 29, 57-70. Harris, D., and Paik, T., 2006. Subspace detectors: Efficient implementation, Natl. Lab. Rep. UCRL-TR-222758, 36 pp., Lawrence Livermore National Laboratory, Livermore, California. Houser, C., Masters, G., Shearer, P., and Laske, G., 2008. Shear and compressional velocity models of the mantle from cluster analysis of long period waveforms, Geophys. J. Int., 174, 195-212. ISEE 2011. Blaster's Handbook 18th Edition. International Society of Explosives Engineers. Cleveland OH. Johnson, C., Bittenbinder, A., Bogart, B., Dietz, L., and Kohler, W., 1996. Earthworm: A flexible approach to seismic network processing. Iris newsletter. Joswig, M., and Schulte-Theis, H., 1993. Master-event correlations of weak local earthquakes by dynamic waveform matching, Geophys. J. Int., 113, 562-574. 49 Kubacki, T., Koper, K. D., Pankow, K. L., and McCarter, M. K., 2014. Changes in mining induced seismicity before and after the 2007 Crandall Canyon Mine collapse, J. Geophys. Res. Solid Earth, 119, 4876-4889. Maceira, M., Rowe, C. A., Beroza, G., and Anderson, D., 2010. Identification of low‐frequency earthquakes in non‐volcanic tremor using the subspace detector method, Geophys. Res. Lett., 37, L06303, doi:10.1029/2009GL041876. McLaughlin, K. L., Bonner, J. L., and Barker, T., 2004. Seismic source mechanisms for quarry blasts: modelling observed Rayleigh and Love wave patterns from a Texas quarry, Geophys. J. Int., 156, 79-93. Megies, T., Beyreuther, M., Barsch, R., Krischer, L., & Wassermann, J. 2011. ObsPy - What can it do for data centers and observatories? Ann. Geophys. 54, 47-58. Pankow, K. L., Moore, J. R., Hale, J. M., Koper, K. D., Kubacki, T., Whidden, K. M., and McCarter, M. K., 2014. Massive landslide at Utah copper mine generates wealth of geophysical data, GSA Today, 24, 4-9. Peng, Z., and Zhao, P., 2009. Migration of early aftershocks following the 2004 Parkfield earthquake, Nature Geoscience, 2, 877-881. Poupinet, G., Ellsworth, W. L., and Frechet, J., 1984. Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California, J. Geophys. Res., 89, 5719-5731. Richards, P. G., Waldhauser, F., Schaff, D., and Kim, W.-Y., 2006. The applicability of modern methods of earthquake location, Pure Appl. Geophys., 163, 351-372. Rubinstein, J. L., and Ellsworth, W. L., 2010. Precise estimation of repeating earthquake moment: Example from Parkfield, California, Bull. Seism. Soc. Am., 100, 1952-1961. Schaff, D. P., 2008. Semiempirical statistics of correlation-detector performance. Bull. Seism. Soc. Am., 98, 1495-1507. Schaff, D., 2010. Improvements to detection capability by cross correlation for similar events: A case study of the 1999 Xiuyan, China, sequence and synthetic sensitivity tests, Geophys. J. Int., 180, 829-846. Schaff, D. P., and Richards, P. G., 2004. Repeating seismic events in China, Science, 303, 1176-1178. Schaff, D. P., and Richards, P. G., 2014. Improvements in magnitude precision, using the statistics of relative amplitudes measured by cross correlation, Geophys. J. Inter., 197, 335-350. 50 Shearer, P. S., 2002. Parallel fault strands at 9-km depth resolved on the Imperial Fault, Southern California, Geophys. Res. Lett., 29, doi:10.1029/2002GL015302. Shelly, D. R., Beroza, G. C., and Ide, S., 2007. Non-volcanic tremor and low-frequency earthquake swarms, Nature, 446, 305-307. Slinkard, M., Schaff, D., Mikhailova, N., Heck, S., Young, C., and Richards, P. G., 2014. Multistation validation of waveform correlation techniques as applied to broad regional monitoring, Bull. Seism. Soc. Am., 104, 2768-2781. Smith, A. T., 1993. Discrimination of explosions from simultaneous mining blasts, Bull. Seism. Soc. Am., 83, 160-179. Song, F., Warpinski, N. R., Toksöz, M. N., and Kuleli, H. S., 2014. Full‐waveform based microseismic event detection and signal enhancement: An application of the subspace approach, Geophysical Prospecting, 62, 1406-1431. Spottiswoode, S. M., and Milev, A. M., 1998. The use of waveform similarity to define planes of mining-induced seismic events. Tectonophysics, 289, 51-60. Stein, J. R., Pankow K. L., Koper, K. D., and McCarter, M. K., 2015. Discriminating mining induced seismicity from natural tectonic earthquakes in the Wasatch Plateau region of Central Utah, 34th Int. Con. on Ground Control in Mining. Stein, S., and Wysession, M. 2009. An introduction to seismology, earthquakes, and earth structure. John Wiley & Sons. Stump, B. W., Hayward, C. T., Hetzer, C., and Zhou R. M., 2001. Utilization of seismic and infrasound signals for characterizing mining explosions, Proceedings of the 23rd Seismic Research Review on Worldwide Monitoring of Nuclear Explosions, 2-5 October 2001, Jackson, Wyoming. Swanson, P., Steward, C., and Koontz, W., 2008. Monitoring coal mine seismicity with an automated wireless digital strong-motion network, 27th Int. Con. on Ground Control in Mining. Waldhauser F., and Ellsworth, W. L., 2000. A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, Bull. Seism. Soc. Am., 90, 1353-1368. Wesson, R. L., 1987. Modelling aftershock migration and afterslip of the San Juan Bautista, California, earthquake of October 3, 1972, Tectonophysics, 144, 215-229. Wüster, J., 1993. Discrimination of chemical explosions and earthquakes in central Europe-a case study, Bull. Seism. Soc. Am., 83, 1184-1212. 51 Zhang, J., Song, X., Li, Y., Richards, P. G., Sun, X., and Waldhauser, F. 2005. Inner core differential motion confirmed by earthquake waveform doublets. Science, 309, 5739, 1357-1360. Zhang, M., and Wen, L., 2015. Seismological evidence for a low-yield nuclear test on 12 May 2010 in North Korea, Seism. Res. Lett., 86, 138-1.