Seismology [S]

S43B  MS:Exh Hall B   Thursday
Advances in Signal Processing Methods for Seismology I Posters
Presiding: P Ma, University of Illinois; A Ferris, Weston Geophysical Corporation

S43B-1299 

Automated Source Depth Estimation

* Junek, W N (wjunek@yahoo.com), Air Force Technical Applications Center, 1030 S. Hwy A1A, PAFB, FL 32925-3002, United States Nieves, J R (jorge@aftac.gov), Air Force Technical Applications Center, 1030 S. Hwy A1A, PAFB, FL 32925-3002, United States Kemerait, R C (kemerait@aftac.gov), Air Force Technical Applications Center, 1030 S. Hwy A1A, PAFB, FL 32925-3002, United States Woods, M T (mwoods@aftac.gov), Air Force Technical Applications Center, 1030 S. Hwy A1A, PAFB, FL 32925-3002, United States Creasey, J P (jon@aftac.gov), Air Force Technical Applications Center, 1030 S. Hwy A1A, PAFB, FL 32925-3002, United States

Source depth estimation is a key process in the discrimination of earthquakes and explosions for nuclear treaty monitoring. The lack of observable depth phases does not mean an event occurred at or near the surface. Shallow events can have closely spaced depth phases that are imperceptible to human analysts and regional events are often complicated by the simultaneous arrival of multiple phases, which makes the observation of depth phases even more problematic. We have developed an automated signal processing algorithm that estimates the depth of an event directly from the observed seismograms. It requires a series of observations of a single event by a network of seismic arrays and calculates site-specific depth phase delay times and ray parameters using cepstral processing and frequency-wavenumber analysis. These measurements are used to generate a series of depth estimates for each station in the observing network which are stacked to identify a consistent result. An iterative cepstral processing technique, and an adaptive, site specific, gamnitude threshold are used to reduce the high false alarm rate inherent to cepstral processing. A series of events having shallow free depth solutions, built from a mixture of regional and teleseismic observations, were used to test our algorithm.

S43B-1300 

Geophysical Wavelet Library: Applications of the Continuous Wavelet Transform to the Polarization and Dispersion Analysis of Signals

* Kulesh, M (mkulesh@math.uni-potsdam.de), Institute for Mathematics, University of Potsdam, Am Neuen Palais 10, Potsdam, 14469, Germany Holschneider, M (hols@math.uni-potsdam.de), Institute for Mathematics, University of Potsdam, Am Neuen Palais 10, Potsdam, 14469, Germany

Surface wave propagation in heterogeneous media can provide a valuable source of information about the subsurface structure and its elastic properties. The processing of experimental seismic data sets related to the surface waves is computationally expensive and requires sophisticated techniques in order to infer the physical properties and structure of the subsurface from the bulk of available information. Most of the previous studies related to these problems are based on Fourier analysis. However, the frequency- dependent measurements, or time-frequency analysis offer additional insight and performance in any applications where Fourier techniques have been used. This analysis consists of examining the variation of the frequency content of a signal with time and is particularly suitable in geophysical applications. The continuous wavelet transform gives a suitable general framework for solving these types of problems; this approach is powerful and elegant, but is not the only available for the practical applications. Other methods such as the Gabor transform, the S-transform or bilinear transforms can be used as well. The relative performance of time-frequency analysis from different approaches is primarily controlled by the frequency resolution capability. To perform the time-frequency analysis of digital seismic data, we propose in this contribution a new free software package developed by the authors and based on the continuous wavelet transform. This package allows to perform the direct and inverse continuous wavelet transform, 2C and 3C polarization analysis and filtering, modeling the dispersed and attenuated wave propagation in the time-frequency domain and optimization in signal and wavelet domains. The aim of these operations is to extract polarization properties, velocities and attenuation parameters from a seismogram. The novelty of this package is that we incorporate the continuous wavelet transform into the library where the kernel is the time-frequency polarization and dispersion analysis. This library has a wide range of potential applications that we illustrate with the analysis of synthetic and real data. http://users.math.uni- potsdam.de/~{}gwl

S43B-1301 

Extraction of P and S Wave Propagations and Measurement of Shear-Wave Splitting From Noise Data

* Miyazawa, M (linen@eqh.dpri.kyoto-u.ac.jp), Disaster Prevention Research Institute, Kyoto University, Uji, Kyoto, 611-0011, Japan Snieder, R (rsnieder@mines.edu), Center for Wave Phenomena, Colorado School of Mines, 1500 Illinois Street, Golden, CO 80401-1887, United States Venkataraman, A (anupama.venkataraman@exxonmobil.com), ExxonMobil Upstream Research Company, PO Box 2189, URC-GW3-917A, Houston, TX 77252-2189, United States

We extracted downward propagating P and S waves from industrial noise propagating down a borehole at Cold Lake, Alberta, Canada and also found shear-wave splitting from these measurements. The continuous seismic data are recorded at 8 sensors along a down-hole well during steam injections into 420 to 470 m-deep oil reservoir. We used cross-correlation of the waveforms observed at the bottom sensor (depth = 370m) and other 7 sensors (depths = 190, 210, 235, 262, 287, 307, 340 m) to extract the Green's function that accounts for wave propagation between sensors. Fast high-frequency and slow low-frequency signals propagating vertically from the surface to the bottom are found for the vertical and horizontal components of the wave motion, which are identified with P and S waves, respectively. The detected signals agree well with the travel time curves calculated for the seismic velocity model obtained by travel time tomography. The fastest S wave, polarized in the northeast-southwest direction, is about 1.9% faster than the slowest S wave polarized in the northwest- southeast direction. The direction of polarization of the fast S wave is rotated by about 20 degrees in clock-wise direction from the maximum principal stress axis as estimated from the regional stress field. This study demonstrates the useful application of seismic interferometry to field data to determine structural parameters.

S43B-1302 

Accurate earthquake location using automatically determined later phases

Gischig, V (gischigv@gmail.com), ETH Zurich, ETH-Hoenggerberg, Zurich, 8093, Switzerland * Leonard, M (Mark.Leonard@ga.gov.au), Geoscience Australia, Cnr Jerrabomberra Ave & Hindmarsh Drive, Canberra, 2609, Australia Gorbatov, A (Alexei.Gorbatov@ga.gov.au), Geoscience Australia, Cnr Jerrabomberra Ave & Hindmarsh Drive, Canberra, 2609, Australia

Within the development of the Australian Tsunami Warning System at Geoscience Australia attempts are being made to automatically detect the later phases pP, sP and S and to use them for more accurate locations. Therefore, the autoregressive method described in GSE/JAPAN/40 (1992) is included in the automatic earthquake location algorithm provided within the real-time system Antelope (developed by BRTT, Boulder, Colorado). This algorithm was applied to a database of waveforms recorded between January and February 2006. The arrival times determined with this method are compared with the results provided by the International Seismological Centre, ISC. For a mb=6.5 event it is shown, that it is possible to detect arrival times of later phases with this method. However, the onset times of the P-waves, mostly classified as emergent, are on average 0.2±0.4s later then reported by the ISC catalogue. For the S-arrivals the deviation is on average 1±1.8s.

S43B-1303 

The estimation of source position of long period seismic wave by Complex Meyer matching pursuit

* Matsubayashi, H (hiro_m@bosai.go.jp), National Research Institute for Earth Science and Disaster Prevention, 3-1 Ten'nodai, Tsukuba, 3050006, Japan Matsumoto, T (mtakumi@bosai.go.jp), National Research Institute for Earth Science and Disaster Prevention, 3-1 Ten'nodai, Tsukuba, 3050006, Japan

We have developed the method of long period source position estimation which uses Complex Meyer matching pursuit (CMMP). Our method can estimate the source position from seismic array data, and the estimation can avoid the influence of source mechanism. We generally analyze the seismic array data by semblance method. The method is useful for the detection of wave propagating direction and the wave source position estimation. But if you analyze array data by the method, the seismic wave has to be coherent wave and longer wavelength than array span, because the method calculates the index of seismic wave coherency. Therefore, the estimation has the effect of the non-isotropical earthquake mechanism. Our method also can estimate wave propagating direction and the wave source position. Furthermore, our method can avoid the influence of earthquake mechanisms. First, we use CMMP for seismic signal picking from array data. Then, we get the reciprocal number distribution of variance of the picking time difference, and estimate the largest number position as the source position. The CMMP picking has 3 merits. The picking can be automatic. Our method can avoid the influence of source mechanism, because Complex Meyer wavelet function varies as seismic signals. And then, the analysis time window of our method is defined by Complex Meyer wavelet function, automatically. We can expect that the analysis results by our method are more equality than semblance. We applied our method to estimation of the source position of the long period (16sec) tremor at Aso volcano. We used the Hi-net tiltmeters data. The distances between tiltmeters are 8 – 27km. We could get the estimated epicenter position of the tremor which was 2.8km southwest from the active crater of Aso volcano. And the depth of the estimated source was 1km under sea level. The large number distribution was on small region. The north-south length of the region was about 1km, and east-west length of the region was about 3km. We could not apply the semblance method to the estimation, because of incoherency of the array data. But we could use our method for the estimation. Therefore, we suggest that our method is a useful method for array data analysis.

S43B-1304 

Application of Seismic Interferometry to Natural Earthquake Records

* Torii, K (torii@earth.kumst.kyoto-u.ac.jp), Kyoto University Faculty of Engineering, Room 118, C-1 cluster, Katsura Campus Kyoto Daigaku, Kyoto, 615-8530, Japan Matsuoka, T (matsuoka@earth.kumst.kyoto-u.ac.jp), Kyoto University Faculty of Engineering, Room 109, C-1 cluster, Katsura Campus Kyoto Daigaku, Kyoto, 615-8530, Japan Aizawa, T (aizawa@suncoh.co.jp), Suncoh Consultants Co., Ltd., 1-8-9, Kameido Koto-ku, Tokyo, 136-8522, Japan

Recently, seismic interferometry has been one of the hottest topics in the exploration geophysics since this can be applied to reflection seismology. Seismic interferometry constructs Greenfs functions between arbitrary two points by taking cross-correlation of records observed at two locations. These Green's functions correspond to the wavefields as if an impulsive source was set at one location and seismic wave propagates from this source to the other receiver. Therefore, if we chose two surface receivers, we can reconstruct reflection seismic data. In case of using many receivers in a survey line, taking cross-correlation of all observed data at receivers generates pseudo shot-gather data for arbitrary locations. This technique does not need information of time 0 as long as all receivers measure wavefields synchronously. Therefore, there is no limitation with regard to the cause of the seismic vibrations. Natural earthquakes may be very good seismic sources for seismic interferometry. In our study, we adopt data provided by Hi-net system for applications of seismic interferometry. In 1995, gThe Great Hanshin Earthquakeh struck around Osaka and Kobe in Japan. After that, Japanese Government decided to construct high-density and high-sensitivity sensor network all over Japan in order to accumulate effective information of earthquakes and understand the mechanism of earthquakes. This seismic network is called gHi- net systemh in Japan. Hi-net system provides us much effective information, origin time, epicenter, depth, magnitude and waveform, etcc We used waveforms provided by Hi-net and generated pseudo shot-gather data by using seismic interferometry. Then, we applied these data to conventional reflection survey and analyzed underground structures in Japan.

S43B-1305 

Receiver Function Multiple Suppression and Separation

* Bostock, M G (bostock@eos.ubc.ca), Department of Earth and Ocean Sciences, The University of British Columbia, 6339 Stores Rd, Vancouver, BC V6T 1Z4, Canada

Interrogation of the continental lithospheric mantle using P-receiver functions is rendered difficult by the contaminating influence of free-surface multiples originating from the crust-mantle boundary. Recently, the use of S-receiver functions has been advanced as a better means of achieving this objective due to the fact that conversions and multiples arrive in distinct temporal windows. However, S-receiver functions suffer from a number of other disadvantages, notably, high levels of signal generated noise and significantly smaller data sets for a given recording period. A means of suppressing/separating free-surface multiples within P-receiver functions could provide a valuable tool for corroborating S-recever function results and might allow more general characterization of the continental mantle using data from portable experiments. We examine the use of the single-scattering Kirchhoff approximation for 1-D media to this end. By assuming a correlation between material properties (through e.g. constant Poisson's ratio, Birch's law), we may construct operators that accomplish approximate separation of direct conversions and free-surface multiples for a single receiver function. We demonstrate the application of this approach through a series of numerical examples and assess its efficacy when the underlying assumptions (i.e. known material property correlations, 1-D structure) are not fully honoured. Applications of the method to data from permanent stations of the Global Seismic Network and Canadian National Seismograph Network will be presented.

S43B-1306 

New Codes for Ambient Seismic Noise Analysis

* Duret, F (fduret@usgs.gov), Polytech Paris UPMC, Bat. Esclangon 4 place Jussieu Case Courrier 135, Paris, 75005, France * Duret, F (fduret@usgs.gov), US Geological Survey, 345 Middlefield Rd. MS 977, Menlo Park, CA 94025, United States Mooney, W D (mooney@usgs.gov), US Geological Survey, 345 Middlefield Rd. MS 977, Menlo Park, CA 94025, United States Detweiler, S (shane@usgs.gov), US Geological Survey, 345 Middlefield Rd. MS 977, Menlo Park, CA 94025, United States

In order to determine a velocity model of the crust, scientists generally use earthquakes recorded by seismic stations. However earthquakes do not occur continuously and most are too weak to be useful. When no event is recorded, a waveform is generally considered to be noise. This noise, however, is not useless and carries a wealth of information. Thus, ambient seismic noise analysis is an inverse method of investigating the Earth's interior. Until recently, this technique was quite difficult to apply, as it requires significant computing capacities. In early 2007, however, a team led by Gregory Benson and Mike Ritzwoller from UC Boulder published a paper describing a new method for extracting group and phase velocities from those waveforms. The analysis consisting of recovering Green functions between a pair of stations, is composed of four steps: 1) single station data preparation, 2) cross-correlation and stacking, 3) quality control and data selection and 4) dispersion measurements. At the USGS, we developed a set of ready-to-use computing codes for analyzing waveforms to run the ambient noise analysis of Benson et al. (2007). Our main contribution to the analysis technique was to fully automate the process. The computation codes were written in Fortran 90 and the automation scripts were written in Perl. Furthermore, some operations were run with SAC. Our choices of programming language offer an opportunity to adapt our codes to the major platforms. The codes were developed under Linux but are meant to be adapted to Mac OS X and Windows platforms. The codes have been tested on Southern California data and our results compare nicely with those from the UC Boulder team. Next, we plan to apply our codes to Indonesian data, so that we might take advantage of newly upgraded seismic stations in that region.

S43B-1307 

Trained automatic pickers used to extract the full-recorded seismicity of the ANCORP network

* Nippress, S E (nippress@liverpool.ac.uk), Univeristy of Liverpool, Department of Earth and Ocean Sciences, Jane Herdman Laboratories, Liverpool, L69 3GP, United Kingdom Rietbrock, A (ariet@liverpool.ac.uk), Univeristy of Liverpool, Department of Earth and Ocean Sciences, Jane Herdman Laboratories, Liverpool, L69 3GP, United Kingdom

The developments in data acquisition, data storage and data transmission technologies in the last decade allow seismologist to collect new data at an unprecedented rate. On large datasets the traditional manual picking requires significant amounts of both man-power and time with sometimes variable results. In order to process these large recorded datasets reliable and accurate automatic event detection and phase picking has become essential. The ANCORP temporary seismic network (34 stations) ran continuously between November 1996 to March 1997 in northern Chile and southern Bolivia. Due to the high seismicity rate in this region only the largest events during the observation period were analyzed (~859 events), thus using only ~25% of the actual recorded seismicity. To fully exploit the full-recorded seismicity in this dataset we use a number of trained automatic picking algorithms. Manual picks are used to train a STA/LTA, a Tp and a statistical automatic picking algorithm at each station and establish the optimum picking parameters that produce the greatest accuracy for each algorithm at each station. The optimum picking parameters for each station show regional variability and have to be adjusted individually to get the best performance. We find that the STA window length is larger in the north (0.5sec compared to 0.2sec in the south) and the LTA window length is shorter to the east (27sec compared to 30sec in the west). The Ds and a parameters for the Tp picker are larger in the west (Ds is 10-6 rather than 10-5and a is 10sec rather than 8-9sec), while the trigger level is constant across the network. Using the optimum picking parameters the automatic pickers are set running on the continuous ANCORP continuous dataset. We locate the events initially using HYPO71 before refining the locations in a 3D velocity model. Additionally, local magnitudes for all events based on automatic determined peak-to-peak values are calculated. We find that after using the trained automatic picking and location procedure the number of events located increases from 118 to 408 (a factor of 4 increase) in just a 2 week period. We locate a large number of new events at 60-130km depth range as well as new events in the overlying continental crust and close to the trench. This new picture of the local seismicity at intermediate depths during the period that the ANCORP network was operational will help to us to better understand the dynamic and complex systems that generate this observed seismicity.

S43B-1308 

Using Small-Aperture Arrays to Characterize the Far-Regional P-Wavefield in Central Asia

* Ferris, A (aferris@westongeophysical.com), Weston Geophysical Corp., 181 Bedford St. suite 1, Lexington, MA 02420, United States Reiter, D (delaine@westongeophysical.com), Weston Geophysical Corp., 181 Bedford St. suite 1, Lexington, MA 02420, United States

Accurate knowledge of waveform characteristics observed at far-regional and near-teleseismic distances is critically important for nuclear monitoring in regions in which observations are sparse. However, seismograms of far-regional earthquakes are often difficult to interpret because of complex wave interactions with upper mantle structure. At these distances (14-22°), seismic discontinuities near 220, 410 and 660 km depth, as well as other mantle heterogeneity, produce interference phenomena that result in phase triplication, amplitude modulations and travel-time anomalies. For example, in central Asia, travel-time bulletins exhibit large residuals (1 - 8 sec) for the direct P arrival from far-regional events. These residuals are likely caused by complicated multi-pathing along the upper mantle propagation path and frequent phase misidentification due to the complexity of the seismograms. To improve the usefulness of far-regional arrivals in seismic event location and discrimination, we have modified some well-known array-processing techniques to increase the accuracy of primary and early coda phase identification recorded by small aperture (2 - 3 km) borehole arrays. To date the high fidelity recordings provided by these small aperture arrays have been primarily utilized for accurate picks of primary phase onset times and local earthquake studies. However, we have found that adequate determinations of other arrival characteristics are possible by modifying two well- known τ-p filtering (or velocity spectra analysis) methods: Nth-root and phase-weighted stacking. We have applied our modified methods to vertical-component data recorded at distances between 14- 18° by 9-element small-aperture arrays in central Asia. Our results indicate that in many cases we can identify closely-spaced arrivals by their slowness values in a time window that includes the direct P arrival, depth phases and arrivals from upper- mantle discontinuities.

S43B-1309 

Demonstration of the "Point Seismic Array" Concept Using Co-Located Rotational and Translational Sensors

* Abbott, R E (reabbot@sandia.gov), Sandia National Laboratories, Energetic Experiments and Dynamic Materials Dept., P.O. Box 5800, Albuquerque, NM 87185-1168, United States Aldridge, D F (dfaldri@sandia.gov), Sandia National Laboratories, Geophysics Dept., P.O. Box 5800, Albuquerque, NM 87185- 0750, United States Hart, D (dhart@sandia.gov), Sandia National Laboratories, Ground Based Monitoring R & E Dept., P.O. Box 5800, Albuquerque, NM 87185-0404, United States

Spatially distributed arrays of multi-component sensors are commonly used to infer the speed, direction, and type (compressional or shear) of incident seismic waves. Recent developments in seismic wave propagation theory indicate that coincident (in space and time) measurements of both three-component translational and three- component rotational motions provide all the necessary information to calculate these parameters. Thus, the spatially extended array, with all its attendant deployment, operational, and post-acquisition data processing issues and problems, may be dispensed with. We present the results of theoretical, numerical, and field data acquisition experiments (the latter using an Eentec R-1 rotational seismometer and vibroseis truck energy source) to demonstrate the concept of the "point seismic array". Preliminary observations suggest that the back- azimuth to the vibroseis source may be inferred reasonably accurately via cross-correlating the measured horizontal acceleration and vertical torsion. Scaling these two measurements leads to an estimate of the wave speed. 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.

S43B-1310 

Detecting Noisy Events Using Waveform Cross-Correlation at Superarrays of Seismic Stations

* Von Seggern, D H (vonseg@seismo.unr.edu), Nevada Seismological Laboratory, U. Nevada, MS 174 1664 N. Virginia St., Reno, NV 89557, United States Tibuleac, I M (ileana@seismo.unr.edu

Cross-correlation using master events, followed by stacking of the correlation series, has been shown to dramatically improve detection thresholds of small-to-medium seismic arrays. With the goal of lowering the detection threshold, determining relative magnitudes or moments, and characterizing sources by empirical Green's functions, we extend the cross-correlation methodology to include "superarrays" of seismic stations. The superarray concept naturally brings further benefits over conventional arrays and single-stations due to the fact that many distances and azimuths can be sampled. This extension is straightforward given the ease with which regional or global data from various stations or arrays can be currently accessed and combined into a single database. We demonstrate the capability of superarrays to detect and analyze events which lie below the detection threshold. This is aided by applying an F-statistic detector to the superarray cross-correlation stack and its components. Our first example illustrates the use of a superarray consisting of the Southern Great Basin Digital Seismic Network, a small-aperture array (NVAR) in Mina, Nevada and the Earthscope Transportable Array to detect events in California-Nevada areas. In our second example, we use a combination of small-to-medium arrays and single stations to study the rupture of the great Sumatra earthquake of 26 December 2004 and to detect its early aftershocks. The location and times of "detected" events are confirmed using a frequency- wavenumber method at the small-to-medium arrays. We propose that ad hoc superarrays can be used in many studies where conventional approaches previously used only single arrays or groups of single stations. The availability of near-real-time data from many networks and of archived data from, for instance, IRIS makes possible the easy assembly of superarrays. Furthermore, the continued improvement of seismic data availability and the continued growth in the number of world-wide seismic sensors will increasingly make superarrays an attractive choice for many studies.

S43B-1311 

Improved depth estimation using integrated small-aperture array and network processing

* Tibuleac, I M (ileana@seismo.unr.edu), University of Nevada, Reno, Laxalt Mining Eng. Bldg., # 174, Reno, NV 89557, United States Anderson, J G (jga@seismo.unr.edu), University of Nevada, Reno, Laxalt Mining Eng. Bldg., # 174, Reno, NV 89557, United States Biasi, G P (glenn@seismo.unr.eu), University of Nevada, Reno, Laxalt Mining Eng. Bldg., # 174, Reno, NV 89557, United States Seggern, D v (vonseg@seismo.unr.edu), University of Nevada, Reno, Laxalt Mining Eng. Bldg., # 174, Reno, NV 89557, United States

We are testing a new approach to estimate the depth of earthquakes reliably and rapidly. Accurate estimates of earthquake depth are very important for nuclear monitoring, as well as for hazard assessment. The use of teleseismic recordings to determine this parameter is particularly important for sparsely monitored regions. Depth can be estimated by measuring the time separation between the primary arrival (P) and associated depth phases (pP, sP). Depth phases, however, are difficult to identify and thus they are used infrequently by institutes monitoring global seismicity such as USGS, IDC, and USNDC. As an example, less than 15% of the events located by USGS have associated depth phases. Using small or medium-aperture (< 25 km) array processing integrated with network processing, we can recognize secondary phases that are not visible on individual stations. Our approach works because, for stations in the vicinity of each array, the arrival-time difference between primary and secondary phases varies slowly with increasing epicentral distance. We estimate P-arrival parameters at small-aperture arrays using crosscorrelation in the time domain. From these parameters, we derive a time-variable set of weights. We beam the weighted envelopes of autocorrelated waveforms from nearby network stations to form the integrated network beam (INB). The resulting time series should have a central pick (P) and symmetrical side picks (pP), similar to autocorrelation of waveforms with ghost arrivals in exploration geophysics. To estimate the pP arrival time, we apply an F-detector to the INB and its components. Using a dataset of well-located events recorded at calibrated arrays and the surrounding networks, we determine whether this depth-phase identification methodology is widely applicable at local, regional and teleseismic distances. We also determine whether the methodology will work for real-time processing, and whether it will provide reliable depth estimates for earthquakes as shallow as 10 km.

S43B-1312 

Towards tomographic imaging with microtremor arrays: the application of Hilbert-Huang Transform

* Chen, Y (yoc03001@engr.uconn.edu), University of Connecticut, 261 Glenbrook Road, Storrs, CT 06268, Liu, L (lanbo@engr.uconn.edu), University of Connecticut, 261 Glenbrook Road, Storrs, CT 06268, Mehl, R (Robert.Mehl@uconn.edu), University of Connecticut, 261 Glenbrook Road, Storrs, CT 06268,

Using microtremor measurements to study subsurface structure is very attractive for both earthquake engineering and earth-scale seismology. Measurement of microtremors allows for is relatively easy and can be widely applied. Research in this area has gained great momentum in the last few years. Nevertheless, on the microtremor signal processing front, the conventional fast Fourier transform (FFT) remains the foremost tool for fundamental time-frequency analysis. We introduce a novel approach for microtremor analysis by employing the so-called ~{!0~}Hilbert-Huang transform (HHT)~{!1~} to extract more information from the microtremor array records. HHT is a combination of the process of empirical mode decomposition and the Hilbert transform to construct the Hilbert-Huang spectrum that contains the time-frequency information of the recorded signals. HHT does not require the signals under analysis to be linear and stationary. Unlike FFT, the HHT does not assume the time-domain signal has harmonic components over the entire frequency-domain; the energy only concentrates in a few intrinsic mode functions (and this is essentially true for many natural signals or noises). In this study, we applied HHT to the empirical Green~{!/~}s functions extracted from cross-correlation between microtremor station pairs and obtained the group and phase velocity dispersion information for the fundamental and higher mode Rayleigh waves. The inversion of shear-wave velocity profile based on these dispersion curves are in agreement with geological structure revealed by other methods with higher order accuracy. With a dense array it is possible to tomographically image the shear wave velocity structure at a given site.

S43B-1313 

Reconstructing impedance from noisy data using Bayesian inversion: preliminary results

* Sanguinetti, R (rrsangui@mail.uh.edu), University of Houston, 312 Science & Research Building I, Houston, TX 77204, United States

Acoustic impedance inversion is a technique for estimating acoustic impedance of subsurface layers from stacked seismic data. If we consider the earth as a series of discrete impedance layers, then each interface has associated with it a reflection coefficient. A seismic trace is just this series of reflection coefficients convolved with a seismic wavelet. Inversion attempts remove the seismic wavelet from the reflection series. Once we have the reflection coefficients it is relatively simple process to reconstruct the impedance earth layers. The task is the recovery of the acoustic impedance from a band-limited reflection seismogram. In terms of inverse theory, it assumed that impedance constraints are provided at different time levels. Also, noise in the trace is modeled by means of the usual Gaussian assumption. The uncertainties of the constraints are also assumed Gaussian. Thus, Bayes's rule is used to combine the data likelihood (by data I mean trace and constraints) with the prior probability of input reflectivity to construct the posteriori distribution of the output reflectivity. Thus, the maximum a posteriori (MAP) solution will be computed by minimizing the associated cost function. On the other hand, a reference model is used. For reaching it, it is necessary to perform the deconvolution of the desired reflectivity, from the observed data. This simple model must be modified to include the presence of random noise, where it emphasizes that r represents the primary reflectivity, i.e., the reflectivity of the layered earth assuming that multiple events have been attenuated and true amplitudes have been recovered (as well as possible) by an initial processing. This term is included into the cost function, and used as a reference model. Once the reference model, wavelet and impedance constraints, and the cost function are set, the inverse modeling is performed: the hyper-parameter is set, and the Cauchy prior is calculated, and added to the cost function. The minimization algorithm is based on the Conjugate Gradient (CG) extension for nonlinear functionals: the CG strategy can find the minimum of the cost function in a reasonable number of iterations. However, as well known, an incorrect selection of the hyper-parameter may yield a solution that is unreasonable. Technically speaking, if it is too small the data misfit is small, which indicates that over-fitting occurred. On the other hand, if it is too large, the regularization term will dominate the cost function, and the problem reduces to find the minimum. From the inversion, it is possible to see recovered models with (or without) using the reference model. However, these show some distortion around the spike, which could be telling us that the solution is pretty unstable. After that, an impedance estimate is reached: this solution is affected by the remaining noise which is seen as smooth distortions of the recovered reflectivity function in both sides of the spike model. So, the use of a reference model improves the solution. In our case, the use of this reference model helped to improve the performance of the inversion algorithm, giving us better results.

S43B-1314 

Time-Reversal Location of the 2004 M6.0 Parkfield Earthquake Using the Vertical Component of Seismic Data.

* Larmat, C S (carene@lanl.gov), Geophysics Group EES-11 Los Alamos National Laboratory of the University of California, MS D443, Los Alamos, NM 87545, United States Johnson, P (paj@lanl.gov), Geophysics Group EES-11 Los Alamos National Laboratory of the University of California, MS D443, Los Alamos, NM 87545, United States Huang, L (ljh@lanl.gov), Geophysics Group EES-11 Los Alamos National Laboratory of the University of California, MS D443, Los Alamos, NM 87545, United States Randall, G (grandall@lanl.gov), Geophysics Group EES-11 Los Alamos National Laboratory of the University of California, MS D443, Los Alamos, NM 87545, United States Patton, H (patton@lanl.gov), Geophysics Group EES-11 Los Alamos National Laboratory of the University of California, MS D443, Los Alamos, NM 87545, United States Montagner, J (jpm@ipgp.jussieu.fr), Institut de Physique du Globe de Paris, Case 89 4, place Jussieu, Paris, 75252-05, France

In this work we describe Time Reversal experiments applying seismic waves recorded from the 2004 M6.0 Parkfield Earthquake. The reverse seismic wavefield is created by time-reversing recorded seismograms and then injecting them from the seismograph locations into a whole entire Earth velocity model. The concept is identical to acoustic Time-Reversal Mirror laboratory experiments except the seismic data are numerically backpropagated through a velocity model (Fink, 1996; Ulrich et al, 2007). Data are backpropagated using the finite element code SPECFEM3D (Komatitsch et al, 2002), employing the velocity model s20rts (Ritsema et al, 2000). In this paper, we backpropagate only the vertical component of seismic data from about 100 broadband surface stations located worldwide (FDSN), using the period band of 23-120s. We use those only waveforms that are highly correlated with forward-propagated synthetics. The focusing quality depends upon the type of waves back- propagated; for the vertical displacement component the possible types include body waves, Rayleigh waves, or their combination. We show that Rayleigh waves, both real and artifact, dominate the reverse movie in all cases. They are created during rebroadcast of the time reverse signals, including body wave phases, because we use point-like-force sources for injection. The artifact waves, termed "ghosts" manifest as surface waves, do not correspond to real wave phases during the forward propagation. The surface ghost waves can significantly blur the focusing at the source. We find that the ghosts cannot be easily eliminated in the manner described by Tsogka&Papanicolaou (2002). It is necessary to understand how they are created in order to remove them during TRM studies, particularly when using only the body waves. For this moderate magnitude of earthquake we demonstrate the robustness of the TRM as an alternative location method despite the restriction to vertical component phases. One advantage of TRM location is that it does not rely on a prior picking of specific phases (Larmat et al, 2006). In future work will be conducted TRM backpropagation using the horizontal displacement components of seismic data as well as study the source complexity (double couples). Our ultimate goal is to determine whether or not Time Reversal offers information about the source that cannot be obtained from other methods, or that complements other methods.

S43B-1315 

Calculation of focal loci used for locating an earthquake in complex media

* Zhao, A (ahzhao123@yahoo.com), Institute of Geophysics, China Earthquake Administration, 5 Minzu Xueyuan Nanlu Street,Haidian District, Beijing, 100081, China Ding, Z (ding@cdsn.org.cn), Institute of Geophysics, China Earthquake Administration, 5 Minzu Xueyuan Nanlu Street,Haidian District, Beijing, 100081, China Wang, C (wangcy@cea-igp.ac.cn), Institute of Geophysics, China Earthquake Administration, 5 Minzu Xueyuan Nanlu Street,Haidian District, Beijing, 100081, China Sun, W (sunwg3s@yahoo.com.cn), Institute of Geophysics, China Earthquake Administration, 5 Minzu Xueyuan Nanlu Street,Haidian District, Beijing, 100081, China

Focal loci are often required in earthquake location. However, It is extremely difficult to analytically express them when the earthquake lies in a complex model. Therefore, the calculation of focal loci is usually limited to simple media. In this paper, we put forward a method for calculating focal loci in complex media by means of a minimum travel time tree algorithm for tracing rays. The focal locus is constrained with differential arrival time so that the problem of origin time is avoided. From all the model nodes, we select a small part with smaller absolute residuals between observed and calculated travel time differences (or double-differences) as representative points of the focal locus. The representative point with minimal double-difference is assigned as an initial point. Ray paths from the initial point to the other selected representative points in the double-difference field actually represent the focal locus, and are calculated with a minimum travel time tree algorithm. When the obtained focal locus is rather rough due to the excessive amount of selected representative points, it can be improved by removing some ray paths to the representative points that there are less rays go through. In addition, we modified the stop condition in the minimum travel time tree algorithm in order to reduce computational time. Our method is applied to calculating focal loci of an earthquake in a complex model in different cases including velocity perturbations and noisy arrival data. The calculation results show that the method is feasible.

S43B-1316 

Accuracy of the Discontinuous Galerkin Method for Computational Seismology

* Hermann, V (hermann@geophysik.uni-muenchen.de), Department of Earth and Environmental Sciences, Theresienstr. 41, Munich, 80333, Germany Kaeser, M (kaeser@geophysik.uni-muenchen.de), Department of Earth and Environmental Sciences, Theresienstr. 41, Munich, 80333, Germany

We present a detailed error analysis for the ADER-DG method on tetrahedral meshes. We consider several different aspects responsible for the accuracy of the scheme: the order of the method due to the degree of the approximation polynomials within each element, the spacial sampling, i.e. the mesh width, and the number of propagated wavelengths. These factors affect the deviation of the numerically computed signal from the reference signal. Furthermore, it is crucial to choose an error norm that describes the accuracy of a synthetic seismogram quantitatively and separates amplitude and phase misfits in the time as well as the frequency domain. The results of this accuracy study help to select an appropriate approximation order for a given mesh spacing, which typically is determined by the geometrical complexity of the problem. Especially with respect to p-adaptive large-scale simulations, where the mesh width might be strongly heterogenous and the approximation order of the scheme is allowed to vary locally, we provide guidelines for the choice of the applied approximation to ensure a desired accuracy. Finally, we apply the proposed method to a test problem suggested by the SCEC code validation project and compare the results with those of other well-established schemes.

S43B-1317 

ADER-DG Seismic Wave Propagation on Unstructured Quadrilateral Meshes

* Castro, C E (cristobal.castro@geophysik.uni-muenchen.de), Munich University, Theresienstr. 41, Munich, 80333, Germany de la Puente, J (jdelapuente@geophysik.uni-muenchen.de), Munich University, Theresienstr. 41, Munich, 80333, Germany Käser, M (martin.kaeser@geophysik.uni-muenchen.de), Munich University, Theresienstr. 41, Munich, 80333, Germany Dumbser, M (iagmidu@iag.uni-stuttgart.de), University of Trento, via Messiano 77, Trento, 38050, Italy

The ADER-DG method for seismic wave propagation has been developed in the recent years for unstructured triangular and tetrahedral meshes. In the present work we apply the same numerical approach to unstructured quadrilateral meshes instead. As the physical domains where we solve the seismic wave propagation equations are in general very complex, we are forced to use unstructured meshes with deformed elements. The geometrical flexibility of triangular and tetrahedral meshes are already well known; nevertheless we consider unstructured quadrilateral meshes as they are computationally less expensive when solving the governing equations and also enable easier comparisons with other existing numerical methods based on this mesh topology. Considering this approach we can reduce the computational time because we need less quadrilateral elements than triangular elements to discretize a particular physical domain considering the same edge lengths. Furthermore, the quadrilateral elements are bigger than the triangular elements and thus computations remain stable for larger time steps. Another advantage is that for quadrilateral meshes we can use a nodal basis based upon Gauss-Lobatto-Legendre integration points, instead of a modal basis, which further reduces the computational costs. Finally, applying the ADER-DG method on unstructured quadrilateral meshes we can compare the results and performance with the well-established Spectral Element Method utilizing exactly the same mesh. A further goal, currently under development, is to carry out simulations on hybrid meshes. The aim it to mesh the geometrically complex parts of a domain with triangles or tetrahedral while the simpler parts are meshed with the more efficient quadrilateral or hexahedral meshes, both using the same high-order ADER-DG method. As a consequence, for realistic setups, the computational costs could be optimized without compromising both the fidelity to the model and the accuracy of the numerical solution.

S43B-1318 

Seismic Sensor Characterization Using Three-Sensor Coherence Analysis

* Hart, D M (dhart@sandia.gov), Sandia National Lab Ground Based Monitoring R & E (Org. 5736), PO Box 5800 MS0404, Albuquerque, NM 87185-0404, United States Merchant, B J (bjmerch@sandia.gov), Sandia National Lab Ground Based Monitoring R & E (Org. 5736), PO Box 5800 MS0404, Albuquerque, NM 87185-0404, United States Chael, E P (epchael@sandia.gov), Sandia National Lab Ground Based Monitoring R & E (Org. 5736), PO Box 5800 MS0404, Albuquerque, NM 87185-0404, United States

Sleeman, van Wettum, and Trampert (2006) recently introduced a method for measuring the intrinsic noise spectra of seismic sensors which relies on determining the mutual signal coherence among three similar, collocated instruments. This in contrast to the standard two-channel coherence tests where it is assumed that the incoherent noise can be equally distributed between the two sensors, or lumped to one of them. The three- sensor coherence technique uses the available comparisons among three sensors to uniquely distribute total incoherent noise power among the individual sensors. This technique also allows us to compare the relative phase and amplitude response between the sensors under analysis. Due to the wealth of sensor characteristics obtained from the three sensor coherence we are looking to add this technique to the suite of analysis tools we use for evaluating sensors at Sandia's Facility for Acceptance, Calibration and Testing (FACT) site. In this poster, we briefly describe the method and demonstrate its effectiveness using synthetic signals with known amounts of coherent and incoherent power. Next we apply the procedure to characterize broadband STS2 low-gain seismometers set up in a vault at the USGS's Albuquerque Seismological Laboratory. We also looked at determining self-noise of a Q330HR digitizer with this technique by feeding a single seismometer component (vertical) into three channels of a common digitizer electronics board. These test configurations allowed us to investigate the use of seismic background as the input signal versus a controlled input signal source. This provided insight into linearity versus SNR of the tested sensors and digitizer. Overall the technique appears promising for providing self-noise, and response characteristics.

S43B-1319 

Hipocentral Inversion Using Spatial Proximity Constraint

Santana, F L (flaviolemos@ufrnet.br) Medeiros, W E (walter@dfte.ufrn.br) * do Nascimento, A F (aderson@dfte.ufrn.br)

We developed a hypocentral inversion methodology which incorporates the spatial proximity constraint among hypocentres. This constraint is suitable to approach this inverse problem because it is expected that these hypocentres define the fault surface. The usage of these constraints stabilises the inverse problem by introducing bias with geological meaning. As consequence the location of the events is improved, as we demonstrate with synthetic and real data. The inverse problem was formulated through an optimisation algorithm which does not require derivatives, allowing the use of robust data fitting as well as the spatial proximity between hypocentres. http://www.dfte.ufrn.br

S43B-1320 

Recent Advances in Seismic Wavefront Tracking Techniques and Their Applications

* Sambridge, M (malcolm@rses.anu.edu.au), Research School of Earth Sciences, Australian National University, Canberra, ACT 0200, Australia Rawlinson, N (nick@rses.anu.edu.au), Research School of Earth Sciences, Australian National University, Canberra, ACT 0200, Australia Hauser, J), Research School of Earth Sciences, Australian National University, Canberra, ACT 0200, Australia

In observational seismology, wavefront tracking techniques are becoming increasingly popular as a means of predicting two point traveltimes and their associated paths. Possible applications include reflection migration, earthquake relocation and seismic tomography at a wide variety of scales. Compared with traditional ray based techniques such as shooting and bending, wavefront tracking has the advantage of locating traveltimes between the source and every point in the medium; in many cases, improved efficiency and robustness; and greater potential for tracking multiple arrivals. In this presentation, two wavefront tracking techniques will be considered: the so-called Fast Marching Method (FMM), and a wavefront construction (WFC) scheme. Over the last several years, FMM has become a mature technique in seismology, with a number of improvements to the underlying theory and the release of software tools that allow it to be used in a variety of applications. At its core, FMM is a grid based solver that implicitly tracks a propagating wavefront by seeking finite difference solutions to the eikonal equation along an evolving narrow band. Recent developments include the use of source grid refinement to improve accuracy, the introduction of a multi-stage scheme to allow reflections and refractions to be tracked in layered media, and extension to spherical coordinates. Implementation of these ideas has led to a number of different applications, including teleseismic tomography, wide-angle reflection and refraction tomography, earthquake relocation, and ambient noise imaging using surface waves. The WFC scheme represents the wavefront surface as a set of points in 6-D phase space; these points are advanced in time using local initial value ray tracing in order to form a sequence of wavefront surfaces that fill the model volume. Surface refinement and simplification techniques inspired by recent developments in computer graphics are used to maintain a fixed density of nodes as the wavefront evolves. In addition to being computationally efficient and robust, the new WFC scheme can also be used to track multi-arrivals in complex media, and thus may lead to new developments in seismic tomography.

S43B-1321 

Seismic wave injection - a general boundary condition for the Alterman and Karal hybrid approach

Matyska, C (cm@karel.troja.mff.cuni.cz), Charles University, Faculty of Mathematics and Physics, Department of Geophysics, V Holesovickach 2, Prague, 18000, Czech Republic * Oprsal, I (io@egmdpri01.dpri.kyoto-u.ac.jp), Disaster Prevention Research Institute, DPRI, Kyoto University, Gokasho, Uji, Kyoto, 611- 0011, Japan Irikura, K (irikura@geor.or.jp), Aichi Institute of Technology, Disaster Prevention Research Center, 1247 Yachigusa, Hachikusacyo, Toyota, Aichi, 470-0392, Japan

Originated by Alterman and Karal (AK) in 1968 as a domain coupling realized by method-tailored technical algorithm, the generalized hybrid approach of wave injection is described by binding two sub-volumes treated per partes by arbitrary wave-propagation methods. The generalized AK two-step procedure possibly combines the source and path effects computed by one arbitrary method, and local site effects computed by another (arbitrary) method using the first method's wave field as input. Advantage of the approach arises from the fact that the connection between the methods keeps the formal wave-injection boundary perfectly permeable for any part of the wave field and can be applied to a variety of hybrid formulations. This hybrid approach leads to more effective modelling of combined source, path, and site effects (e.g., by multiple second step computations with varying structure, using single first step input), saving computer memory and time. The main innovation of the paper is generalization of the boundary condtion acting between two complementary sub- volumes of AK hybrid wave-propagation methods.

S43B-1322 

Analytical Computation of Effective Grid Parameters for the Finite-Difference Seismic Waveform Modeling With the PREM, IASP91, SP6, and AK135

* Toyokuni, G (toyokuni@geo.kyushu-u.ac.jp), Department of Earth and Planetary Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi- ku, Fukuoka, 812-8581, Japan Takenaka, H (takenaka@geo.kyushu-u.ac.jp), Department of Earth and Planetary Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi- ku, Fukuoka, 812-8581, Japan

We propose a method to obtain effective grid parameters for the finite-difference (FD) method with standard Earth models using analytical ways. In spite of the broad use of the heterogeneous FD formulation for seismic waveform modeling, accurate treatment of material discontinuities inside the grid cells has been a serious problem for many years. One possible way to solve this problem is to introduce effective grid elastic moduli and densities (effective parameters) calculated by the volume harmonic averaging of elastic moduli and volume arithmetic averaging of density in grid cells. This scheme enables us to put a material discontinuity into an arbitrary position in the spatial grids. Most of the methods used for synthetic seismogram calculation today receives the blessing of the standard Earth models, such as the PREM, IASP91, SP6, and AK135, represented as functions of normalized radius. For the FD computation of seismic waveform with such models, we first need accurate treatment of material discontinuities in radius. This study provides a numerical scheme for analytical calculations of the effective parameters for an arbitrary spatial grids in radial direction as to these major four standard Earth models making the best use of their functional features. This scheme can analytically obtain the integral volume averages through partial fraction decompositions (PFDs) and integral formulae. We have developed a FORTRAN subroutine to perform the computations, which is opened to utilization in a large variety of FD schemes ranging from 1-D to 3-D, with conventional- and staggered-grids. In the presentation, we show some numerical examples displaying the accuracy of the FD synthetics simulated with the analytical effective parameters.