{'fraction_non_alphanumeric': 0.028118087459588225, 'fraction_numerical': 0.017738642164369576, 'mean_word_length': 4.693630418987649, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Automatic Method for Determination of Lg Arrival Times Using Wavelet Transforms

The regional phase Lg is used to estimate location and magnitude for sources closer than 1500 km. The complexity of Lg waveforms makes it difficult to consistently determine Lg arrival time, thus affecting source location with a single station or array. This study tests an automatic method for timing Lg arrivals using wavelet transforms to decompose the Lg signal into its components localized both in time and scale. A Continuous Wavelet Transform (CWT) using a Daubechies order two (db2) wavelet is applied to 10 seconds of raw data, containing the start of Lg. Initial positioning of the window is obtained using the standard Lg travel time tables. The coefficients at scale 8 from the db2 decomposition are squared and the resulting time series is represented by an approximation of the 4'th level Discrete Wavelet Transform (DWT) using a Haar wavelet. A threshold detector is then applied to the resulting time series to determine the Lg arrival time. The method was tested using well located earthquakes (USGS) and explosions from known mines (mb less than 4.0), recorded on the vertical components at TXAR (Lajitas, Texas) and PDAR (Pinedale, Wyoming) arrays. The Lg arrival time was automatically picked with a standard deviation of less than 1.5 seconds (less than 10 km location error) for well known locations. Location errors are larger with the increase in distance and smaller with the increase in signal to noise ratio of events.


INTRODUCTION
Lg is one of the most used regional phases, both in location and discrimination, and dominant for many events up to 1500 km epicentral distance. Lg is made up of short period guided waves including scattering contributions from crustal heterogeneity, modeled as a sequence of multiply reflected postcritical shear waves trapped in the crustal waveguide (Xie and Lay, 1994), or as a superposition of higher mode surface waves (Herrmann and Kijko, 1983;Bollinger, 1979). The properties of the phase include a period of 1 to 7 seconds, a tendency to horizontal polarization, and a mean surface velocity of about 3.5 km/s (Herrin and Minton, 1960). The physical mechanism of Lg generation is still under study (Gupta et al.,1997;Patton and Taylor, 1995;Zhang and Lay, 1995;Xie and Lay, 1994) and there is no widely accepted automatic way of picking its arrival. Different types of seismic phase detectors, in time-domain and frequency-domain and using particle motion processing or pattern matching were previously developed (see Withers et al., 1998, for a good review). None of them have proved to be acceptable for all situations. In their paper, containing the study of an automated, near-real-time, waveform correlation event-detection and location system, Withers et al., (1998), use the idea of detecting transients in the frequency content of seismograms for one of the digital algorithms which they test in order to detect seismic phases for teleseismic events. Their solution was to use adaptive STA/LTA processing (variable length window function of the dominant frequency of the signal) and they were able to enhance a large variety of phase arrivals over a wide range of source, path, receiver and noise scenarios. However, the length of the windows used affected the precision of their estimates, even in the case of adaptive processing. Envelope processing, suitable for cases when original waveforms are very complex, was used by Husebye et al., (1998), in local event semi-automatic locations, for a frequency band between 3 and 6 Hz. Their procedure, while simple and suitable for real-time monitoring purposes, was also affected by the lack of resolution due to windowing. Wavelet transform methods were used for P and S phase identification in three -component seismograms by Anant and Dowla (1997). Assuming that strong features of the seismic signal appear in the wavelet coefficients across several scales, they used a multi-scalar representation to construct "locator" functions that identify the P and S arrivals and developed the first detector based on an alternative method of signal processing using wavelet decomposition of the signal.
We present an automatic method of picking Lg phases based on detecting transients in frequency content of seismograms and obtain the envelopes of the waveforms with a discrete wavelet decomposition method. The wavelet analysis is a windowing technique, with variable sized regions, that allows the use of long time intervals in instances for which we want precise low frequency information and shorter intervals where we want high frequency information. The method is based on the capacity of the wavelet transform for detecting changes in the nature of a process and to express frequency (scale) content changing in time, which is not possible with Fourier analysis or even windowed Fourier analysis (Misiti et al., 1996;Strang and Nguyen, 1996). We make the assumption that Lg, being composed of surface waves following the P and S coda, constitutes a change in process and is likely to be seen as a change of pattern in the time -scale representation of the continuous wavelet transform (CWT). We also assume that the Lg duration on the seismogram is of the order of at least 3 seconds, based on the fact that the formation of Lg requires a significant amount of energy propagating as shear waves trapped in the crustal wave -guide (see also Xie and Lay, 1994). The automatic detection method was tested on events recorded at two arrays: TXAR (Lajitas, Texas) and PDAR (Pinedale, Wyoming).
The location and station configuration of the arrays are presented in Figure 1.

THE DATA
The data set of small magnitude events ( m bLg < 4.0 ) consists of explosions from four mining districts and aftershocks of the April 14, 1995, mw = 5.8 Alpine (West Texas) earthquake, recorded at TXAR, and of a set of explosions from Black Thunder Mine (Powder River Basin, Wyoming), recorded at PDAR. The data from the vertical component of one of the elements TX08 or TX09 at TXAR and PD03 at PDAR were analyzed for each event. Locations of the earthquakes and of the mines, the mine status (multiple mines or single mine), the distance and backazimuth to the respective array, as well as the number of events used in this study are presented in Table 1. Locations relative to TXAR and PDAR of the events used to assess the Lg detector are presented in Figure 2.
Red squares are mine locations, the blue square represents the Alpine aftershocks locations.
The 16 Alpine aftershocks were located by USGS (Preliminary Determination of Epicenters Catalog) northeast of TXAR, at a distance of 108 km, backazimuth 19 degrees.
They were used because they are well located relative to the main shock.
Two clusters of mines were found in Mexico south and southeast of TXAR by Sorrells et al. (1997) The seven events from the Black Thunder mine, at 360 km distance from PDAR, backazimuth 73 degrees, were well documented blasts (Stump and Pearson, 1997) located in the same pit; one was a simultaneous explosion, the others were delay fired explosions. Bonner et al, (1997), observed that a relatively smaller yield in a simultaneous explosion from Tyrone produced a similar Lg magnitude to a larger delay fired explosion.
For this reason we did not use the magnitude to characterize the events, but the signal to noise ratio (SNR) instead. For this study, the SNR is defined as the ratio of the sum of the squares of the coefficients of the continuous wavelet transform (CWT) at scale 16 in a 5 seconds window after and immediately before the Lg arrival.

Wavelets and Fourier Analysis
Fourier analysis assumes that the spectral content of the data does not change with time within the analysis window. Thus frequencies present at one time are presumed to be present with the same intensities at all times. In fact, the time varying spectrum is undefined in Fourier analysis. However, seismic data is clearly time-varying in frequency, a fact that has been realized for some time. Previous researchers have addressed this problem by computing Fourier spectra with overlapping windows and piecing them together. Such an approach is not altogether satisfactory, is costly in computation time and difficult to optimize, although it is headed in the right direction -the direction of wavelet analysis. In wavelet analysis, unlike the Fourier transform which is global in frequency and independent of time, the wavelet transform is a localized transform in space (time) and scale (generalized frequency). Figure 3 is an example of both transforms applied to a signal containing the same frequencies superposed as in Figure 3 a, or successive as in Daubechies order 10 (db10) continuous wavelet transform ( see Misiti et al., 1996;Strang and Nguyen, 1996), in subplots 3a and respectively 3b, shows the frequency content changing in time, which is not possible to identify in the Fourier analysis. The y axis in subplots 3a and 3b represents the scale; low scales correspond to high frequencies, larger scales correspond to low frequencies.
A first approach to the problem of obtaining the frequency content of a process locally in time was a windowing technique applied to the signal, called the Short Time Fourier Transform (STFT) that maps a signal into a two-dimensional function of time and frequency (Chui, 1997;Kumar and Foufoula-Georgiou, 1996). However, the time/frequency window of the STFT is fixed and hence unsuitable for detecting signals with both high and/or low frequencies.
In essence, wavelet analysis is an extension of the STFT in such a way that the time/frequency window is no longer fixed, but can be varied. The use of long time intervals where we want more precise low frequency information and shorter intervals where we want high frequency information is accomplished through the introduction of a location parameter and a scale parameter as opposed to the Fourier method which has only a scale parameter. One of the two parameters is the translation or location parameter as in the windowed Fourier transform case. The other parameter is a dilation or scale parameter, the latter corresponding to the frequency ω in the Fourier case.
A major advantage of analyzing a signal with wavelets as the analyzing kernels is that it allows for the study of signal features locally with a detail matched to their scale, i.e.
broad features at a large scale and fine features on small scales. This property is especially useful for signals that are non-stationary, have short-lived transient components, have features at different scales or have singularities, discontinuities in higher derivatives and self-similarity (Kumar and Foufoula-Georgiou, 1996).

Basic wavelets theory
The essence of wavelet analysis is to expand a given function f(t) as a sum of 'elementary' functions called wavelets. Wavelets are functions that satisfy certain requirements. The name wavelet comes from the two requirements that they should integrate to zero, waving above and below the x-axis and that the function has to be well localized. A wavelet is a waveform of effectively limited duration that has an average value of zero. They tend to be irregular and asymmetric (see Strang and Nguyen, 1996).
The wavelets are themselves derived from a single function ψ , called the mother wavelet, by translations and dilations (Priestley, 1996). The mother wavelet is chosen so that it satisfies the following conditions: is the Fourier transform of ψ (t) and Ψ( ) 0 0 = . The mother wavelet is 'well localized' and decays 'quickly' to zero at a suitable rate (has compact support).
A classical mother wavelet is the Haar function defined by:  Figure 4, subplot 1. The Haar wavelet is not continuous, therefore the choice of the Haar basis for representing smooth functions, for example, is not natural and economic. However, when looking for a sudden change in the signal, the Haar wavelet is recommended (Misiti et al., 1996).
Given a mother wavelet, a doubly infinite sequence of wavelets can be constructed by applying varying degrees of translations and dilations to the mother wavelet.
so that the parameter a, usually restricted to positive values, represents the scale parameter and b the translation parameter. The normalizing factor a −1 / 2 is included so that This equality is always true for a > 0 . Figure 4, subplot 2 presents the compactly-supported wavelet db2 of Ingrid Daubechies (named Daubechies order two wavelet). The names of the Daubechies wavelets ( see Misiti et. al, 1996) are written dbN, where N is the order and db the ìsurnameî of the wavelet.
Selection of the wavelet shape is an important decision. Each may show a part of the reality with specificity, and each may reveal something that the others had concealed.
Choosing the wavelet is an advantage of the wavelet application and it allows for flexibility to adjust the analysis depending on the signal type, but there is no generally accepted procedure for doing so. Anant and Dowla (1997), for P and S arrivals, matched the wavelets with the arrival shape to help locate the arrival in the seismic signal. They also used a pseudo-wavelet constructed from the P arrival shape to locate the S arrival. The Lg arrival is too complicated for this procedure, so that we used a very simple wavelet, db2, composed of only one cycle, for our study. Also, for the best characterization of the beginning of the Lg arrival, the Haar wavelet (db1) was considered the most suitable. .

Wavelet analysis
As Fourier analysis consists of breaking up a signal into sine waves of various frequencies, wavelet analysis is the breaking up of a signal into shifted (translated) and scaled (dilated or compressed) versions of the mother wavelet, each multiplied by an appropriate coefficient. Any function f (t) ∈L 2 , the space of square integrable functions, can be expressed as a linear superposition of ψ a,b (t ) { }.
The process of Fourier analysis is represented by the Fourier transform: which is the sum over all time of the signal f(t) multiplied by a complex exponential. The Fourier coefficients F(ω), when multiplied by a sinusoid of appropriate frequency ω, yield the constituent sinusoidal components of the original signal.
Similarly, the Continuous Wavelet Transform (CWT) is defined as the sum over time of the signal multiplied by scaled, shifted versions of the wavelet function ψ : where C are the wavelet coefficients, functions of scale and position, which multiplied by the appropriated scaled and shifted wavelet yields the constituent wavelets of the original signal ( Misiti et al., 1996). The continuous wavelet transform can operate at every scale and the analyzing wavelet is shifted smoothly over the full domain of the analyzing function. This distinguishes it from the Discrete Wavelet Transform (DWT). A good review of theoretical aspects related to DWT was presented by Anant and Dowla, (1997) and Priestley, (1996).
The advantage of the DWT is that, using dyadic scales and positions (based on powers of two) the analysis is more efficient and not redundant. For fixed values a 0 and b 0 , .. , and for a 0 = 1/2 and b 0 = 1, the Daubechies discrete wavelet family is: This class of wavelets forms a complete orthonormal basis for the space of square integrable functions, so f(t) admits the representation:  Misiti et al., 1996;Strang and Nguyen, 1996) yields a fast wavelet transform. Two filters, high pass and low pass, are used to obtain respectively the detail (D) and the approximation (A) of the signal.
Downsampling ( throwing away every second data point) is introduced to keep the total number of data points constant. The decomposition process can be iterated, with successive approximations (A) being decomposed in turn, so that the signal is broken down into many lower resolution components (the wavelet decomposition tree). A disadvantage of DWT, from a seismological point of view, is the downsampling steps that result in the length of detail j+1, j≥1, j ∈ N being half the length of detail j (see also Anant and Dowla, 1997). As an example, for a signal with a sampling rate of 40 samples/second, detail 1 and approximation 1 will have 20 samples/sec, and detail 4 and approximation 4 will have 2.5 samples /second. When looking for the overall trend of the signal, equivalent to its behavior at low frequencies, the trend is more clear with each approximation, when j increases, and there is a trade-off between the resolution (sampling rate) and the smoothness of the approximation. The concept of "frequency" is related purely to the complex exponential function and has no precise meaning when applied to other families of functions, unless a more general concept of "frequency" can be defined (Priestley, 1996). In a loose sense, the wavelet parameters a, b, or m, n may be identified with generalized frequency and with time: when parameter a decreases, the "oscillations" become compressed in the time domain, i. e., they exhibit "high frequency" behavior, and as the parameter a increases, the "oscillations" exhibit ìlow frequencyî behavior.

The method
A plot of typical raw waveforms from all locations is presented in Figure 5. The first five (5) subplots represent events recorded at TXAR. Subplot 6 represents an event recorded at PDAR. The search window of 10 seconds (400 samples) is also represented and the amplitude in nanometres @ 1 second is counts x 0.0073 for TXAR and counts x 0.005 for PDAR.
A continuous wavelet transform is applied to the raw data, for a scale from 1 to 20, with step 1, using the db2 wavelet. We used the db2 CWT in order to maintain constant sample rate of the signal. The method is illustrated in Figures 6-9 which show the processing steps for two Tyrone events recorded at the station TX08. Figure 6 shows the October 11, 1996 simultaneous explosion and Figure 8 shows the September 9, 1997 delay fired explosion. In Figures 6 and 8 The coefficients at scale 16 were empirically chosen as the best to be analyzed for all the events. The range of frequency they correspond to is 0.9 to 2.1 Hz, centered on 1.5 Hz (see Chui, 1997) , which is a band with significant energy in the Lg Fourier spectrum for all the events in this study as well as for events from other regions. As an example, the corner frequency for events with magnitude less than 4.0 in Eastern Canada, for distances between 100 and 900 km, was calculated by Hasegawa, (1983) to be higher than 2 Hz.
Also, 1 sec period Lg waves were used by Nuttli (1986) for estimating seismic magnitudes because their anelastic attenuation is small in shield and other geologically stable regions.
Studies of the decay of spectral energy of Lg waves with distance for different frequencies and different regions generally show a less rapid decay in the range of 0.5 Hz to 2 Hz (Campillo et al., 1985;Herrmann et al., 1997). Subplot 3 of both Figures 6 and 8 presents the coefficients at scale 16 ( the sampling rate of the initial signal is preserved).

RESULTS AND DISCUSSION
The standard deviation of the Lg and P arrival time difference, expressed in km, for each location as a function of the distance to the respective array is presented in Figure 10. automatically picked with a standard deviation of less than 10 km for all the single mines.
The Mexico mines, considered to be multiple mines, have standard deviations greater than 10 km and less than 15 km.
The plot of the standard deviation as a function of Lg SNR threshold, in Figure 11, for all the single mine locations generally shows smaller errors with increase in SNR.
However for larger SNR fewer events were considered. There is no significant variation in the standard deviation with SNR for the Alpine aftershocks.
The raw data with the mean removed for the October 11, 1996 simultaneous explosion at Tyrone are presented in Figure 12, subplot 1. Subplot 2 shows the db2 CWT coefficients at scale 16 and the waveforms filtered between 0.9 and 2.1 Hz are presented in subplot 3. A zero phase, forward and reverse, 4'th order Butterworth filter was applied on the raw data in subplot 1 and, due to the very narrow frequency range, the resulting waveforms loose features. Also, the amplitude difference between the Lg signal and the ìnoiseî before Lg is much less than in the case of wavelet coefficients. We conclude that the wavelet transform is a more efficient way of filtering data in narrow frequency bands.
Comparing the shapes of the Lg envelopes for the simultaneously fired and for the delay fired explosion at Tyrone (Figures 7 and 9, subplots 3) it is seen that the Lg signal from the simultaneous explosion is impulsive, followed by an exponential decay, like a delta function convolved with an exponential, while the delay fired Lg is emergent, like a boxcar convolved with an exponential. The same difference in shape is observed for the approximations of Lg squared CWT coefficients for the events at Black Thunder Mine, not only at scale 16, but also for the raw data. We compare in Figure 13 two Black Thunder events, a simultaneous explosion and a delay fired explosion. The raw data with the mean removed for the August 24, 1995 simultaneous explosion is presented in subplot 1 and the raw data with the mean removed of the delay fired explosion on June 23, 1995 is presented in subplot 3. The Haar approximation at level 4 of the squared CWT coefficients at scale 16, in subplots 2 and 4, is different for each type of explosion; i.e., a sharp increase in Lg for the simultaneous explosion and a less rapid rise of Lg followed by a slower decay of amplitude for the delay fired blast. If the length of the boxcar could be considered proportional to the duration of the source, then it is possible that this type of analysis could be used as a discriminant between simultaneous explosions and delay fired explosions, provided that both occur in the same location. It could also be an additional tool for discrimination between simultaneous and delay fired explosions besides the study of long period regional surface waves, described by Stump and Pearson (1997), applicable for the large cast shots in same data set.

CONCLUSIONS
The automatic method for determination of Lg arrival times presented here results in a location error of about 10 km or less for events at epicentral distances from 360 km up to 700 km. Larger location errors at the Micare coal mines and the Minas de Hercules iron mines are probably caused by the existence of clusters of mines. A future study will focus on the statistical discrimination of individual mines in such a cluster.
The method described here allows extraction of meaningful signal attributes from single station vertical records using a set of conditions of which many are empirically established. The use of different wavelets and, specifically, different scales will be investigated in future studies. An automatic method to choose the threshold for the threshold detector as well as a better way of recognizing Lg than linear interpolation are under study.   Subplot 1. Contains the raw data, with the mean removed, recorded at the station TX08.

BIBLIOGRAPHY
The y axis represents amplitude in nanometres @ 1 second is counts x 0.0073. The x axis represents samples (40 samples/sec) for all subplots.

Subplot 2
A continuous wavelet transform is applied to the raw data, for a scale from 1 to 20, with step 1, using the db2 wavelet. Absolute values of the CWT coefficients are