S13B-1794
Application of Frequency-dependent Multi-Channel Wiener Filters to Event Detection in 2D Three-component Arrays
We present a frequency-dependent Multi-Channel Wiener Filtering technique with linear constraints. This filter removes coherent noise from seismic array data. We implement and test first an unconstrained version of this filter. This filter is generally far more effective than stacking in removing coherent noise, and ---unlike F-k-filtering --- does not depend critically on regular array spacing, or uniformity of the surface response of the individual sensors. The unconstrained version of the filter achieves maximal noise suppression but can result in signal distortion. Two methods are then introduced to achieve signal preservation, with an adjustable parameter to control the balance between signal preservation and noise suppression. We present versions of the filter for vertical components only, and for three-component data. We apply this technique to a dense array (95 three-component seismometers, within an area of 130 m×56 m) located near the village of Cavola, Italy. This array suffers from relatively high noise levels caused by industrial activities (a nearby ceramics factory). The Multi- Channel Wiener filter applied to this data set results in an improvement in the Signal-to-Noise Ratio by up to 22 dB (at 8 - 15 Hz), using 35 stations. The unconstrained version of the filter yielded the best improvement in the Signal-to-Noise Ratio, but the constrained filter is useful when waveform distortion is not acceptable.
S13B-1795
Gaussian-Beam S-wave Polarization-based Location Method for Hydraulic Fracturing Induced Seismicity
We propose a new method to locate hydraulic fracturing induced seismicity using three-component waveforms based on the Gaussian-beam S-wave polarization modified from Rentsch et al. (Geophysics, 2007) who used P-wave polarization. As a highly automated method, the method only requires a window around the S-wave peak amplitude and does not depend on arrival-time picking as used in standard location routines. Therefore, it has the potential to locate hydro-fracturing induced microseismic events which have clear S-waves but whose P-waves are difficult to be accurately picked or are immerged in background noise. This is a common situation in hydraulic fracturing due to the high amplitude ratio of S-wave to P-wave generally found in the recorded microseismic signals. The initial direction of wave propagation from a receiver is obtained from the eigenvector associated with the minimum eigenvalue by the S-wave polarization analysis. Along ray paths, the S-wave energy of all available receivers are back propagated by weighting a modified Gaussian-beam factor around the rays, and event locations correspond to the regions with the maximum summed energy. We have successfully applied the approach to two synthetic data sets with white noise under different recording networks and analyzed the sensitivity of the signal-to-noise ratio (SNR) and the percentage of S-wave splitting (PSS) on location errors for a homogeneous velocity model. From these tests, we conclude that our location method offers the successful location determination for S-wave SNR down to 2.7dB (which corresponds to P-wave SNR of -5dB) for an arbitrary receiver layout and 14.7dB (P-wave SNR of 7dB) for a linear receiver geometry comparing with P-wave SNR 6.8dB and 9.1 dB respectively from Rentsch et al, and the variation of PSS is up to 30% without a large influence on the location accuracy for both cases under each SNR limit. Furthermore, we also applied the technique to real data of a hydraulic fracturing experiment with a linear receiver array, and all of the locations obtained by the presented procedure are in good agreement with corresponding locations obtained by arrival-time-based location routines when S-wave SNR is higher than 20dB, which also approximately agrees with results from synthetic data.
S13B-1796
Hidden Semi-Markov Models and Their Application
In the framework of detection and classification of seismic signals there are several different approaches. Our choice for a more robust detection and classification algorithm is to adopt Hidden Markov Models (HMM), a technique showing major success in speech recognition. HMM provide a powerful tool to describe highly variable time series based on a double stochastic model and therefore allow for a broader class description than e.g. template based pattern matching techniques. Being a fully probabilistic model, HMM directly provide a confidence measure of an estimated classification. Furthermore and in contrast to classic artificial neuronal networks or support vector machines, HMM are incorporating the time dependence explicitly in the models thus providing a adequate representation of the seismic signal. As the majority of detection algorithms, HMM are not based on the time and amplitude dependent seismogram itself but on features estimated from the seismogram which characterize the different classes. Features, or in other words characteristic functions, are e.g. the sonogram bands, instantaneous frequency, instantaneous bandwidth or centroid time. In this study we apply continuous Hidden Semi-Markov Models (HSMM), an extension of continuous HMM. The duration probability of a HMM is an exponentially decaying function of the time, which is not a realistic representation of the duration of an earthquake. In contrast HSMM use Gaussians as duration probabilities, which results in an more adequate model. The HSMM detection and classification system is running online as an EARTHWORM module at the Bavarian Earthquake Service. Here the signals that are to be classified simply differ in epicentral distance. This makes it possible to easily decide whether a classification is correct or wrong and thus allows to better evaluate the advantages and disadvantages of the proposed algorithm. The evaluation is based on several month long continuous data and the results are additionally compared to the previously published discrete HMM, continuous HMM and a classic STA/LTA. The intermediate evaluation results are very promising.
S13B-1797
Wavelet Based Seismogram Denoising: Application to Receiver Function Estimation
Computation of Receiver Functions from observed seismograms involves the process of deconvolution which is highly sensitive to the presence of noise. Earthquake waveforms are non-stationary signals that contain colored Gaussian noise. The frequency content of the noise generally overlaps with the frequency content of the original signal and hence it is not possible to simply frequency filter the Gaussian noise. Techniques like the use of optimum filters are not suitable for seismic signals as seismic signals are non-stationary. Previous attempts to improve the signal-to-noise ratio in Receiver Functions involved stacking at the expense of path averaging of the crustal structure. Here we demonstrate that denoising of the seismic signal using the wavelet transform to sparsely represent the data, thresholding the amplitude and reconstructing the signal from remaining coefficients leads to a significant improvement in the signal to noise ratio. Hence it is useful in stabilizing the deconvolution process used to compute Receiver Functions. We have conducted tests on synthetic seismograms by adding different percentages of noise and trying to retrieve the original signal from the noisy traces. Cross-correlation coefficients between the denoised signal and the synthetic trace have been observed to improve over the cross-correlation coefficients between the noisy signal and the synthetic trace. It has been shown that this method of denoising seismograms does not introduce any noticeable distortions in the original signal [Galiana-Merino et al, 2003]. Receiver Function computed from the denoised seismograms showed a better correlation with the Receiver Function computed from synthetic seismograms than the Receiver Function computed from noisy seismograms. The signal-to-noise ratio was observed to have improved significantly. The improved stability of deconvolution due to denoising aided a more robust estimation of Receiver Functions from the retrieved signals. We show that this technique works well on real data without introducing any noticeable artifacts into the data and thus leads to a better Receiver Function estimation.
S13B-1798
Ground Truth Comparison of Sensor Orientation determined from Polarization Analysis of Large-eatthquake Body-was at USArray Seismic Stations
Polarization of teleseismic body-waves from large earthquakes (M > 7.3) recorded from April 2004 to August 2008 at USArray stations are used to estimate the sensor orientation. Three-component broadband recordings at over 500 sites are demeaned and band-pass filtered between 5 and 100 sec to determine the polarization of the arriving wave at all USArray sites. We compute the three-dimensional single value decomposition of say the arriving P-wave and compare the resulting orientation to that of the theoretical plane wave approximation (i.e. azimuth and incidence angle) to obtain the sensor orientation. We do this for body waves of events arriving from multiple directions to obtain an average sensor orientation at each station. In addition, to take advantage of the USArray deployment geometry, we determine an empirical polarization estimate for each event from the robust stack of body-waves recorded at sites located within 150 km of the reference station. This empirical estimate three dimensional single value decomposition is compared to that of the reference station to determine the orientation of the sensor at the reference station for multiple events. Given that since September 2007 ground truth sensor orientation is known at USArray sites because the use of an Octans gyroscope to a level of +/- 0.5 degree, it is possible to test the accuracy of the single station theoretical estimate and that of the nearest neighbor estimate for different body-waves. We present the results of the polarization estimates with body waves of large events that vary within +/-2 degrees for most stations when using the nearest neighbor estimate compared to those of the measured ground truth sensor orientation. The single station estimates variation is larger.
S13B-1799
Broad-Scale Applicability of Correlation Detectors to China and Parkfield, California, Seismicity
We investigate the potential use of correlation detectors on a broad scale to improve seismic monitoring and reduce magnitude detection thresholds. Previous work has indicated from semi-empirical analysis and a case study in Xiuyan, China, that an order of magnitude improvement is possible comparing a correlation detector for similar events with a standard STA/LTA detector. This unit reduction in detection threshold is achieved with acceptably low false alarm rates of about one per day. Semi-similar events due to less than perfect matches arising from location and mechanism differences or source complexities can provide useful detections. Synthetic tests on 78,028 focal mechanisms indicate that statistically significant detections are still triggered for strike, dip, and rake variations as large as 55 degrees. The correlation techniques were then applied on a larger scale to 5,000 events at Parkfield, California, and 19,000 events in and near China. We are attempting to see how broadly applicable correlation methods can be applied to different tectonic settings and for what percentage of the seismicity. 111 million correlations were performed on Lg-waves for the events in China at 363 stations. Final results indicate two thirds of the 19,000 events can be detected by cross correlation using this relatively sparse regional network. For Parkfield 82% of the events studied can be detected by cross correlation. Correlation detection is able to find additional events beyond what standard processing detects for China (70% increase) and for Parkfield (factor of 10 increase like Gutenberg-Richter predicts). Most event separation distances for events that correlate at Parkfield are less than 1 km. The distribution of magnitude differences for events that correlate at Parkfield is not distinguishable from the input magnitude distribution. Detection magnitude threshold reduction of about 1 unit holds for large scale application to the 19,000 events in China and 5,000 events in Parkfield with false alarm rates of a few percent.
S13B-1800
Sensor Technology Comparison Method and Case Study, Part II Time Domain Analysis
This study compares an accelerometer with geophone impulse response, then with deconvolution, to
differentiate sensors' performance in the time domain. Part I (ref. 3) of this analysis covered frequency
domain comparison only. Digital synthesis, which can closely resemble the whole geophysical data
processing steps, provided further quantifying evidence of performance differentiators by statistical
measures in terms of system time delay, correlation, and coherence. From a synthesized Ricker wavelet, the
simulation revealed salient features of the two types of sensors. It discusses the interaction of noise level and
deconvolution processing. There are a number of ways to evaluate sensor performance with mathematics.
Traditionally, researchers used sensor transfer function for this purpose. One of the shortcomings of transfer
function approach is that both amplitude and phase response are discussed separately but the impact in
real world, time domain, takes into consideration both aspects. In frequency domain, it is difficult to establish
an absolute scale (normalize in mathematical sense) for the total difference. Therefore it is logical to discuss
the issues in time domain directly. Impulse response reveals hardware in time domain, instead of frequency
domain as transfer function and reveals the intrinsic features of each sensor, making it easier to analyze and
understand. Seismic wavelet interval in survey is normally in the range of 10s to 100s milliseconds and
reflection wavelet is shorter. This requires a very quick sensor with rise/fall time less than ms without
overshot/undershot. Therefore, the sensor's temporal resolution needs to be improved by deconvolution
which makes the total filter as a δ(t), delta function. In this research, we propose to adopt other tools
such as lagged maximum correlation in time domain and coherence function in spectrum (coherence
discussion is omitted for part III). Lagged maximum correlation method is a statistical process which inherently
addresses noise which is not easily capable of transfer function method. The output is sensitive to sensor
design, signal, as well as noise. In fact, lagged max correlation fully incorporates the impact of both amplitude
and phase response in one number – time delay which is directly correlated with locating of a reservoir,
source, or a targeting structure. This study demonstrates that lagged cross correlation function is able to
detect system time shift in seismic applications. By using Ricker wavelet and convoluting with impulse
response of the geophone, the system time delay can be as large as 100 μs. This large amount of time
delay should be of concern. The advantage of this method is that it is based upon statistical methods as well
as versatile. In impulse response analysis, this study shows drastic dissimilarity between the two sensor
types, an accelerometer and a geophone. After applying typical data processing filters and deconvolution to
the geophone, both look fairly similar but the residue in geophone response needs to be carefully removed
by further data processing.
http://www.colibrys.com
S13B-1801
Implementation of Near Real-Time Methods Using Surface Waves to Determine Earthquake Source Characteristics at the National Earthquake Information Center
We discuss the implementation of two near real-time methods for determining earthquake source characteristics using long-period surface waves at the US Geological Survey's National Earthquake Information Center (USGS/NEIC). Long-period (60 to 300 s) seismic waveforms are well suited for this purpose because they can be well modeled using simple propagation models and are less sensitive to source complexity than short period (1 s) waves that are commonly used for earthquake monitoring. A prototype system for Surface wave Location and Association in Quasi Real time (SLAQR) that employs very long period (> 60 s) vertical-component surface waves has been implemented in test mode using data from the Global Seismographic Network (GSN). SLAQR continuously back-projects waveform envelopes on a global grid using surface wave dispersion relations. Preliminary results show that this method, as currently in operation at the NEIC, can consistently locate global earthquakes down to a magnitude 5.5. The magnitude determination, which is based on a simple empirical relationship, is generally accurate to within 0.2 magnitude units. More importantly, SLAQR can provide reliable locations and magnitudes for very large earthquakes, such as the 2004/2005 Sumatra events, within 30-40 minutes of their origin time. Furthermore, since the magnitude calculation is based on long period data, magnitudes for slow earthquakes such as some ridge and tsunami earthquakes are not underestimated as commonly occurs in shorter period analysis. Future development will focus on a reliable triggering algorithm for automated event detection and the continuous calculation of moment tensors and earthquake depths from the spectral amplitude and phase measurements already produced by the system. A fully automatic system to determine centroid moment tensors using three component surface waves with periods between 150-300 s is also running at the NEIC. Two versions are currently operational: one in a research/evaluation mode and another fully incorporated into the NEIC Hydra system. New improvements in the area of reliability assessment are currently being tested. Future work will investigate the finite fault information contained in the centroid time and centroid location parameters, how to incorporate a priori knowledge of the fault orientation, and the use of noise characteristics in the automatic selection of channels.
S13B-1802
Processing Aftershock Sequences Using Waveform Correlation
For most event monitoring systems, the objective is to keep up with the flow of incoming data, producing a bulletin with some modest, relatively constant, time delay after present time, often a period of a few hours or less. Because the association problem scales exponentially and not linearly with the number of detections, a dramatic increase in seismicity due to an aftershock sequence can easily cause the bulletin delay time to increase dramatically. In some cases, the production of a bulletin may cease altogether, until the automatic system can catch up. For a nuclear monitoring system, the implications of such a delay could be dire. Given the expected similarity between a mainshock and aftershocks, it has been proposed that waveform correlation may provide a powerful means to simultaneously increase the efficiency of processing aftershock sequences, while also lowering the detection threshold and improving the quality of the event solutions. However, many questions remain unanswered. What are the key parameters for achieving the best correlations between waveforms (window length, filtering, etc.), and are they sequence-dependent? What is the overall percentage of similar events in an aftershock sequence, i.e. what is the maximum level of efficiency that a waveform correlation could be expected to achieve? Finally, how does this percentage of events vary among sequences? Using data from the aftershock sequence for the December 26, 2004 Mw 9.1 Sumatra event, we investigate these issues by building and testing a prototype waveform correlation event detection system that automatically expands its library of known events as new signatures are indentified in the aftershock sequence (by traditional signal detection and event processing). Our system tests all incoming data against this dynamic library, thereby identify any similar events before traditional processing takes place. In the region surrounding the Sumatra event, the NEIC EDR contains 4997 events in the 9 months following the mainshock, and only 265 events during the same period for the previous year, so this sequence represents a formidable challenge for any automatic processing system. Preliminary results suggest that a waveform correlation-based system can detect on the order of 10% or more of the aftershocks for this event. Results published in the recent literature suggest that significantly larger proportions may be achievable for other aftershock sequences with smaller fault ruptures; we investigate and report encouraging results from one such sequence. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy's National Nuclear Security Administration under Contract DE-AC04- 94AL85000.
S13B-1803
Finite difference calculation of traveltime in rectangular topology mesh
Finite difference methods have been widely employed in solving the eikonal equation, to calculate first arrival time. Most previous researchers used regular grids, which did well in small scale problems, but the regular grids required denser meshes to sample the irregular interfaces, such as the Moho discontinuity, the 660km discontinuity, and it is difficult to solve large scale problems that the spherical surface effect must be taken into consideration. We propose a new finite difference method to solve the eikonal equation in rectangular topology grids. This method can accept rectangular topology mesh with slight distortion. And we also provided a reverse timing procedure to deal with head waves. After adopting appropriate mesh, we tested this method to calculate the first arrival times in a two layer model with an irregular interface, and the result showed it has similar accuracy with ten times grid spacing in regular mesh. Then we compared result obtained by our method with ray tracing¡¯s result in IASP91 model, and the traveltime curves at free surface are very close. If we add a rolling 660 km discontinuity, the traveltime curve shows an obvious perturbation while delta ranges from 24 degree to 34 degree, which indicates the effect caused by topography of 660 km discontinuity. This algorithm is also suitable to quickly calculate the first arrival travel time at local distances where a laterally varying Moho plays an important role or at upper mantle triplication distances where the 410km and 660 km discontinuities show substantial topographic variation. But our method may fail when the topography is too difficult to be defined with a mesh consisting of well-behaved quadrangles (aspect ratio too large or the angles too different from 90 degree).
S13B-1804
An Energy Rate Magnitude for Large Earthquakes
The ability to rapidly assess the approximate size of very large and destructive earthquakes is important for early hazard mitigation from both strong shaking and potential tsunami generation. Using a methodology to rapidly determine earthquake energy and duration using teleseismic high-frequency energy, we develop an adaptation to approximate the magnitude of a very large earthquake before the full duration of rupture can be measured at available teleseismic stations. We utilize available vertical component data to analyze the high-frequency energy growth between 0.5 and 2 Hz, minimizing the effect of later arrivals that are mostly attenuated in this range. Because events smaller than M~6.5 occur rapidly, this method is most adequate for larger events, whose rupture duration exceeds ~20 seconds. Using a catalog of about 200 large and great earthquakes we compare the high-frequency energy rate (· Ehf) to the total broad- band energy (· Ebb) to find a relationship for: Log(· Ehf)/Log(Ebb)≈ 0.85. Hence, combining this relation to the broad-band energy magnitude (Me) [Choy and Boatwright, 1995], yields a new high-frequency energy rate magnitude: M· E=⅔ log10(· Ehf)/0.85-2.9. Such an empirical approach can thus be used to obtain a reasonable assessment of an event magnitude from the initial estimate of energy growth, even before the arrival of the full direct-P rupture signal. For large shallow events thus far examined, the M· E predicts the ultimate Me to within ±0.2 units of M. For fast rupturing deep earthquakes M· E overpredicts, while for slow-rupturing tsunami earthquakes M· E underpredicts Me likely due to material strength changes at the source rupture. We will report on the utility of this method in both research mode, and in real-time scenarios when data availability is limited. Because the high-frequency energy is clearly discernable in real-time, this result suggests that the growth of energy can be used as a good initial indicator of the ultimate shaking.
S13B-1805
Locating point diffractors in layered media by spatial dynamics
We present a new approach to the problem of detecting point diffractors from active source surface seismic data. We formulate an optimization problem in the configuration space of possible collections of scatterers, seeking to find one collection that best matches an available shot-gather. We solve this problem by constructing a birth-and-death spatial dynamic, which converges to the optimal solution. This spatial random process tends to preserve diffractors that produce a response similar to the observed data, and it tends to remove ones that do not. This dynamic by design does not have a resolution limit typical of migration based techniques, which allows its use for sub-wavelength sensing.
S13B-1806
Time-reversal imaging of Earthquake and Seismic hum
The time-reversal technique is based upon spatial reciprocity and time invariance. This method was successfully applied in the past to acoustic waves in many fields such as sound waves in water or air, ultrasonic waves in human bodies, and electromagnetic waves in free space and recently to seismic waves in seismology. We present here, applications of time-reversal method in Seismology to synthetic and real tests, by using normal mode theory in the PREM model (Dziewonski and Anderson, 1981). We back-propagate 3 components of seismic data at very long period (T > 120s) (complete seismogram and one-bit seismogram). We show that the focusing is primarily dependent on the phase rather than the amplitude of seismogram. An excellent focusing at location and time of earthquake is usually obtained. Ten years ago a few groups reported existence of Earth's background free oscillations even on seismically quiet days (the "Hum")(Suda et al.,1998; Kobayashi and Nishida, 1998; Tanimoto, 1998). We started a systematic investigation of station located worldwide (FDSN) data during quiet periods of time. In this work we show that the excited modes are almost exclusively fundamental spheroidal modes and time-reversal experiment of seismic hum data (2-6 mHz) is attempted for the first time (only the vertical component of seismic data).
S13B-1807
Automatic P- and S-Phase Picker for High-Resolution Local Source Tomography: Application to the Alpine Region
Resolution and reliability of tomographic velocity models strongly depend on quality and consistency of available travel-time data. Arrival times routinely picked by network analysts on a day-to-day basis often yield a high level of noise due to mispicks and other inconsistencies, particularly in error assessment. Furthermore, tomographic studies at regional scales require merging of phase picks of several networks. Since a common quality assessment is usually not available for phase data provided by different networks, additional inconsistencies are introduced by the merging process. Considerable improvement in the quality of phase data can only be achieved by massive re-picking of seismograms. Considering the amount of data necessary for regional high-resolution tomography, algorithms combining accurate picking with an automated error assessment represent the best tool to derive large suitable data sets. In this work, we present algorithms for automatic P- and S-phase picking at local to regional distances including consistent picking error assessment. These algorithms are tested and applied to a waveform data set of the greater Alpine region. The MPX software is used to derive high-quality P-phase arrival times and the quality-attributed automatic picks are inverted for regional 1-D and 3-D P-wave models of the greater Alpine region. The comparison with a tomographic model based on standard routine phase data extracted from the ISC Bulletin illustrates effects on tomographic results due to consistency and reliability of our high- quality data set. In our proposed S-wave picking approach, we combine three commonly used phase detection and picking methods to a robust automated picking procedure. Information from the different techniques provides an 'in- situ' estimate of timing uncertainty and phase identification of the automatic S-phase pick. The average accuracy of automatic picks and their classification is comparable with manually picked reference picks, even if a higher number of S-phase picks is downgraded to lower quality classes by the automatic picker. Together with the MPX picking tool, our approach offers the possibility to generate sets of high-quality P- and S-phase data suitable for local and regional tomography.