Improving regional tsunami early warning with seismic array techniques

Data collected via seismic and infrasound array deployments are leveraged in the geosciences to detect and characterize a myriad of natural and anthropogenic sources. These deployments consist of numerous sensors placed in a predetermined configuration to amplify signal strength and improve the efficacy of array processing techniques used to measure signal directionality and waveform coherence. High-fidelity feature extraction is often predicated on interstation distance as well as the frequency content and wavelength of an incident signal. Numerous array processing softwares analyze data in sequential frequency bands to obtain a more detailed characterization of a signal. However, current algorithms are limited in their ability to determine optimal array configuration for each band. We introduce an open-source Python code, called Cardinal, to process seismic and infrasound array data in discretized time–frequency space with the option of applying an adaptive array design to determine optimal subarray configuration for each frequency band. To reduce computational time, the array processing step can be run in parallel using multithreading. Moreover, the software has the capability to aggregate array processing results from different time–frequency pixels to produce separate sets of detections, or families, with added utility via the application of an adaptive semblance threshold, which aids in isolating signals-of-interest from coherent background noise. Upon appropriate configuration, Cardinal exhibits the potential to combine distinct seismic and infrasound phases into separate families.

Seismological studies are traditionally based on the observation of ground motions recorded by translation sensors. However, to assess a comprehensive description of any wavefield produced by seismic sources, the three components of rotation are as important as the three translation components. Due to improvements in instrumentation through the last two decades, the ground rotational motion is now an observable. Several recent publications show that 6 degrees of freedom (-dof) seismic station, recording the 3 translation and 3 rotation components of the ground motion, provides valuable information to locate events or to infer their source mechanism.Rotation motions can be recorded directly via dedicated sensors, or numerically derived via a dense array of seismometers. The formers primarily use technologies based on fiber optic gyroscopes, liquid-based systems, or mechanical principles. However, these broadband instruments do not have the high sensitivity required to detect "weak" ground movements. Conventional sensor arrays, on the other hand, use finite differences techniques to estimate reliable indirect rotation measurements, but these are limited at high frequencies (wavelength > 4 × array aperture).Since the end of 2025, a temporary experimental campaign is set at the Low Noise Underground Laboratory (LSBB) in Rustrel, France. Two rotational rate sensors, namely a BlueSeis 3A and a Eentec R3, are deployed in the underground galleries, jointly with five seismometers which complete the permanent seismic array. The dense seismic array co-located around the rotation sensors is used to compute array derived rotations (ADR), and validate the direct observations of rotational ground motions. Furthermore, six seismometers are deployed at the surface to form an array with an aperture of around fifty kilometers. It has been designed for long-term observations of regional and global seismicity and array processing analysis, such as beamforming techniques. This campaign allows to compare the performance of classical seismic array processing with innovative gradiometric approaches based on a single 6-dof station, focusing on detection, location, and characterization of seismic events within a common frequency band. A sensitivity study on rotation signals, in terms of instrumental conditions (array geometry, station quality) and processing parameters (signal duration, filtering), is carried out with the help of numerical full waveform modelling in order to quantify the uncertainty of the estimated source parameters. We present the benefit of such multi-component seismic wavefield recording, illustrated on several events of interest such as the recent (2025/07/29) Kamchatka event (Mw 8.8).

Earthquake early warning (EEW) systems that issue warnings prior to the arrival of strong shaking are essential in mitigating earthquake hazard. Currently operating EEW systems work on point‐source assumptions and are of limited effectiveness for large events, for which ignoring finite‐source effects result in magnitude underestimation. Here, we explore the concept of characterizing rupture dimensions in real time for EEW using small‐aperture seismic arrays located near active faults. Back tracing array waveforms allow estimation of the extent of the rupture front (as a proxy of the rupture size) and directivity in real time, providing complementary EEW capabilities for M>7 earthquakes to existing EEW systems. We implement it in a simulated real‐time environment and analyze the 2004 M 6 Parkfield, California, earthquake recordings by the U.S. Geological Survey Parkfield dense Seismograph ARray (UPSAR) array and the 2010 M 7.2 El Mayor–Cucapah earthquake recordings by strong‐motion sensors in San Diego, California. We find it important to correct for the bias in back azimuth induced by dipping structures beneath the UPSAR array, based on data from smaller events. Our estimated rupture length is 30% shorter than those inferred from other studies but still reasonable for EEW purposes. We attribute this difference to rupture directivity effects and the limited field of view of a single array. The accuracy of the approach may be improved with a network of arrays with overlapping fields of view. We demonstrate this by tracking the 2011 Tohoku earthquake rupture with two clusters of Hi‐net stations in Kyushu and northern Hokkaido. The obtained results are consistent with teleseismic back‐projection results and yield reasonable estimates of rupture length and directivity. Compared with other proposed finite‐fault EEW approaches, the array method is less affected by the coarseness of a Global Positioning System or seismic network and provides a high‐frequency characterization of the rupture that yields more suitable predictors of ground shaking for certain structures.

Array processing is an integral part of automatic seismic event detection pipelines for measuring apparent velocity and backazimuth of seismic arrivals. Both quantities are usually measured under the plane-wave assumption, and are essential to classify the phase type and to determine the direction toward the event epicenter. However, structural inhomogeneities can lead to deviations from the plane-wave model, which must be taken into account for phase classification and back-azimuth estimation. We suggest a combined classification and regression neural network, which we call ArrayNet, to determine the phase type and backazimuth directly from the arrival-time differences between all combinations of stations of a given seismic array without assuming a plane-wave model. ArrayNet is trained using regional P- and S-wave arrivals of over 30,000 seismic events from reviewed regional bulletins in northern Europe from the past three decades. ArrayNet models are generated and trained for each of the ARCES, FINES, and SPITS seismic arrays. We observe excellent performance for the seismic phase classification (up to 99% accuracy), and the derived back-azimuth residuals are significantly improved in comparison with traditional array processing results using the plane-wave assumption. The SPITS array in Svalbard exhibits particular issues when it comes to array processing in the form of high apparent seismic velocities and a multitude of frost quake signals inside the array, and we show how our new approach better handles these obstacles. Furthermore, we demonstrate the performance of ArrayNet on 20 months of continuous phase detections from the ARCES array and investigate the results for a selection of regional seismic events of interest. Our results demonstrate that automatic event detection at seismic arrays can be further enhanced using a machine learning approach that takes advantage of the unique array data recorded at these stations.

The statistical theory of signal detection and estimation has been applied to problems in large array seismology. Using this theory the structure of the optimum detector for a known signal and long observation time in additive Gaussian noise is derived. The array processing filter employed by the optimum detector is known as the maximum-likelihood filter. This filter also has the property that it provides a minimum-variance unbiased estimate for the input signal when it is not known, which is the same as the maximum-likelihood estimate of the signal if the noise is a multidimensional Gaussian process. A series of experiments was performed using data from the large aperture seismic array to determine the effectiveness of the maximum-likelihood method relative to simpler methods such as beam-forming. These results provide significant conclusions regarding the design and processing of data from large seismic arrays. The conventional and high-resolution estimation of the frequency-wavenumber spectrum of the background microseismic noise is also presented. The diffuse structure of this spectrum is shown to aid in explaining the relative performance of array processing methods.

Distributed or decentralized estimation of covariance, and distributed principal component analysis have been introduced and studied in the signal processing community in recent years, and applications in array processing have been indicated in some detail. Inspired by these, this paper provides a detailed development of several distributed algorithms for array processing. New distributed algorithms are proposed for DOA estimation methods like root-MUSIC, total least squares-ESPRIT, and FOCUSS. Other contributions include distributed design of the Capon beamformer from data, and distributed implementation of the spatial smoothing method for coherent sources. A distributed implementation of a recently proposed beamspace method called the convolutional beamspace (CBS) is also proposed. The proposed algorithms are fully distributed – an average consensus (AC) is used to avoid the need for a fusion center. The algorithms are based on a recently reported finite-time version of AC which converges to the exact solution in a finite number of iterations. Numerical examples are given throughout the paper to show the effectiveness of the proposed algorithms.

Since the decommission of the underground launching control room of the ground-based component of the French nuclear missile system, the whole installation has been turned into a cross-disciplinary underground laboratory. The LSBB is a unique low-noise underground laboratory because of its initial military conception and its location in the regional park of Luberon far from large cities, industry and heavy traffic. The deepest point is 500 m below the surface. At this depth a huge and non-conventional shielded cylindrical capsule is installed with no μ-metal, 1268 m3 in volume, with a residual electromagnetic noise lower than 2 fT Hz−1/2 above 10 Hz. As a result, fluctuations of the Earth's magnetic field under 10 Hz can be recorded at a very low-noise level with a low-Tc SQUID 3-D magnetometer. Taking advantage of the main gallery topology, a broad-band underground seismic array has been deployed since 2001. An analysis of data recorded simultaneously by the seismic underground array and by the magnetometer sensors during the Indian earthquake of 2001 January 26 is presented. Evidence of a magnetic field perturbation induced by the seismic waves at teleseismic distance (6250 km) is supported by a polarization analysis of seismic and magnetic signals. Spectral analysis shows specific frequency bands of perturbation related to physical processes such as ground water flow acceleration within the mountain structure.

Abstract. While flowing downhill, a snow avalanche radiates seismic waves in the ground and infrasonic waves in the atmosphere. Seismic energy is radiated by the dense basal layer flowing above the ground, while infrasound energy is likely radiated by the powder front. However, the mutual energy partitioning is not fully understood. We present infrasonic and seismic array data of a powder snow avalanche, which was released on 5 February 2016, in the Dischma valley above Davos, Switzerland. A five-element infrasound array, sensitive above 0.1 Hz, and a seven-element seismic array, sensitive above 4.5 Hz, were deployed at a short distance (<500 m) from each other and close (<1500 m) to the avalanche path. The avalanche dynamics were modelled by using RAMMS (rapid mass movement simulation) and characterized in terms of front velocity and flow height. The use of arrays rather than single sensors allowed us to increase the signal-to-noise ratio and to identify the event in terms of back-azimuth angle and apparent velocity of the recorded wave fields. Wave parameters, derived from array processing, were used to identify the avalanche path and highlight the areas, along the path, where seismic and infrasound energy radiation occurred. The analysis showed that seismic energy is radiated all along the avalanche path, from the initiation to the deposition area, while infrasound is radiated only from a limited sector, where the flow is accelerated and the powder cloud develops. The recorded seismic signal is characterized by scattered back-azimuth angle, suggesting that seismic energy is likely radiated by multiple sources acting at once. On the contrary, the infrasound signal is characterized by a clear variation of back-azimuth angle and apparent velocity. This indicates that infrasound energy radiation is dominated by a moving point source, likely consistent with the powder cloud. Thanks to such clear wave parameters, infrasound is revealed to be particularly efficient for avalanche detection and path identification. While the infrasound apparent velocity decreases as the flow moves downhill, the seismic apparent velocity is quite scattered but decreases to sound velocity during the phase of maximum infrasound radiation. This indicates an efficient process of infrasound to seismic energy transition, which, in our case, increases the recorded seismic amplitude by ∼20 %, at least in our frequency band of analysis. Such an effect can be accounted for when the avalanche magnitude is estimated from seismic amplitude. Presented results clearly indicate how the process of seismo-acoustic energy radiation by a powder avalanche is very complex and likely controlled by the powder cloud formation and dynamics, and the process is hence affected by the path geometry and snow characteristics.

3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses

Crustal structure of the central Cathaysia block in South China derived from receiver functions

We present results from a study conducted to characterize the P and S wave trains and codas for local earthquakes observed near the San Jacinto fault zone, southern California. The data were recorded by a small‐aperture, three‐component seismic array and analyzed using a statistical array processing algorithm that provides estimates of the propagation direction, apparent wavenumber, and polarization characteristics for plane wave arrivals observed at the array. The P codas for these local earthquake data are similar in that each is composed of predominantly P waves traveling in a direction subparallel to the path of the direct P arrivals. These P coda exhibit a character consistent with propagation in a waveguide as opposed to propagation in a medium with randomly distributed isotropic scattering sources. The S codas differ from the P codas in two dramatic ways: (1) the S codas are composed of predominantly S waves, and (2) the S coda waves arrive from more random/spatially uniform directions. The dramatic differences between the P and S codas provide compelling evidence that P and S coda are the product of distinctly different scattering mechanisms and therefore that the P and S codas provide information about different source characteristics and/or structural features. A third important observation concerning the S coda is that the direction‐of‐arrival spectra for the S codas differ from event to event considerably more so than do those of the P codas. This observation is at odds with the conventional wisdom and suggests that the character of the S coda reflects a potentially complex interaction among several factors including, but not limited to, the source depth, the source distance, and the source radiation pattern.

: In this subject, we have considered a number of possible approaches to resolving mixtures of propagating signals and noises over seismic arrays. The problem is important because there will often be interfering noise sources, extraneous signals, or mixtures of phases that will be organized as plane waves. Analyzing such mixtures using conventional beam-forming and other methods will produce distorted velocities and azimuths and may even give incorrect magnitudes. One interesting result of this study is that beam-forming, single-signal F-Statistics as well as typical high resolution multiple-signal estimators from engineering literature, such as the Multiple Signal Characteristic (MUSIC) algorithm all suffer from the above problems for some very common situations involving regional and teleseismic data. We have approached the problem by showing that the multiple-signal plane-wave model is essentially in the form of a nonlinear regression in the frequency domain so that a sequential approach to isolating component signals analogous to stepwise linear regression can be applied. This approach provides a sequence of tests showing the power contribution of each component of a mixture and yields estimators for velocities and azimuths and their uncertainties. In the first year, we used this regression model, coupled with a corrected form of Akaike's Information Criterion (AICC) to setfie on a final mixture. Applying the AICC as a model selection criterion yielded three signals known to have occurred on this long-period mixture and a decomposition of a regional earthquake into a known signal and a propagating noise.

Synthetic seismograms produced by the finite difference method are used to study the scattering of elastic and acoustic waves in two‐dimensional media with random spatial variations in seismic velocity. The results of this study provide important insights about the propagation of short‐period (< 1 s) seismic waves in the earth's crust and place significant constraints on the fluctuation spectrum of crustal heterogeneity on length scales from tens of kilometers to tens of meters. The synthetic seismograms are analyzed to determine the variation in travel times and waveforms across arrays of receivers. The apparent attenuation caused by scattering and the time decay and amplitude of the seismic coda are also quantified with the numerical simulations. Random media with Gaussian and exponential correlation functions are considered, as well as a self‐similar medium with equal variations in seismic velocity over a broad range of length scales. These media differ in the spectral falloff of their velocity fluctuations at wavelengths smaller than 2π times the correlation distance a. The synthetic seismograms demonstrate that a random medium with self‐similar velocity fluctuations at length scales less than about 50 km (a ≥ 10 km) can explain both travel time anomalies reported for teleseismic arrivals across large‐scale seismic arrays (e.g., LASA and NORSAR) and the presence of seismic coda at frequencies of 30 Hz and greater commonly observed in microearthquake waveforms. Media with Gaussian and exponential correlation functions in velocity do not account for both sets of observations for reasonable standard deviations in velocity (≤10%). The scattering attenuation (Q−1) observed in the simulations for Gaussian media is peaked at ka between 1 and 2, where k is the seismic wave number. The observed attenuation in exponential media increases with frequency for ka < 1 and remains about constant for 1 ≤ ka ≤ 5.6. At high frequencies (ka > 5), the self‐similar medium is characterized by a scattering Q that is constant with frequency, whereas theory predicts that the apparent Q in an exponential medium is proportional to frequency. These alternative models of crustal heterogeneity can thus be tested by improved measurements of the frequency dependence of crustal Q at frequencies greater than about 1 Hz, assuming that scattering is responsible for most of the attenuation at these frequencies. Measurements of the time decay of the synthetic coda waves clearly show that the single scattering model of coda decay is not appropriate in the presence of moderate amounts of scattering attenuation (scattering Q ≤ 200). In these cases, Q values derived from the coda decay rate using the single scattering theory do not correspond to the transmission Q of the medium. The cross correlation of synthetic waveforms observed for an array of receivers along the free surface is observed to be dependent on the correlation distance of the medium. The self‐similar random medium proposed here for the crust produces waveform variations at high frequencies (15–30 Hz) similar to those reported for actual small‐scale seismic arrays with apertures of hundreds of meters.

Continuous wavelet transform and linear time-frequency transforms are coefficients of continuous unitary group representations of the affine and Heisenberg groups. Many properties of these transforms that are important for wideband radar and sonar signal processing follow directly from group representation theory. These properties include volume invariance and variance of narrowband and wideband ambiguity functions and wavelet transform domain implementations of detects and signal estimators. For several radar, sonar, and array processing applications, the basic definition of wavelet and time-frequency representations must be generalized by using unitary representations of other groups and using reproducing kernel Hilbert space (RKHS) inner products in the definition of the linear transforms. The general definition then leads to weighted continuous wavelet transforms where the RKHS is determined by a nonstationary covariance function; generalized wideband ambiguity functions also follow from this general definition along with other important generalizations that arise in wideband array processing and model based signal processing in complex scattering and propagation media. This paper presents the generalized wavelet transform along with the weighted wavelet transforms. The classical narrowband and wideband ambiguity functions are then special cases. The application of generalized transform to multidimensional transforms for space-time processing are also presented, along with the application to conformal groups for detections and estimation of accelerating scatterers.© (1998) COPYRIGHT SPIE--The International Society for Optical Engineering. Downloading of the abstract is permitted for personal use only.

This paper addresses the problem of source localization and waveform estimation in array processing applications. We present two nonparametric user parameter free algorithms, namely the iterative adaptive approach (IAA) and the maximum likelihood based IAA (IAA-ML). Both IAA and IAA-ML can work with arbitrary array geometries and uncorrelated as well as coherent sources. We extend IAA and IAA-ML to yield sparse results by using the Bayesian information criterion (BIC). Further improvements in performance can be achieved by employing the last step of the parametric RELAX algorithm initialized by the estimates of IAA and IAA-ML with BIC. Our simulations demonstrate that IAA and IAA-ML show better performance than the most competitive sparse signal representation approaches in the literature: M-FOCUSS, M-SBL and lscr <sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">1</sub> -SVD.