Single Seismogram Research Articles

Measuring Fundamental and Higher Mode Surface Wave Dispersion on Mars From Seismic Waveforms

AbstractOne of the goals of the Interior Exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) mission is to constrain the interior structure of Mars. We present a hierarchical transdimensional Bayesian approach to extract phase velocity dispersion and interior shear‐wave velocity (VS) models from a single seismogram. This method was adapted to Mars from a technique recently developed for Earth (Xu & Beghein, 2019, https://doi.org/10.1093/gji/ggz133). Monte Carlo Markov Chains seek an ensemble of one dimensional (1‐D) VS models between a source and a receiver that can explain the observed waveform. The models obtained are used to calculate the phase velocities of fundamental and higher modes at selected periods, and a subsequent analysis is performed to assess which modes were reliably measured. An advantage of our approach is that it can also fit unknown data noise, which reduces the risk of overfitting the data. In addition, uncertainties in the source parameters can be propagated, yielding more accurate model parameter uncertainties. In this study, we first present our technique and discuss the challenges stemming from using a single station to characterize both structure and the source and from the absence of a Mars reference model. We then demonstrate the method feasibility using the Mars Structure Service blind test data and our own synthetic data, which included realistic noise levels based on the noise recorded by InSight.

Travel Time and Waveform Measurements of Global Multibounce Seismic Waves Using Virtual Station Seismogram Stacks

AbstractWe construct geographically localized bin stacks of waveforms, called virtual stations, to enhance signal‐to‐noise ratios (SNRs) for travel time and waveform measurements of multibounce S and ScS phases (S up to S6 and ScS up to ScS5), as well as direct S, ScS, and Sdiff, on tangential component data. Major arc S and ScS multibounce waves were also measured. Virtual station data are referenced to empirical wavelets constructed from direct S waves for each event. The virtual station approach is useful for low SNR data, bolstering wave path coverage in the southern hemisphere. Goodness of fit measurements between the adapted empirical wavelet and virtual station waveforms are documented, as well as SNRs, allowing for objective definition of travel time measurement quality. From a data set of 360 earthquakes and 8,407 seismographic stations, nearly 4 million records were utilized to construct 248,657 virtual station stacked seismograms, which were compared to best‐fitting empirical wavelets. After human inspection of virtual station results, 8,871 travel time measurements were retained from 19 different minor and major arc seismic wave types. Higher multibounce data improve sampling of the southern hemisphere. From 188,003 single seismograms, 3,331 multibounce wave measurements were also made. Comparisons of single seismogram and virtual station stack measurements show a consistent bias: Virtual stack onset times are systematically early due to a broadening effect from stacking records with arrival time differences, which we correct for. The travel time and waveform measurements are publicly available.

Numerically Calculated 3D Space-Weighting Functions to Image Crustal Volcanic Structures Using Diffuse Coda Waves

Seismic coda measurements retrieve parameters linked to the physical characteristics of rock volumes illuminated by high frequency scattered waves. Space weighting functions (SWF) and kernels are different tools that model the spatial sensitivity of coda envelopes to scattering and absorption anomalies in these rock matrices, allowing coda-wave attenuation ( Q c o d a ) imaging. This note clarifies the difference between SWF and sensitivity kernels developed for coda wave imaging. It extends the SWF previously developed in 2D to the third dimension by using radiative transfer and the diffusion equation, based on the assumption that variations of Q c o d a depend solely on variations of the extinction length. When applied to active data (Deception Island, Antarctica), 3D SWF images strongly resemble 2D images, making this 3D extension redundant. On the other hand, diffusion does not efficiently model coda waveforms when using earthquake datasets spanning depths between 0 and 20 km, such as at Mount St. Helens volcano. In this setting, scattering attenuation and absorption suffer tradeoffs and cannot be separated by fitting a single seismogram energy envelope for SWF imaging. We propose that an approximate analytical 3D SWF, similar in shape to the common coda kernels used in literature, can still be used in a space weighted back-projection approach. While Q c o d a is not a physical parameter of the propagation medium, its spatially-dependent modeling allows improved reconstruction of crustal-scale tectonic and geological features. It is even more efficient as a velocity independent imaging tool for magma and fluid storage when applied to deep volcanism.

Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements

Resolution in seismic tomography intimately depends on data coverage, with different parts of seismograms sensitive to different parts of Earth's structure. In classical seismic tomography, the usable amount of data is often restricted because of approximations to the wave equation. 3-D numerical simulations of wave propagation provide new opportunities for increasing the amount of usable data in seismograms by choosing appropriate misfit functions which have direct control on Frechet derivatives. We propose new misfit functions for full waveform tomography based on instantaneous phase differences and envelope ratios between observed and synthetic seismograms. The aim is to extract as much information as possible from a single seismogram. Using the properties of the Hilbert transform, we separate phase and amplitude information in the time domain. To gain insight in the advantages and disadvantages of chosen misfit functions, we make qualitative comparisons of the corresponding finite-frequency adjoint sensitivity kernels with those from commonly used misfit functions based on cross-correlation traveltime, amplitude and waveform differences. The major advantages of our misfit functions are: (1) working in the Hilbert domain reduces non-linear behaviour of waveforms due to interaction of phase and amplitude information, and (2) we show with noise-free synthetic seismograms that it is possible to use a complete seismogram without losing information from low-amplitude phases. Complementary to instantaneous phase measurements, envelope measurements provide a way of using amplitude information of waveforms, which may also easily be extended to constrain anelastic properties. The properties of the kernels allow us to simplify the tomography problem by separating elastic and anelastic inversions. First indications are that the kernels remain well behaved in the presence of noise.

Automated multimode phase speed measurements for high-resolution regional-scale tomography: application to North America

SUMMARY A fully automated method for obtaining multimode phase speed measurements from a single seismogram has been developed and applied to a large data set of three-component longperiod seismograms in North America, constructing high-resolution phase speed maps on a continental scale. The method of our phase speed estimation is based on a fully non-linear waveform inversion by Yoshizawa & Kennett working with a global search method (the Neighbourhood Algorithm). The entire process of waveform fitting and the evaluation of the estimated phase speed have been fully automated employing several empirical quantitative measures, assessing the quality of waveform fit and the relative contributions of each mode in a chosen time window. The measured phase speed data undergo automatic screening for quality control, comprising the threshold evaluation of their reliability and outlier detection and removal. This new automated method has been applied to a large data set recorded at North American stations, including the latest transportable stations of USArray. Using long-period three-component seismograms recorded during the past eight years, we have successfully retrieved large numbers of regional surface wave paths, including over 20000 paths for the fundamental-mode Rayleigh waves over a wide range of frequencies, and over 10000 paths for the higher mode Rayleigh as well as the fundamental-mode Love waves. The consistentresultsoftheautomatedmeasurementproceduresuggestthatthemethodworkswell at regional distances, allowing us to perform a high-resolution mapping of multimode phase speeds in North America. The results of the automated waveform analysis also indicate some intrinsic limitations in the higher mode phase speed measurements from a single seismogram particularly in the short period range, mainly due to the overlapping of higher mode arrivals as well as coupling between mode branches. Despite such an innate difficulty in the higher mode dispersion measurements, the automated method allows us to construct reliable multimode phase speed maps. The current data set of ray paths is significantly biased towards the western half of North America, resulting in non-uniform horizontal resolution across the continent. This issue will be resolved by the future migration of the USArray stations to cover the central and eastern United States. The new automated method can be a useful tool for high-resolution mapping of regional 3-D shear wave structure including possible anisotropy.

A broad 660 km discontinuity beneath northeast China revealed by dense regional seismic networks in China

We examined the P wave velocity structure around the 660 km discontinuity at the tip of the subducting Pacific slab beneath northeastern China by forward modeling waveform triplication data. A total of 742 broadband seismograms were recorded by dense regional seismic networks in China from a deep earthquake that occurred near the border of east Russia and northeast China, providing an unprecedented density of ray coverage near the front edge of the subducting Pacific slab. Multiple P waves were observed on single seismograms in the distance range of ∼14°–29°. The P wave triplication shows the following two features: (1) the direct arrival traveling above the 660 km discontinuity (AB branch) extends as far as ∼29°, approximately 6° further than the prediction from velocity model of International Association of Seismology and the Physics of the earth's interior (iasp91); (2) the refracted wave propagating through the lower mantle (CD branch) appears at a distance a few degrees greater than that of the iasp91 synthetics. Forward waveform modeling suggests that they are best explained by a high‐velocity transition zone underlain by a ∼50–70 km thick 660 km discontinuity. The broadened discontinuity is likely caused by multiple phase transitions associated with the dissolutions of olivine and garnet components.

Crustal thickness in the central Andes from teleseismically recorded depth phase precursors

SUMMARY Although crustal thickness is a critical constraint for geodynamic models of Andean orogenesis, relatively few measurements exist from local seismic studies. In this study, we use reflections from the underside of the Moho as illuminated by intermediate depth earthquakes to provide new estimates of crustal thickness for the central Andes (16–34 ◦ S). These reflected signals (pmP) are identified as precursors to the depth phase pP when recorded at teleseismic distances (35–85 ◦ ). Although relatively small in amplitude, the pmP phase is often clear even on single seismograms. Less obvious pmP phases were enhanced by stacking traces from arrays. This method was most effective for events of M w > 6 and depth > 100 km. Crustal thickness determined in this study ranges from 59 to 70 km in southern Peru, from 49 to 80 km across the Puna-Altiplano and from 50 to 60 km above the Pampean flat slab. The lack of such phases for certain events may be evidence of heterogeneity at the Moho, and some precursors appear to correspond to intracrustal (magma?) and lithospheric mantle discontinuities. The results obtained using the pmP technique combined with those from other studies show the crustal thicknesses under the Central Andes are more variable than expected from purely isostatic considerations. We infer that one or a combination of geodynamic processes involving lower crustal flow, lower crustal delamination, and mantle–crustal lithospheric coupling can explain the highly variable topography on the Moho beneath the Central Andes.

Control of tall building vibrations by sealed tuned liquid column dampers

Sealed liquid column damper (TLCD) are tuned in the low-frequency range of application (say below 5 Hz) with respect to a selected natural frequency of the main structure by means of a geometrical transformation in analogy to the classical tuned mass damper (TMD). To improve the performance in MDOF buildings even further, the influence of neighboring modes is taken into account by fine-tuning in state space, rendering the optimal parameters modified. Final adjustments are easily made in the course of in situ testing. The passive sealed TLCD making use of the gas-spring effect and, since turbulent damping of assigned strength is present, is perfectly suited to counteract steady-state vibrations. Reduction of transient vibration peaks, observed within the initial period of the strong motion phase of earthquakes, requires active control, realized by controlled pressure input into the gas volume. The TLCD is commonly sufficiently damped and the vertical component of the seismic excitation does not cause parametric resonance. A sufficient condition based on the cut-off damping must be checked. Computer simulations and small scale experimental testing under time-harmonic-, single seismogram- and random- forcing confirmed the excellent performance and robustness of the TLCD. Optimal solutions for benchmark MDOF structures with multiple TLCD are presented when excited either by wind gusts or by earthquakes. A base-isolated structure is considered with the effective damping supplied by a TLCD in the basement. Copyright © 2005 John Wiley & Sons, Ltd.

Non-linear waveform inversion for surface waves with a neighbourhood algorithm-application to multimode dispersion measurements

SUMMARY A new technique for multimode dispersion measurement has been developed by employing fully non-linear waveform inversion for a path-specific 1-D profile using a neighbourhood algorithm (NA). One-dimensional models derived from waveform inversion are quite sensitive to the model parametrization and the reference model used to start the inversion. With different approaches to the parametrization of the shear wave speed profile, we can find models with significant differences in velocity variation with depth, which provide similar levels of fit to the observed waveforms. Although the models differ, the calculated phase dispersion for the first few modes of the surface waves are very close indeed. We therefore regard the 1-D models derived from the multimode waveform inversion as an implicit description of the path-specific dispersion for each of the modes. The inversion procedure using NA samples a substantial number of models and we select the 1-D wave speed model that achieves minimum misfit. Phase speeds for each mode branch are then calculated from the model parameters with the minimum misfit. Because we use the 1-D models as a summary of the average multimode dispersion along the path, we do not need to make assumptions concerning the nature of anisotropy. We can employ simple isotropic models and the method can be applied to both Rayleigh and Love waves independently. The new method is applied to sets of paths in the Australian region and provides stable measurements of multimode dispersion from a single seismogram. The multimode phase speeds measured using the NA inversion can be used to retrieve phase velocity maps as a function of frequency for higher-mode branches as well as the fundamental mode, which will be crucial constraints on 3-D Earth models.

Polarization of the Earth's teleseismic wavefield

The average polarization of the Earth's teleseismic wavefield is imaged by applying a vector stacking technique to ∼32 000 three-component broad-band seismograms. Polarization record sections capture the average radial-vertical particle motion in two narrow frequency bands: one above the microseism peak (0.5–2.0 Hz) and one below (0.03–0.07 Hz). Within these record sections, colour indicates polarization angle and brightness denotes particle motion linearity. These parameters are estimated using the time-domain method of Jurkevics (1988), enhanced to allow the simultaneous analysis of multiple sources. The low- and high-frequency record sections highlight different features in the teleseismic wavefield. The low-frequency stack shows P, S and their surface multiples as well as more exotic phases originating from combinations of reflections, refractions and conversions. The high-frequency stack detects fewer arrivals. However, the attenuation of the multiple S waves allows the detection of later-arriving core phases. Both the high- and low-frequency stacks illuminate regions containing consistent waveforms and warn of contaminating arrivals. The coherence of lower-mantle-turning P waves is indicated by the agreement between the observed and predicted polarization angles. In contrast, a large discrepancy between the observed and predicted low-frequency S angles evidences strong SPL contamination between 20° and 60°. Both polarization stacking and conventional amplitude stacking techniques allow the detection of seismic phases that are buried within the noise on single seismograms, but polarization stacking additionally yields polarization angle and particle motion linearity.

A mathematical theorem useful for the direct estimation of seismic source spectra

Abstract In this note, we prove a mathematical theorem for the direct estimation of the source spectra of earthquakes using a single seismogram. Considering a combination of Fourier spectra, the effects of the “source term” and the “attenuation term” from a single seismogram can be separated. Using signal-processing techniques, such an approach may be expressed in terms of the estimation of the Wigner distribution.

Inversion of sonic waveforms with unknown source and receiver functions

SUMMARY A seismogram that is several times the length of the source-receiver wavelet is windowed into two parts—these may overlap—to obtain two seismograms with approximately the same source function but different Green's functions. A similarly windowed synthetic seismogram gives two corresponding synthetic seismograms. The spectral product of the window 1 data with the window 2 synthetic is equal to the spectral product of the window 1 synthetic with the window 2 data only if the correct earth model is used to compute the synthetic. This partition principle is applied to well-log sonic waveform data from Ocean Drilling Project hole 806B, a slow formation, and used there to estimate Poisson's ratio from a single seismogram whose transmitter and receiver functions are unknown. A multichannel extension of the algorithm gives even better results. The effective borehole radius Rb, was included in the inversion procedure, because of waveform sensitivity to Rb. Inversion results for Rb agreed with the sonic caliper, but not the mechanical caliper; thus if Rb is not included in the inversion its value should be taken from the sonic caliper.

Variability of P660s phases as a consequence of topography of the 660 km discontinuity

Abstract Records from the broad-band NARS-NL array (deployed in the Netherlands, western Germany and Belgium) have been analyzed for P to S converted phases from the 660 km discontinuity (P660s). No individual seismogram from this array shows a clear P660s phase. Stacks of (up to 113) crosscorrelograms show P660s energy that is weak in comparison with similar stacks elsewhere. Very high amplitude observations of P660s in single seismograms from the same region have been made in a previous study ( Paulssen, 1985, Geophys. Res. Lett., 12: 709–712; Paulssen, 1988, J. Geophys. Res., 93: 10489–10500). Such variations may be due to relief on the 660 km discontinuity. In this paper we quantify the effects of different topographies on the time delays, amplitudes, and waveforms of the direct P wave as well as on the P660s converted wave. Synthetic P and P660s waveforms are calculated by applying a Kirchhoff-Helmholtz integration over the wave fields at the 660 km discontinuity interface. For topography on scales of a few hundreds of kilometers, the effects on the P wave are very small compared with those on the P660s phases. Time delays, focused and defocused amplitudes, and waveform distortion of P660s phases influenced by such topography diminish the coherence of P660s in (synthetic) stacks. Our results show that topography of the 660 km discontinuity with dominant scale lengths of a few hundreds of kilometers and variations in depth of ±15 km to ±25 km, is consistent with the pattern of P660s observations in single seismograms and stacks.

Incorporating Socorro magma body reflections into the earthquake location process

When network station spacing is large relative to event depths, constraining focal depths using direct arrivals is difficult. To help overcome this problem, we incorporated the SzS, SzP, and PzP Socorro magma body reflections into the earthquake location process. To find a velocity model appropriate for combining direct and reflected arrivals, we jointly inverted for 75 hypocenters and a one-dimensional velocity model. We used over 1000 direct arrivals and over 400 reflected arrivals, solving for P velocity ( Vp ) and Poisson's ratio ( v ) above the approximate base of the seismogenic zone, Vp and v below the base of the seismogenic zone, depth to the upper surface of the magma body, and station corrections. High-velocity and low-velocity initial models converged to the same solutions without damping. We found that Vp and v may both decrease slightly below about 10-km depth. Magma body depth is near 18.75 km. A close fit between observed and modeled reflected phase arrival times and an even lateral distribution of positive and negative reflected-phase residuals indicate that the magma body upper surface is nearly flat. In the Socorro area, combining direct and reflected phases reduces average focal depth error by more than a factor of 3 and average origin time error by a factor of 2 compared with using only direct phases. This is because reflected arrivals are more sensitive to changes in focal depth than direct arrivals, and combining upgoing (direct) and downgoing (reflected) times reduces the trade-off between focal depth and origin time. For this method to be successful, phases must be carefully identified. To help identify secondary phases, we used relative arrival time curves, plotting phase-minus- P times versus S -minus- P times. Picking reflection pairs (SzP and SzS, SzP and PzP, or SzS and PzP identified on a single seismogram) significantly reduces the chances of phase misidentification. The SzP-SzS and SzP-PzP pairs also help prevent trade-offs between reflector depth and focal depths.

Measurement ofQ−1 forSwaves in mudstone at Chikura, Japan: Comparison of incident and reflected phases in borehole seismograms

AbstractSeismograms from small events in the southern Kanto region have been recorded by a seismometer installed in a 732-m-deep borehole at Chikura observatory, Japan, where mudstone of early Pliocene age is found from the ground surface to the bottom of the borehole. Strong phases occurring 1.5 ∼ 1.6 sec after the S-wave arrival were interpreted as reflections at the ground surface. Ray path directions of the incident S waves were determined from the minimum principal axis of the three-dimensional trajectory ellipsoid, which represents the particle motion for a head wave of S phase. Transverse horizontal SH components were used to estimate the Q−1 value for S waves. Assuming the free surface acted as an 100% reflector for the SH waves and the incident SH phase as an input and the reflected phase as an output, we calculated system functions. The Q−1 value was measured from the transfer characteristics between the incident and the reflected phases in a single seismogram, thus no corrections were necessary for source or site (including instrument) effects. If we applied a power law model, the following relationship was obtained from the regression analyses of 20 events: log 10 ( Q − 1 ) = ( − 0.52 ± 0.48 ) ⋅ log 10 ( f ) − ( 1.28 ± 0.22 ) ( 1.0 < f < 5.0 Hz ) , where f is frequency in Hz and error values are the 95% confidence levels of the regression coefficients.

Pattern recognition for earthquake detection

The detector algorithms in use at date rely on negative decision logic: based on a model of the ambient noise process they detect all deviations, but many of them are false alarms. The principal alternative to this approach is pattern recognition, which tests on positive correlation with some known signal patterns. The Sonogram-detector realizes this scheme for single seismogram traces. Sonograms display spectral energy versus time. Suitably scaled, these images display only information which is significant to the detection process. Patterns of known earthquakes and noise signals are defined by means of these images. Event detection is performed by recognizing one of the patterns in the actual sonogram. The overall processing scheme is similar to the visual inspection of seismograms by the human observer. An off-line test installation for detecting local earthquakes proves the expected low false alarm rate, high timing accuracy and good detection probability of the Sonogram-detector.

A method of searching for radial velocity anomalies in the lower mantle using interfering signals

Radial velocity anomalies in the lower mantle that give rise to triplications in the travel-time curve for short-periodP waves will produce arrivals havingdT/dΔ values that differe by roughly 0.2–0.5 s/deg. The first two arrivals associated with such triplications will be separated by less than one second over a distance range of 4°–10° they may not, therefore, be separable visually on single seismograms, so that their presence can only be inferred from some measurable property that depends on their mutual interference. If there are lateral variations in the regions of anomalous velocity gradients, the interfering signals will also have different azimuths of arrival. Using two synthetic wavelets we have investigated the effect of interference on bothdT/dΔ and azimuth measurements at the Yellowknife Array. We found that if the interfering pulses have a dominant frequencyv, there is a range of time separations (0.30/v→0.55/v) over which the measureddT/dΔ and azimuth values may fluctuate by much more than the differences indT/dΔ and azimuth between the interfering signals. We have evaluated the following empirically defined functions for three different primary signals, and for three different relative amplitudes of the interfering signals:f τ(t), the drift function, which expresses how the measured slownesses,p, and azimuths, ϕ, differ from the slownesses and azimuths of the primary wavelets; Δf(τ), the range function, which describes the behaviour of the upper and lower bounds ofp and ϕ as a function of the difference in arrival times of the signals, τ andf ′ τ , studied the properties of these functions, and have outlined how these properties provide criteria based on the numerical and statistical characteristics of the arrival vectors, and on the waveform of the signal that will enable small radial velocity anomalies to be more clearly delineated.