An Interpretable and Domain-Informed Real-Time Hybrid Earthquake Early Warning for Ground Shaking Intensity Prediction

Engineering ›› 2025, Vol. 49 ›› Issue (6) :190 -204. DOI: 10.1016/j.eng.2025.03.009

Research

Article

Author information +

History +

Published Online
2025-03-18

PDF (3293KB)

Abstract

In the face of the unrelenting challenge posed by earthquakes—a natural hazard of unpredictable nature with a legacy of significant loss of life, destruction of infrastructure, and profound economic and social impacts—the scientific community has pursued advancements in earthquake early warning systems (EEWSs). These systems are vital for pre-emptive actions and decision-making that can save lives and safeguard critical infrastructure. This study proposes and validates a domain-informed deep learning-based EEWS called the hybrid earthquake early warning framework for estimating response spectra (HEWFERS), which represents a significant leap forward in the capabilities to predict ground shaking intensity in real-time, aligning with the United Nations’ disaster risk reduction goals. HEWFERS ingeniously integrates a domain-informed variational autoencoder for physics-based latent variable (LV) extraction, a feed-forward neural network for on-site prediction, and Gaussian process regression for spatial prediction. Adopting explainable artificial intelligence-based Shapley explanations further elucidates the predictive mechanisms, ensuring stakeholder-informed decisions. By conducting an extensive analysis of the proposed framework under a large database of approximately 14 000 recorded ground motions, this study offers insights into the potential of integrating machine learning with seismology to revolutionize earthquake preparedness and response, thus paving the way for a safer and more resilient future.

Graphical abstract

Keywords

Domain-informed neural networks / Physics-informed neural networks / Earthquake early warning / Variational autoencoder / Bayesian updating / Spatial regression / Interpretable artificial intelligence

Cite this article

Download citation ▾

Jawad Fayaz, Rodrigo Astroza, Sergio Ruiz. An Interpretable and Domain-Informed Real-Time Hybrid Earthquake Early Warning for Ground Shaking Intensity Prediction. Engineering, 2025, 49(6): 190-204 DOI:10.1016/j.eng.2025.03.009

登录浏览全文

4963

注册一个新账户忘记密码

1. Introduction

The unpredictable nature of earthquakes poses a formidable challenge to global safety, often resulting in significant loss of life, destruction of infrastructure, and profound economic impacts. Between 1990 and 2019, earthquakes caused over 1.3 million fatalities globally [1] and natural catastrophes caused 76 billion USD of global insured losses in 2020 out of which earthquakes accounted for 34% (directly or indirectly) [2]. Furthermore, earthquakes can result in significant downtime for affected areas, as buildings, infrastructure, and businesses may need to be repaired or rebuilt, and environmental impacts due to construction waste and debris [3], [4], [5], [6], [7], [8], [9], [10].

The quest to mitigate these devastating effects has led to the development of earthquake early warning systems (EEWSs), designed to provide crucial advance notice of seismic events [11], [12]. These systems hold the potential to significantly reduce the damage caused by earthquakes, enabling pre-emptive actions that can save lives and safeguard critical infrastructure [13]. Despite considerable advancements in technology and methodology, the continuous evolution of EEWSs is imperative to address the complex dynamics of seismic activity and the diverse needs of communities worldwide [14].

The United Nations (UN) has set several goals to attain community resiliency against natural disasters. Resiliency, in the context of natural disasters, refers to a community’s ability to prepare for, respond to, and recover from such events. The UN’s Sendai Framework for Disaster Risk Reduction 2015–2030 [15] underscores the critical importance of enhancing community resilience against natural disasters. Among them are earthquakes, which are catastrophic events capable of generating great destruction, loss of human life, and economic loss if the necessary safeguards are not implemented [16], [17], [18]. This comprehensive strategy advocates for the earthquake risk reduction in national development planning, emphasizing the need for robust and reliable EEWSs as a cornerstone of earthquake preparedness efforts helping to prevent the destructiveness of civil structures and critical infrastructure [11], [19].

Recent advancements in the field of seismology and information technology have paved the way for innovative approaches to earthquake detection and warning. Studies by Gleeson [11] and Chung [12] have laid the foundation by reviewing the evolution of EEWSs, stressing the significance of leveraging new technologies to enhance system efficiency and accuracy. Further research [13] on the potential effectiveness of EEWSs across different regions exemplifies the global applicability and critical importance of these systems in reducing earthquake-induced casualties and damage. Such systems can provide advanced warning to communities in the path of an earthquake, allowing them to take necessary actions to protect themselves and their property. This can include actions such as evacuating buildings, slowing down fast-moving trains, stopping gas, or shutting down critical infrastructure.

Conventional EEWSs developed in the last three decades have been primarily based on the inferred physics of the seismic ruptures and wave propagation theories [20], [21], [22], [23], [24]. With the recent advancements in computational resources and the availability of big data, data-driven models have become widely recognized as effective and applicable alternatives. Utilizing search engine techniques and refinements [25], [26] offered foundational methods that enhance the speed and accuracy of earthquake detection using computational tools. The incorporation of deep learning (DL) and artificial intelligence (AI) into EEWSs has opened new avenues for the real-time determination of earthquake characteristics such as magnitude (

M)

, focal mechanisms, ground motion intensity, and so forth, offering enhanced predictive capabilities and more accurate alerts [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37]. However, most of the existing data-driven EEWSs are not highly reliable in their predictions of ground-shaking intensity or provide minimal practical time for decision-making [38], [39], [40].

There are three main types of EEWSs: regional, on-site, and hybrid [28], [41]. Regional EEWS uses a network of seismic sensors distributed throughout a region to detect the initial waves of an earthquake and generally estimate the location and

M

of the event. This information is then relayed to a central processing facility, which determines the expected ground shaking, and issues alerts to users in affected areas. On-site EEWS uses sensors installed at individual locations, such as buildings or critical infrastructure, to detect earthquake shaking and issue alerts to occupants. These systems are particularly useful for protecting vulnerable and critical structures and equipment, such as hospitals, power plants, and bridges. Hybrid EEWS integrates both regional and on-site sensors to offer a more comprehensive method for earthquake detection and warning. The working principle of EEWSs may vary within each category. Some of them directly provide information on the expected ground shaking intensity measures (IMs), such as peak ground acceleration (

P G A

) [29], [32], [34], [36], [42]; others estimate earthquake source parameters, such as location and

M

. These parameters are subsequently employed in combination with pre-calibrated ground motion models (GMMs) to issue warnings [31], [35], [43], [44], [45], [46].

In the last three decades, studies focusing on seismic analysis have been exploring several critical parameters [47], [48], [49], [50], [51], [52], and among them, the spectral acceleration (

S_{a}

) for a period (

T

) collectively denoted as

S_{a} (T)

which represents the acceleration response spectrum has attracted significant attention. This is because it uniquely combines the ground motion waveform characteristics, such as amplitude and frequency content, with the dynamic behaviour of structural systems [53] through the utilization of an idealized single degree of freedom (SDOF). Over time, the utility of

S_{a} (T)

as an IM has expanded in structural and earthquake engineering applications, including ground motion simulation [54], [55], [56], ground motion selection [57], [58], and EEWSs [28], [36].

In an effort to enhance EEWSs, this study proposes a novel data-driven DL-based hybrid EEWS referred to as hybrid earthquake early warning framework for estimating response spectra (HEWFERS) which employs early primary waves (p-waves) information and site information to estimate the ground motion

S_{a} (T)

in real-time and provides a powerful tool for community resilience. HEWFERS involves three main components: ① domain-informed variational autoencoder (VAE) [59], used to obtain a physics-based two-dimensional latent variables (LVs) space (denoted as LV in vector form) representation of ground motion

S_{a} (T)

; ② feedforward neural network (FFNN) [60], employed to estimate the on-site LVs in real-time by utilizing vectors of site characteristics (SC) and IMs computed from the first 10 s of arriving ground motion waves (IM_10s); and ③ Gaussian process regression (GPR) [61] based spatial regression model used to obtain regional estimates of

S_{a} (T)

at the target sites. The VAE and FFNN are further illuminated through explainable artificial intelligence (XAI) to provide insights into the LVs. By offering accurate and timely warnings, HEWFERS not only contributes to saving lives and protecting infrastructure but also significantly reduces economic impacts, showcasing the indispensable role of technology-driven solutions in combating the seismic challenge. This aligns with the UN’s disaster risk reduction goals, showcasing the potential for technology-driven solutions to address the challenges posed by natural hazards, especially for earthquakes leading to large numbers of fatalities around the world [62].

2. Conceptualization of the proposed HEWFERS framework

The HEWFERS framework is an innovative EEWS designed to operate during an earthquake, delivering rapid and accurate estimates of ground motion intensity at specific locations both before and after the seismic waves arrive. The framework is divided into two phases.

Phase I (left-hand side of Fig. 1) activates immediately after the onset of an earthquake rupture when seismic waves reach the station closest to the rupture (first station). A pre-trained FFNN is used to estimate two LVs at the first station, which are denoted as

{L V}^{1}

. This estimation is based on the

SC

value recorded at first station (denoted as

{S C}^{1}

) and IM_10s recorded at the first station (denoted as

{I M}_{10 s}^{1}

). The estimated

{L V}^{1}

is then input into a pre-trained, domain-informed VAE decoder to estimate the on-site

S_{a} (T)

spectrum at the first station (i.e.,

S_{a}^{1} (T)

). Next, the framework uses GPR-based spatial regression model to estimate the

S_{a} (T)

spectrum at the nth target station (i.e.,

S_{a}^{n} (T)

), providing a regional early warning. For this, the estimated

{L V}^{1}

is combined with

{S C}^{1}

{S C}^{n}

(SC at the nth station),

{I M}_{10 s}^{1}

, and the distance between the first and nth station

{(d}_{n, 1}

). These combined inputs are used in GPR model to build spatial correlations and estimate the LV at the nth station (denoted as

L V^{n}

). The estimated

L V^{n}

is then subsequently fed into the VAE decoder to obtain a prior estimate of

S_{a}^{n} (T)

. This process can be repeated for different target stations throughout the region. Thus, the first phase of the framework provides the prior estimate of the

S_{a} (T)

spectrum at each station in the region before the arrival of seismic waves, offering a regional early warning.

Phase II of the framework enacts during the earthquake event once the seismic waves have reached the target nth station (right-hand side of Fig. 1). Analogous to the Phase I computations at the first station, the pre-trained FFNN and VAE decoder are utilized to calculate the on-site estimate of

{L V}^{n}

using

{S C}^{n}

and

{I M}_{10 s}^{n}

. This on-site estimate of

{L V}^{n}

is then integrated with the prior estimate of

{L V}^{n}

obtained in Phase I through Bayesian updating, resulting in the final posterior estimate of the

S_{a}^{n} (T)

spectrum at the nth station. This updated estimate considers the initially recorded ground motion data (i.e., first 10 s) and provides a more reliable estimation of the ground motion intensity at the target location. The framework is computationally efficient, taking less than three seconds on average to process the necessary calculations.

3. Ground motion database

A comprehensive database of unprocessed bi-directional recorded ground motions from the strong motion seismograph networks K-Net and Kik-Net [63] is used to train and test the HEWFERS framework. More than 85% of ground motions of the database originate from subduction sources. The ground motion component time histories undergo minimal processing, including baseline correction and linear trend removal (only using the samples of the time window of the early raw waves), similar to real-time processing during an earthquake [31], [64]. Ground motion components characterized by PGA value larger than 0.01g are selected for analysis, resulting in approximately 14 000 ground motion components from 1860 earthquake events between 1996 and 2022. The selected earthquakes correspond to events occurring in the megathrust and also within the continental and oceanic plates, thus ensuring the universality of the results. Fig. 2 shows a description of the ground motion database in terms of

M

versus epicentral distance (

R_{e p i})

distribution of the database. The framework is trained using a randomly selected training dataset (80% of the events; shown in green markers in Fig. 2), while evaluations are conducted on the remaining test dataset (shown in red markers in Fig. 2). While a large portion of the data comes from events with

4 < M < 8

, many ground motions are recorded from events with

M > 8,

particularly from the well-known 2003 Tokachi-Oki and 2011 Tohoku earthquakes in Japan. This extensive dataset is thus well-suited for DL and advanced machine learning (ML) approaches.

The detection of the p-waves arrival is a crucial step in the proposed EEWS. Given that the ground motion records do not necessarily start with p-waves and may contain additional noise or zero-padding, an automated p-phase picker algorithm, P_PHASEP_ICKER [65], is used to accurately determine the p-wave arrival time. This method does not require any interval or threshold settings, making it more robust than traditional methods. The algorithm utilizes an SDOF oscillator with a high damping ratio to detect the p-wave onset by tracking changes in the oscillator’s damping energy. The p-wave arrival time, identified through this approach, is then used to define the start of the time history for further analysis, ensuring that only data after the p-wave arrival are used for computing the IMs.

It is important to note, however, that while P_PHASEP_ICKER is used in this study as an aftermath, other p-wave picker algorithms can be employed in real-time applications [66], [67]. These algorithms vary in their approaches and computational requirements. These methods are also designed to accurately identify the p-wave arrival time, and their selection can depend on the specific requirements of the monitoring system, such as processing speed, noise tolerance, and the type of seismic data.

4. HEWFERS training and prediction

The proposed framework’s efficacy is based on the accuracy of three essential components: ① domain-informed VAE (for physics-based surrogacy of

S_{a} (T)

spectrum into LV), ② GPR (for spatial prediction of

{L V}^{n}

using of

{S C}^{1}

{S C}^{n}

{I M}_{10 s}^{1}

, and

d_{n, 1}

), and ③ FFNN (for on-site prediction of

{L V}^{n}

using

{S C}^{n}

and

{I M}_{10 s}^{n}

4.1. Domain-informed variational autoencoder (Phase I)

The log-transformed vector of

P G A

and 84-period

S_{a} (T)

spectrum are transformed into two normally distributed LVs (denoted as

z_{1}

and

z_{2}

with means

μ_{z_{1}}

and

μ_{z_{2}}

and variances

σ_{z_{1}}^{2}

and

σ_{z_{2}}^{2}

, respectively) using a domain-informed VAE [59]. In this study, a domain-informed approach is adopted to enhance the training of the VAE by integrating a Laplacian-based regularization loss. This additional loss component is designed to leverage domain-specific insights, particularly focusing on the intrinsic physical characteristics of earthquake’s

M

and stations’

R_{e p i}

. It penalizes the VAE when

S_{a} (T)

that originate from similar

M

and

R_{e p i}

are mapped far apart in the latent space. This is based on the principle that inputs with closer values of

M

and

R_{e p i}

should have latent representations that are also close to each other, reflecting their physical similarity. This ensures that the geometry of the latent space mirrors real-world correlations and dependencies among the LVs according to the physical attributes of the earthquake event. This physics-informed training thereby improves the interpretability of the latent space and potentially leads to more meaningful embeddings for downstream tasks.

An illustration of the domain-informed VAE is provided in Fig. 3, where

x

represents the true vector of input variables (in this case, 85-point

S_{a} (T)

spectrum) and

\hat{x}

represents the vector of reconstructed input values (in this case, predicted 85-point

S_{a} (T)

spectrum). VAE offers a probabilistic method to represent vectorial observations in their LV space through a neural network-based encoder (recognition model). This encoder is trained alongside a neural network-based decoder (generative model) that utilizes the LV space to reconstruct the observations. As a result, the LV space is designed to have continuous and smooth representations. Consequently, values that are close to each other in the latent space correspond to similar reconstructions by the decoder.

The fundamental concept behind training a VAE is based on Bayes’ theorem, as illustrated in Eq. (1), where

z

denotes the LV space,

p (z)

is the prior distribution, and

p (x | z)

is the conditional probability distribution of

x

given

z

. Given the intractability of probability density of

x

(denoted as

p (x)

), it can be approximated using variational inference [59], which involves estimating conditional probability distribution of

z

given

x

(denoted as

p (z | x)

) with another probability distribution of

z

given

x

(denoted as

q (z | x)

) that is defined to have a tractable distribution. The Kullback–Leibler (KL) divergence (

K L [q |p]

) between any two probability distributions for random variable (

y

) denoted as

p (y)

and

q (y)

is computed using Eq. (2). By designing

q (z | x)

to closely resemble

p (z | x)

, it allows for approximate inference of the intractable distribution. This is achieved by minimizing the KL divergence loss (

{L o s s}_{K L}

) [68] between

q (z | x)

and

p (z | x)

shown in Eq. (3). In this context,

p (z | x)

is the prior distribution

p (z)

which is assumed to be a unit Gaussian distribution

(N (0, 1))

for each LV, where N is the batch size. Eq. (4) represents the reconstruction loss (

{L o s s}_{r e c o n}

), where

E_{q (z | x)} \log p (x | z)

(

E

represents the expectation function) corresponds to the mean squared error (MSE) between the true

x

and the predicted

\hat{x}

values [59].

{L o s s}_{r e c o n}

penalizes the VAE for discrepancies between the true and predicted values while

{L o s s}_{K L}

promotes compact and continuous latent representations.

(1)

p (z | x) = \frac{p (x | z) p (z)}{p (x)}

(2)

K L [q | p] = \sum q (y) \cdot l o g \frac{q (y)}{p (y)}

(3)

{L o s s}_{K L} = | K L [q (z | x) | |p (z | x)] |

(4)

{L o s s}_{r e c o n} = | E_{q (z | x)} \log p (x | z) |

To further make the VAE domain-informed, in this study, the training of VAE is enhanced by incorporating a domain-informed regularization loss (

{L o s s}_{r e g}

), aiming to embed physical insights into the latent representation. This ensures that the latent space not only captures the essential variance and structure for accurate reconstruction but also reflects the physical relationships inherent to the earthquake attenuation problem. This is done by forcing the latent space mappings of

S_{a} (T)

spectra originating from similar

M

and

R_{e p i}

to be close to each other. Consequently, the geometry of the latent space mirrors real-world correlations and dependencies.

Though various regularizations are tested in this case including radial basis function (RBF), Laplacian function, inverse multiquadric, Student’s t-distribution, and so forth. The Laplacian function is used as the final regularization function based on its capacity to regularize the LVs.

{L o s s}_{r e g}

is defined in Eq. (5), where

z_{i}

and

z_{j}

are the LVs of ith and jth data point;

W_{ij}

is weight calculated using the Laplacian kernel, representing the similarity between the ith and jth data points based on

M

and

R_{e p i}

, calculated as per Eq. (6).

M_{i}

and

R_{e p i, i}

represent the

M

and

R_{e p i}

for the ith data point and

M_{j}

and

R_{e p i, j}

represent the

M

and

R_{e p i}

for the jth data point.

τ_{M}^{2}

and

τ_{R_{e p i}}^{2}

represent the variances of

M

and

R_{e p i}

, respectively, that determine the scale of influence for

M

and distance differences respectively. By ensuring that the latent space accurately reflects the physical similarities of

M

and

R_{e p i}

, the model can produce more nuanced and physically meaningful embeddings. Hence the latent space becomes more interpretable with respect to the physical domain. For instance, one can expect that traversing the latent space along certain dimensions would reflect changes in

M

and

R_{e p i}

. Embedding domain-specific knowledge through

{L o s s}_{r e g}

enables the model to produce latent representations that are not only statistically coherent but also meaningful within the specific physical context, facilitating deeper insights and more robust applications of the learned embeddings.

(5)

{L o s s}_{r e g} = \frac{1}{N^{2}} \sum_{i, j} W_{ij} \cdot {‖ z_{i} - z_{j} ‖}^{2}

(6)

W_{ij} = e x p (- \frac{{‖ M_{i} - M_{j} ‖}^{2}}{τ_{M}^{2}} - \frac{{‖ R_{e p i, i} - R_{e p i, j} ‖}^{2}}{τ_{R_{e p i}}^{2}})

(7)

{L o s s}_{T o t a l} = {L o s s}_{r e c o n} + {γ \cdot L o s s}_{K L} + {δ \cdot L o s s}_{r e g}

Hence the model is trained with the total loss (

{L o s s}_{T o t a l}

) as described in Eq. (7). Given the scale of the

{L o s s}_{r e c o n}

(MSE of

S_{a} (T)

spectra in natural logarithm scale) and to allow equal contributions from the three loss terms,

γ

and

δ

of 10 and 10², respectively, are used in this study. The minimization process of Eq. (7) is conducted through stochastic gradient descent and the VAE is optimally trained with hyperparameter tuning [69] (including the number of layers, number of neurons, activation functions, batch sizes, number of epochs, optimizer algorithms, learning rate) and early stopping regularization.

The means of the two LV distributions (

μ_{z_{1}}

and

μ_{z_{2}}

) are presented in Figs. 4(a) and (b), where the colours of the markers represent the

M

of the seismic event (Fig. 4(a)) and

R_{e p i}

of the station site (Fig. 4(b)). It can be observed from the figures that the LVs do seem to have non-linear and interactive trends with

M

and

R_{e p i}

thereby emphasizing the efficacy of the LVs. It should be noted that with this analysis, if the LVs are known one can also deduce the

M

and

R_{e p i}

of the consequential event (perhaps using another regression model). The coefficient of determination (

R^{2}

) between the true

S_{a} (T)

spectrum and constructed

S_{a} (T)

spectrum for 85 periods (covering the important range for the built environment) for both train and test sets using the VAE decoder are presented in Fig. 4(c). For all periods, the

R^{2}

value consistently exceeds 0.85, exhibiting excellent reconstruction power of the surrogate LV with minimal bias and variance. This means that the LVs,

z_{1}

and

z_{2}

, can sufficiently and efficiently reconstruct the

S_{a} (T)

spectrum using the domain-informed VAE-decoder. Hence, for any site k, the VAE-decoder computes

P (S_{a} (T) | {L V}^{k}

) representing the probability density of

S_{a} (T)

conditioned on

{L V}^{k}

(i.e.,

L V

for the site k).

4.2. Feed-forward neural network (Phase I)

As the LVs and VAE-decoder are deemed sufficient and efficient to construct the

S_{a} (T)

spectra, in real-time it is important to estimate the LV quickly. This is done by computing IMs from early recorded p-waves and using them to estimate the LVs. Hence, it requires three key decisions on the parameters: ① which IMs to compute, ② how long of time window to use to compute the IMs after detection of p-waves, and ③ which type of regression model to use.

Based on the previous studies in Refs. [28], [31], [32], [36], [70], [71], [72], the final IMs for computing IM_TW (representing IMs computed for initial time window) include: Arias intensity (

I_{a}

), m·s⁻¹; significant duration (

D_{5 - 95}

), s; mean period (

T_{m}

), s;

P G A

, g; peak ground velocity (

P G V

), m·s⁻¹; peak ground displacement (

P G D

), m; and cumulative absolute velocity (

C A V

), m·s⁻¹. These are described in Eqs. (8), (9), (10), (11), (12), (13), (14), where

a (t)

represents the acceleration time history of the ground motion;

Ti

represents the time instance;

I_{a @ 95 %}

and

I_{a @ 5 %}

represent the 95% and 5% of

I_{a}

, respectively;

C

is the Fourier amplitude spectrum of acceleration at linearly spaced frequencies;

f

is spanning the range 0.25

\leq f \leq

20 Hz. Furthermore, in this study, SC is only quantified using the site’s averaged shear-wave velocity up to 30 m depth (

V_{s 30}

(8)

I_{a} = \frac{π}{2 g} \int {a (t)}^{2} d t

(9)

D_{5 - 95} = T i (I_{a @ 95 %}) - T i (I_{a @ 5 %})

(10)

T_{m} = \frac{\sum C^{2} (\frac{1}{f})}{\sum C^{2}}

(11)

P G A = \max (|a (t)|)

(12)

P G V = \max (|\int a (t) d t|)

(13)

P G D = \max (|\int \int a (t) d t d t|)

(14)

C A V = \int | \int a (t) d t | d t

In order to check the efficacy of initial p-waves IMs to predict the LV during the earthquake event in real-time, various time windows after p-wave detection (including 3, 5, 7, 10, 12, 15, 18, and 20 s) are utilized to compute IMs (denoted as IM_TW). It should be noted that the ground motions used in this study are mostly recorded from subduction earthquakes in Japan. Subduction ground motions tend to have longer duration as compared to crustal sources lasting up to three to five minutes (180 to 300 s) [73], [74]. Fig. 5(a) presents the ground motion duration (

D_{G M}

) for the approximately 14 000 ground motions used in this study where left-hand axis (in green color) presents the histogram and right-hand axis (in red color) presents the cumulative density function (CDF). It is evident that the majority of the ground motions have

D_{G M} >

50 s and 50% of them have

D_{G M} >

120 s.

The computed IM_TW is first utilized to predict LV through a simple five-layered FFNN with 50 neurons in the first layer and a funnel shape with two neurons for LV in the last layer. The FFNN was selected based on prior research [36], where such architectures have been effective for similar tasks. The chosen architecture demonstrated strong predictive accuracy and generalization. The model was designed to provide a balance between expressive power and computational efficiency. The

R^{2}

is recorded for each time window’s FFNN. The results of this analysis are presented in Fig. 5(b) where it can be observed that a time window of 10 s is a good trade-off between having high prediction power for LV and the necessity of a short time window for EEWS application. Table 1 further presents the correlations of the

μ_{z_{1}}

and

μ_{z_{2}}

against the IM_10s for the approximately 14 000 ground motion components time window. It can be observed that most IMs in the IM_10s vector are highly correlated with

μ_{Z_{1}}

and

μ_{Z_{2}}

, emphasizing their importance in the real-time prediction process. Due to the observed collinearity between

C A V

and

P G D

and a low correlation with the

μ_{z_{1}}

P G D

is dropped from the IM_10s vector.

μ_{z_{1}}

and

μ_{z_{2}}

are observed to be almost uncorrelated with each other thereby indicating their unique behaviours.

To estimate LV using SC and IM_10s vectors, four types of regression models were employed: ① linear regression; ② support vector machines [75] (with RBF kernel); ③ XGBoost [76] (with a maximum depth of 10); and ④ FFNN [60]. For all four regressions, the predictors (SC and IM_10s) were transformed to the log domain, while the target variables

μ_{z_{1}}

and

μ_{z_{2}}

underwent a log (

x

+ 5) transformation to accommodate their values, which include both positive and negative values close to zero. The regressions were performed using a training dataset consisting of 80% of the events, randomly selected, with evaluations conducted on the remaining 20% of the dataset. Given the inherent correlation between

μ_{z_{1}}

and

μ_{z_{2}}

, using regression methods ①–③ would result in independent predictions of

μ_{z_{1}}

and

μ_{z_{2}}

, necessitating further postprocessing to explicitly model their correlation. Therefore, an optimized FFNN with two nodes and a linear activation function in the output layer was selected through hyperparameter tuning (considering factors such as the number of layers, number of neurons, activation functions, learning rate, optimizer, and dropout rate) [69] to estimate

μ_{z_{1}}

and

μ_{z_{2}}

using SC and IM_10s.

The true versus predicted

μ_{z_{1}}

and

μ_{z_{2}}

from the trained FFNN for both training and testing sets are presented in Figs. 6(a) and (b), along with their

R^{2}

. The predictions are observed to closely follow the 1:1 purity line, indicating the FFNN’s strong and generalized prediction capability for estimating LV. Beyond high prediction power, another benefit of using FFNNs is the simultaneous estimation of both LVs. The FFNN’s predictions of

μ_{z_{1}}

and

μ_{z_{2}}

have a correlation coefficient of −0.08, compared to the true correlation coefficient of −0.07. This demonstrates the FFNN’s success in maintaining implicit correlations in its estimations. The FFNN provides the on-site LV estimate for any site k (i.e.,

P ({L V}^{k} | {I M}_{10 s}^{k}, {S C}^{k})

) representing the probability density of

{L V}^{k}

conditioned on

{I M}_{10 s}^{k}

(i.e., IM_10s for the site k) and

{S C}^{k}

(i.e., SC for the site k).

Due to the hierarchical nature of the data (i.e., multiple recordings from the same event and multiple recordings from different events), the FFNN residuals are trained with mixed effects [77].

4.3. Gaussian process regression (Phase I)

To provide regional estimates of the LVs and subsequently

S_{a} (T)

, a GPR [61]-based model is employed to develop the spatial relationship of LV. GPR, grounded in the Gaussian process (GP), is a collection of random variables where any finite number of them have a joint (multivariate) Gaussian distribution [61]. GPR is an effective tool for nonlinear regression problems due to its simple structure and its ability to provide non-uniform uncertainty estimations based on feature proximity. Unlike most regression methods, GPR offers efficient inference capabilities for both interpolation (e.g., completing missing data) and extrapolation (e.g., forecasting or prediction), as well as active processes like filtering and smoothing [61]. It accurately captures various relationships between features and targets by utilizing an infinite number of parameters and allowing data to dictate the level of complexity through Bayesian inference. Notably, many Bayesian regression models based on artificial neural networks converge to the Gaussian process with an infinite number of hidden units [78].

The spatial relationship of LV is established using a GPR model trained on earthquake events recorded at more than 10 stations from the training set [61]. The inputs

{S C}^{1}

{S C}^{n}

{I M}_{10 s}^{1}

, and

d_{n, 1}

are utilized to construct the covariance structure of the GPR for spatial prediction of

L V^{n}

. Various kernels were tested to develop this structure, but the final model uses a summation of Matern (

{k (x, x^{'})}_{m a t}

) and White (

{k (x, x^{'})}_{w h i t e}

) kernels [61] as given in Eqs. (15), (17). The Matern kernel, a generalization of RBF kernels, and is parameterized by the length scale (

l

), which must be positive, and the smoothness parameter (

v

). The

l

can either be a scalar (isotropic variant) or a vector with the same number of dimensions as the inputs x (anisotropic variant). Smaller values of

v

lead to less smooth function approximations, and as

v \to \infty

, the kernel approximates the RBF kernel. In Eqs. (15), (17),

d (\cdot)

K_{v} (\cdot)

, and

Γ (\cdot)

denote the Euclidean distance, modified Bessel function, and gamma function, respectively. The White kernel is primarily used to account for noise (

∊

) in the signal, modeled as an independent and identically distributed normal distribution with variance

ϕ^{2}

. Finally,

η

represents the scaling factor for the kernel function

k (x, x^{'})

(Eq. (17)).

(15)

{k (x, x^{'})}_{m a t} = \frac{1}{Γ (v) 2^{v - 1}} {(\frac{\sqrt{2 v}}{l} d (x, x'))}^{v} K_{v} (\frac{\sqrt{2 v}}{l} d (x, x'))

(16)

{k (x, x^{'})}_{w h i t e} = ∊ (N (0, ϕ^{2})); if x = x' else 0

(17)

k (x, x^{'}) = η^{2} \times ({k (x, x^{'})}_{m a t} + {k (x, x^{'})}_{w h i t e})

Given the lognormal nature of seismic demands and IMs, the GPR models are trained in the lognormal domain using a log(x + 5) transformation of the LV [79]. This approach results in an average

R^{2}

of approximately 0.7 for the spatial prediction of both LVs across the training and testing sets. The trained GPR is then utilized in Phase I of the framework to obtain the prior estimate of LV at the target site, leveraging the spatial correlation of LV through distance and SC (i.e.,

P ({L V}^{n} | {L V}^{1}, {S C}^{n}, {S C}^{1}, d_{n, 1})

) representing the probability density of

{L V}^{n}

conditioned on

{L V}^{1}

{S C}^{n}

{S C}^{1}

, and

d_{n, 1}

4.4. Bayesian updating (Phase II)

As discussed in the conceptualization section, Phase II enacts when the waves from the seismic source start arriving at the target nth station. At this point, the framework applies Bayesian updating to the prior LV distribution (

P ({L V}^{n} | {L V}^{1}, {S C}^{n}, {S C}^{1}, d_{n, 1})

) using the on-site LV distribution

(P ({L V}^{n} | {I M}_{10 s}^{n}, {S C}^{n}

)), to obtain the posterior LV and then

S_{a} (T)

at the target nth station.

Bayesian updating is a belief-updating process rooted in Bayes’ theorem of conditional probability [80]. The goal of Bayesian updating is to update prior beliefs about a quantity of interest

θ

(e.g.,

L V^{n} ∣ L V^{1}

) by incorporating new empirical observations

y_{o b s}

. This process involves four key probability density functions (PDFs) termed prior, likelihood, evidence, and posterior, defined as

•Prior distribution (i.e.,

π (θ)

): represents the initial beliefs about the quantity of interest.

•Likelihood function (i.e.,

P (y_{o b s} | θ)

): measures how well the initial beliefs predict the newly observed data.

•Evidence (i.e.,

P (y_{o b s}

)): represents the probability of the model predicting the actual data over the domain of

θ

•Posterior distribution (i.e.,

P (θ | y_{o b s})

): represents the updated beliefs (i.e., PDFs) about

θ

after incorporating the observations.

By incorporating these in Bayes’ theorem, the Bayesian updating framework becomes:

(18)

P (θ | y_{o b s}) = \frac{P (y_{o b s} | θ) π (θ)}{P (y_{o b s})}

which can be expressed as

(19)

P (θ | y_{o b s}) \propto P (y_{o b s} | θ) π (θ)

The posterior distribution is often approximated using simulation algorithms such as Markov chain Monte Carlo (MCMC), which iteratively proposes samples of

θ

and accepts or rejects them based on a probabilistic acceptance criterion [81]. In the context of this work, Bayesian updating provides an idoneous tool to refine the prediction of the LV at the nth target site (i.e.,

P (L V^{n} | L V^{1}, S C^{n}, S C^{1}, d_{n, 1})

) by incorporating new information about the on-site LV distribution (i.e.,

P (L V^{n} | I M_{10 s}^{n}, S C^{n})

). The likelihood function is evaluated through simulation, and conditional updating of the quantity of interest is performed. This process is described as follows:

•Prior beliefs: the prior LV distribution for the target site is derived from the GPR-based prediction, represented by

P ({L V}^{n} | {L V}^{1}, {S C}^{n}, {S C}^{1}, d_{n, 1}) N (μ (μ_{L V}^{G P R}), σ (μ_{L V}^{G P R}))

. Here,

N (μ (μ_{L V}^{G P R}), σ (μ_{L V}^{G P R}))

denotes a normal distribution with mean

μ (μ_{L V}^{G P R})

and standard deviation

σ (μ_{L V}^{G P R})

representing the mean values of the LV estimated by the GPR model, respectively.

•New information: once the initial 10 s of waves reach the target site, the trained FFNN is employed to establish

P ({L V}^{n} | {I M}_{10 s}^{n}, {S C}^{n}) N (μ (μ_{L V}^{F F N N}), σ (μ_{L V}^{F F N N}))

, where

μ (μ_{L V}^{F F N N})

and

σ (μ_{L V}^{F F N N})

represent the mean and standard deviation of the mean LV estimated by the FFNN, respectively.

•The likelihood of observing

μ_{L V}^{T r u e}

(mean value of LV obtained from the initial 10 s of seismic waves; a deterministic quantity) in the prediction is modeled using an error function as

P ({e r r o r (μ}_{L V}^{T r u e}, μ_{L V}^{G P R}) | {e r r o r (μ}_{L V}^{T r u e}, μ_{L V}^{F F N N}) N ({e r r o r (μ}_{L V}^{T r u e}, μ_{L V}^{G P R}), N (0, σ (μ_{L V}^{F F N N}))

. The error function is defined as

{e r r o r (μ}_{L V}^{T r u e}, μ_{L V}^{G P R}) = \sqrt{{(μ_{z_{1}}^{T r u e} - μ_{z_{1}}^{G P R})}^{2} + {(μ_{z_{2}}^{T r u e} - μ_{z_{2}}^{G P R})}^{2}}

, where

μ_{z_{1}}^{T r u e}

and

μ_{z_{2}}^{T r u e}

denote the mean values of both components of LV (i.e.,

z_{1}

and

z_{2}

) obtained from the initial 10 s of seismic waves, and

μ_{z_{1}}^{G P R}

and

μ_{z_{2}}^{G P R}

correspond to the respective counterparts of the GPR model. In this context,

σ (μ_{L V}^{F F N N})

acts as the Gaussian discrepancy model and the observed data or evidence is represented by the error function

{e r r o r (μ}_{L V}^{T r u e}, μ_{L V}^{G P R})

The main drawback of MCMC is the computational time required to approximate the posterior, as the Markov chain needs to be updated sequentially. To mitigate this, certain stopping criteria can be used to terminate the updating process. The first criterion considered in this work is the multivariate potential scale reduction factor (MPSRF) [82], which measures the convergence of multiple independent chains to a unique mean value. Commonly, an MPSRF of 1.02 is acceptable to terminate MCMC. The second criterion considered is the multivariate effective sample size (

m E S S

) [81], which measures the size of uncorrelated samples generated by multiple MCMC chains. After the mESS is greater than a minimum acceptable value (i.e.,

m i n E S S (\hat{θ}, α, ε_{r}))

for a certain confidence level (

α)

at a maximum relative Monte Carlo error (

ε_{r})

, where

\hat{θ}

is the current multivariate posterior sample of all quantities of interest, MCMC can be terminated.

Following this procedure, MCMC is used herein to update the distribution of

L V^{n}

, terminating at a MPSRF of 1.02 and verifying that

m E S S > m i n E S S (\hat{θ}, 0.05, 0.065)

. In this manner, the wall-clock time for MCMC at each station performed on a desktop computer (Intel Core i7 8-core@3 GHz and 16 GB random access memory (RAM)) was about 1.5 s, which is low enough for an EEWS and that could be noticeably improved by using a more powerful machine.

Fig. 7(a) shows the kernel density of the mean LVs of both GPR-based prior and updated posterior distributions normalized by the corresponding true values of LVs for all ground motions. The densities converge around the value of 1.0, indicative of predictions that match the true LVs. The prior distributions exhibit greater variance than their posterior counterparts, pointing to the enhanced precision of predictions after the on-site data assimilation in Phase II. It is observed that the average of this normalized ratio falls very close to one (representing prediction equal to the true value) for both LVs and for both prior and posterior. This indicates the validity of both regional (Phase I) and on-site (Phase II) predictions to provide good estimates for the LV which in turn have high conversion power to the expected

S_{a} (T)

spectrum. The variance of the prior ratios is observed to be higher than the posterior, thereby indicating more accurate predictions made by Phase II of the framework.

Furthermore, Fig. 7(b) presents the coefficient of variation (COV) of the prior and posterior predictive distributions of the prior and posterior estimates, with a notable reduction in COV observed after Bayesian updating. This decrease in COV across the distributions signifies an increase in prediction confidence for both LVs following Phase II.

The COVs are observed to decrease significantly after the Bayesian updating process. This further indicates that Phase II leads to more confident prediction distributions for both LVs. Hence Fig. 7, in general, provides an overview of the prediction power and confidence of the two phases of the proposed framework (i.e., Phase I includes GPR-based prior prediction, and Phase II includes updating through on-site estimate). The proposed framework offers both on-site and regional EEWS capabilities, enabling real-time alerts to be issued to the community. The efficiency allows timely decision-making, facilitating the use of risk-informed EEWS decision support systems and the generation of shake maps.

The unbiased nature of the proposed framework is further assessed by checking the predicted LVs and comparing them against their corresponding

M

and

R_{e p i}

. Fig. 8 showcases the updated posterior means normalized by the corresponding true values of LVs for all ground motions recorded for different events and stations characterized by

M

and

R_{e p i}

. The four sub-figures correspond to the predictive performance of the two LVs against the two source characteristics (i.e.,

M

and

R_{e p i}

). The plots also show bars demonstrating the 5th, 50th, and 90th percentile of the normalized predicted LVs marginalized over

M

and

R_{e p i}

increments of 0.5 and 50 km, respectively. In general, these plots show that the ratios tend to be normally distributed with medians of one across all

M

and

R_{e p i}

. It is further noticed that the prediction of

μ_{z_{2}}

leads to comparatively lower variability than

μ_{z_{1}}

, however, the ratios tend to mostly range between 0.8 and 1.2 indicating strong predictive performance.

Finally, LV computed from the three estimates (i.e., GPR-based prior, VAE + ANN-based on-site, and Bayesian updating-based posterior) is utilized to obtain the

S_{a} (T)

spectra for all approximately 14 000 ground motions. The computed

S_{a} (T)

spectra are compared against the true

S_{a} (T)

spectra through the

R^{2}

metric. This summary is presented in Fig. 9(a). As expected, it can be observed that the spatial regression-based prior leads to relatively low

R^{2}

values with an average of around 0.7. However, with the updating from the on-site estimate, the posterior reaches

R^{2}

values of around 0.9 indicating the high efficiency of the framework. It is further observed in all three cases that the

R^{2}

value drops significantly after PGA for the high-frequency domains (with periods less than 0.4 s). This is adhered to the fact that the ground motions are minimally processed in this case for EEWS, and hence the high noise and variance present in the waveform’s high-frequency domain are not filtered out properly (reflecting the true real-time case). However, this is not expected to critically impact HEWFERS's EEWS applicability since most of the infrastructural systems of interest lie outside this range with periods > 0.5 s. Lastly, to validate LV’s and VAE-decoder’s capability to maintain the cross-correlations within the

S_{a} (T)

spectra, as deemed necessary by several studies [56], [83], the cross-correlation (

ρ

) matrices of the 85-point true and predicted spectra are compared. This is done by computing the difference between the cross-correlations of the true spectra (

ρ_{t r u e}

) and posterior predicted spectra (

ρ_{p o s t e r i o r}

). This is presented in Fig. 9(b). It can be observed that the differences in the cross-correlation lie mainly within the range of −0.1 to +0.1 thereby showcasing the sufficiency of the LV and VAE-decoder to maintain the spectrum cross-correlations.

5. Interpretation of the latent variables

To understand the nature of the developed VAE and FFNN models and to investigate the explainability of LVs, this study utilizes Shapley additive explanations (SHAP) [84]. SHAP is a model-agnostic technique that provides insights into individual predictions by using SHAP values from cooperative game theory. These values are effective in explaining the contribution of each feature to the model’s output. SHAP values uniquely satisfy properties such as efficiency, symmetry, dummy, and additivity [85].

SHAP values essentially quantify the impact of each feature on the respective outputs, similar to how regression model coefficients indicate the influence of features on the target variable. Due to the high computational complexity, SHAP values are approximated using various explainers like kernel-explainer, tree-explainer, and deep-explainer [86]. In this study, a kernel explainer using weighted regression is employed to determine the importance of each feature. These computed importance values function as Shapley values from game theory and as coefficients from a local linear regression, enabling the interpretation of neural network predictions and the effects of feature changes on the target variables [86].

The decoder of the trained VAE is analysed using SHAP with

μ_{z_{1}}

and

μ_{z_{2}}

of the approximately 14 000 ground motions as the inputs (features) and the corresponding predictions of

S_{a} (T)

spectra as the outputs (targets). SHAP values are computed for each input-output combination, resulting in about 14 000 SHAP values for each of the 85 output spectra (

S_{a} (T)

spectra) and two inputs (

μ_{z_{1}}

and

μ_{z_{2}}

). The values are illustrated for

P G A

S_{a} (T = 0.2 s)

S_{a} (T = 1.0 s)

, and

S_{a} (T = 2.5 s)

in Figs. 10(a) to (d). In these figures, the colour of the data points represents the

M

of the feature values, with “low” indicating values close to −3 and “high” indicating values near +2 for both

μ_{z_{1}}

and

μ_{z_{2}}

(based on Figs. 4(a) and (b)).

It can be observed that for all periods, increase in the value of

μ_{z_{2}}

leads to SHAP values moving from negative to positive values. This implies that as

μ_{z_{2}}

increase, it tends to increase the predicted value of

S_{a} (T)

obtained from the VAE-decoder (indicating positive relation between

μ_{z_{2}}

and

S_{a} (T)

for all periods). Similar behaviour is observed for

μ_{z_{1}}

for short period domain (

P G A a n d S_{a} (T = 0.2 s)

), where increase in

μ_{z_{1}}

increases predicted value of

S_{a} (T)

. However, the trend for

μ_{z_{1}}

flips for mid-to-long period ranges (

S_{a} (T = 1.0 s) a n d S_{a} (T = 2.5 s)

). In these cases, it is observed that an increase in

μ_{z_{1}}

tends to decrease the predicted value of

S_{a} (T)

(indicating negative relation between

μ_{z_{1}}

and

S_{a} (T)

for mid-to-long periods). Hence,

μ_{z_{1}}

fluctuates its behavior for stiff and flexible structures thereby indicating period-specific surrogacy of

μ_{z_{1}}

. It is further observed that the variance of SHAP values for

μ_{z_{2}}

tend to grow with an increase in the period (as observed by the spread of SHAP values) while the variance of SHAP values for

μ_{z_{1}}

tend to decrease with an increase in the period. This means that the range of impact of

μ_{z_{2}}

on the predictions increase with an increase in the period. Also, it is observed that the SHAP values generally tend to be symmetric on both sides of zero for both LVs for all

S_{a} (T)

. This means that for both extreme values of the LVs (i.e., −3 and 2), their absolute impact on

S_{a} (T)

is similar (since the absolute SHAP value is similar). In a nutshell,

μ_{z_{2}}

is observed to have a higher impact on all periods with same trend, while

μ_{z_{1}}

tends to have a higher impact on shorter periods with opposite predictive trends for short and long period structures.

Similar to the interpretation process of the VAE decoder, the trained FFNN is analysed using SHAP analysis to relate LV (

μ_{z_{1}}

and

μ_{z_{2}}

) with IM_10s and SC features of the approximately 14 000 ground motions. Hence, approximately 14 000 SHAP values are computed for each combination of the seven inputs (IM_10s and SC) and two outputs (

μ_{z_{1}}

and

μ_{z_{2}}

). Figs. 11(a) and (b) present the SHAP values for the IM_10s and SC corresponding to the two target mean LVs

μ_{z_{1}}

and

μ_{z_{2}}

in descending order of absolute average contribution. The colour of the data points represents the M of the corresponding feature values. It can be observed from Fig. 11(a),

I_{a}

P G V

, and

P G A

lead to the highest SHAP values for

μ_{z_{1}}

while from Fig. 11(b), it is noted that

P G A, C A V

, and

P G V

lead to the highest SHAP values for

μ_{z_{2}}

. Hence it can be concluded

μ_{z_{1}}

contains information regarding the energy and frequency-content of initial ground shaking (as quantified by

I_{a}

and

P G V

) while

μ_{z_{2}}

contains information of the amplitude and energy of initial ground shaking (as quantified by

P G A

and

C A V

). This also explains the trends observed from the VAE-decoder where

μ_{z_{2}}

leads to higher contribution for prediction of

S_{a} (T

) as compared to

μ_{z_{1}}

. The following IMs’ contributions are observed to be similar for

μ_{z_{1}}

and

μ_{z_{2}}

where duration-based

D_{5 - 95}

, site feature

V_{s 30}

, and frequency-based

T_{m}

leads to similar contributions for estimating

μ_{z_{1}}

and

μ_{z_{2}}

From Fig. 11(a) it can be observed that

μ_{z_{1}}

varies positively with the energy-based IM

I_{a}

and amplitude-based IM

P G A

, and negatively with frequency-based IM

PGV

. This means increase in the initial energy and amplitude of the seismic waves increases

μ_{z_{1}}

while as increase in the frequency content leads to decrease in the

μ_{z_{1}}

. Similarly, from Fig. 11(b) it is observed that the prediction of

μ_{z_{2}}

is related to the amplitude-based IM

P G A

in a negative sense and positively proportional to the energy-based IM

C A V

of the initial seismic waves. It should be noted that here negative sense does not mean lower contribution but rather specifies that the corresponding feature value lowers the prediction from the average prediction value of the FFNN.

For a clearer presentation, Figs. 12(a) and (b) display the relative feature importance of

μ_{z_{1}}

and

μ_{z_{2}}

on the VAE decoder predictions, as well as the importance of IM_10s and SC on the FFNN predictions of

μ_{z_{1}}

and

μ_{z_{2}}

, based on their mean absolute SHAP values (

|

SHAP|). The mean absolute SHAP values are calculated for the approximately 14 000 samples for each target and then normalized by dividing them by the sum of the mean absolute SHAP values for each target. Since SHAP values reflect the contribution of features (i.e.,

μ_{z_{1}}

and

μ_{z_{2}}

for VAE and IM_10s and SC for FFNN) to the model outputs, the relative sum of these values in absolute terms indicates the importance of each feature in predicting the target (i.e.,

S_{a} (T)

for VAE and

μ_{z_{1}}

and

μ_{z_{2}}

for FFNN). Fig. 10(a) reveals that

μ_{z_{2}}

’s influence on

S_{a} (T

) significantly increases from shorter to mid-range spectral periods (PGA, 0.2, 1.0, and 2.5 s). This influence peaks at the 0.5 s period and then stabilizes for longer periods (2.5 to 5.0 s), suggesting that

μ_{z_{2}}

is more indicative of the ground motion behaviour where lower frequencies are more prominent. This observation aligns with the understanding that the energy and amplitude of ground motions control long-period behaviour [87], indicating that

μ_{z_{2}}

reflects energy- and amplitude-controlled characteristics of the ground motion. Conversely,

μ_{z_{1}}

generally leads to lower contributions towards the prediction of

S_{a} (T)

. The observed |SHAP| trends are also different for

μ_{z_{1}}

, with its relative importance increasing for

S_{a} (T

) with very short (i.e.,

P G A

) and very long spectral periods (i.e., 5.0 s). This trend indicates that

μ_{z_{1}}

fluctuates significantly with the frequency content of the ground shaking as they significantly influence the behaviour of flexible and long-period structures [90]. For the 0.5 s period,

μ_{z_{1}}

’s importance peaks significantly compared to

μ_{z_{2}}

, indicating complex interactions between the amplitude, energy, and frequency content of the ground motions. These interactions warrant further analysis, which is beyond the scope of this study.

From Fig. 12(b) it is observed that amplitude-based IM

P G A

and energy-based IM

C A V

of the early detected waves have the most dominant influence for

μ_{z_{2}}

, followed by frequency-based IM

P G V

and energy-based

I_{a}

. On the other hand,

μ_{z_{1}}

is observed to be impacted by energy-based IM

I_{a}

and frequency-based IM

P G V

, followed by amplitude-based IM

P G A

and energy-based

C A V

. Both

μ_{z_{1}}

and

μ_{z_{2}}

tend to be affected similarly by the duration-based IM

D_{5 - 95}

, site feature

V_{s 30}

, and frequency-based IM

T_{m}

. Furthermore, it is observed that the variability of influences of IM_10s and SC is greater in

μ_{z_{2}}

as compared to

μ_{z_{1}}

. The prediction of

μ_{z_{1}}

is observed to be dominated by only two out of the seven p-wave IMs while

μ_{z_{2}}

shows a high relationship with four IMs. The relative mean |SHAP| values of IM_10s and SC for

μ_{z_{1}}

show high values for

P G A

I_{a}

, and

P G V

, indicating that the dependency on these measures is stronger compared to others. On the other hand,

μ_{z_{2}}

is observed to be most sensitive to amplitude- and energy- based IMs

P G A

and

C A V

. The other ground motion p-wave IMs and site properties tend to have a uniform dependency on both

μ_{z_{1}}

and

μ_{z_{2}}

Thus, in general, it is observed that

μ_{z_{1}}

is more effective at capturing the frequency and amplitude effects on stiffer SDOFs. As the flexibility of the SDOF increases,

μ_{z_{2}}

becomes more important, reflecting amplitude and energy-based characteristics. This process of interpreting DL-based models offers insights into the cause-effect relationship between features and targets, thereby enhancing transparency and detailing how LVs capture the behaviour of SDOF systems with different periods.

6. Case assessment on the $M$ = 8.0 (2003 Tokachi-Oki earthquake)

HEWFERS is further tested on the 2003 Tokachi-Oki earthquake (

M =

8.0) from the test set of the database. The database contains 101 station recordings of the event with

R_{e p i} < 300

km. The trained VAE-decoder is used to obtain the corresponding LV (i.e.,

μ_{z_{1}}

and

μ_{z_{2}}

) for each ground motion component spectrum. The true computed

μ_{z_{1}}

and

μ_{z_{2}}

are compared against GPR-based prior predictions and updated posterior predictions. The results for each station recording are presented in Figs. 13(a) and (b) for both

μ_{z_{1}}

and

μ_{z_{2}}

. These figures visually summarize the predictive performance of the model by checking the ratio of the predicted mean

μ_{z_{1}}

and

μ_{z_{2}}

and their true values against the stations’

R_{e p i}

. Blue circles depict the spatial regression GPR-based prior predictions’ ratio with the true values. The circles generally cluster around the ideal value of 1.0, indicating good average predictions. Nevertheless, there is notable variability, as indicated by the spread of the blue predictive distribution at various

R_{e p i}

distances (= 100, 150, 200, and 250 km). Hence it is observed that while the prior predictions provide a good average prediction of LV, they are characterized by large variances and substantial uncertainty.

On the other hand, as evidenced by the red squares of Figs. 13(a) and (b), the posterior predictions, which incorporate on-site data, show a marked improvement in accuracy. The updated posteriors yield better average predictions with narrower variances, implying a higher confidence level in the predictions. Notably, the red squares, which are closer to the ratio of 1.0, suggest that the posterior mean predictions are nearly equivalent to the true LV values. The effectiveness of the HEWFERS is evident as both the posterior and prior predictions align closely with the actual recorded LVs, indicating models’ superior predictive capabilities. Furthermore, it is observed neither the prior predictions nor the posterior predictions tend to have any biases with respect to the

R_{e p i}

or distance from the rupture, thereby showing the efficacy of the proposed framework to conduct real-time EEWS for different target stations.

7. Conclusions

The regional earthquake early warning framework, HEWFERS, presented in this study represents a significant advancement in EEWS through the integration of DL and seismological expertise. HEWFERS utilizes a novel, domain-informed VAE architecture that projects the ground motion

S_{a} (T)

spectrum to two physically meaningful surrogate LVs, an FFNN for on-site prediction of the expected ground motion spectrum, and GPR for spatial prediction of the ground motion spectrum across a region in real time. This approach significantly enhances the precision of ground motion intensity predictions, crucial for timely and effective seismic risk mitigation.

The innovative use of XAI techniques within HEWFERS has enabled a transparent and interpretable system, empowering stakeholders to make informed EEWS decisions. The framework has been rigorously tested on a comprehensive database of approximately 14 000 ground motion records, demonstrating its ability to provide accurate and robust predictions of seismic responses. In particular, the framework is tested on two major earthquake events (which aren’t used in the training process). The good prediction capabilities on such large-scale events underscore its versatility and scalability. These capabilities are critically aligned with the UN's disaster risk reduction goals, emphasizing the framework’s potential to substantially contribute to global earthquake preparedness and resilience initiatives.

Although the model was primarily trained using Japanese data, the framework is designed to be scalable and adaptable to other regions. Similar to ground motion prediction models (i.e., attenuation relations), the proposed framework can be trained independently for different geographic regions to capture their localized latent spaces and FFNN structures. This approach not only facilitates regional adaptation but also enhances the understanding of regional seismicity from a mathematical perspective. Given the data-driven nature of the framework, its performance is not expected to change significantly when applied to other regions, as it can adapt to localized variations in seismic conditions.

The successful deployment of HEWFERS could revolutionize the field of earthquake early warning by providing a scalable, accurate, and user-friendly tool for disaster risk management and response. Thus, HEWFERS not only embodies a significant technological leap in the realm of EEWSs but also offers a scalable solution with profound implications for enhancing community resilience and safeguarding lives and infrastructure against seismic threats worldwide.

CRediT authorship contribution statement

Jawad Fayaz: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Conceptualization. Rodrigo Astroza: Writing – review & editing, Methodology, Data curation. Sergio Ruiz: Writing – review & editing, Visualization, Validation, Methodology.

Data availability

The ground motion records utilized in this study are openly available through the Kyoshin Network operated by the National Research Institute for Earth Science and Disaster Resilience (NIED) in Japan. The data can be accessed and downloaded from the Kyoshin Network’s website at https://www.kyoshin.bosai.go.jp/kyoshin/docs/overview_kyoshin_en.shtml. Researchers and interested parties are encouraged to explore and utilize these openly accessible ground motion records for further analysis and scientific investigations.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The first author, Dr. Jawad Fayaz, acknowledges the Department of Computer Science at the University of Exeter, UK, for financial support provided for the high-performance computing. The second author, Dr Rodrigo Astroza, acknowledges the financial support from the Chilean National Research and Development Agency (Agencia Nacional de Investigación y Desarrollo, ANID) through Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) Regular 1240503 and Fondo de Valorización de la Investigación (FOVI)230030 projects. The third author, Dr Sergio Ruiz, acknowledges the financial support from the ANID through FONDECYT Regular 1240501.

References

Publishing order | Descend order by publishing year | Descend order by cited within

[1]	Ritchie H, Rosado P, Roser M.Natural disasters [Internet].Oxford: Our World in Data; [cited 2025 Feb 18]. Available from: https://ourworldindata.org/natural-disasters.

[2]	Swiss R.Global natural and manmade disasters cost insurers US$83B in 2020 [Internet].San Diego: Wells Media Group, Inc.; 2020 Dec 16 [cited 2025 Feb 18]. Available from: https://www.insurancejournal.com/news/international/2020/12/16/594158.htm.

[3]	Xu S, Dimasaka J, Wald DJ, Noh HY.Seismic multi-hazard and impact estimation via causal inference from satellite imagery.Nat Commun 2022; 13(1):7793.

[4]	Bata MH, Carriveau R, Ting DSK.Urban water supply systems’ resilience under earthquake scenario.Sci Rep 2022; 12(1):20555.

[5]	Zhang S, Zhang L, Lacasse S, Nadim F.Evolution of mass movements near epicentre of Wenchuan earthquake, the first eight years.Sci Rep 2016; 6(1):36154.

[6]	León JA, Ordaz M, Haddad E, Ara IFújo.Risk caused by the propagation of earthquake losses through the economy.Nat Commun 2022; 13(1):2908.

[7]	Yunus AP, Xinyu C, Catani F, Subramaniam SS, Fan X, Jie D, et al.Earthquake-induced soil landslides: volume estimates and uncertainties with the existing scaling exponents.Sci Rep 2023; 13(1):8151.

[8]	Guerrero-Miranda P, Luque GA.Social responsibility, sustainability, and public policy: the lessons of debris management after the manabí earthquake in ecuador.Int J Environ Res Public Health 2021; 18(7):3494.

[9]	Mavroulis S, Mavrouli M, Vassilakis E, Argyropoulos I, Carydis P, Lekkas E.Debris management in turkey provinces affected by the 6 February 2023 earthquakes: challenges during recovery and potential health and environmental risks.Appl Sci 2023; 13(15):8823.

[10]	Mavroulis S, Mavrouli M, Lekkas E, Tsakris A.Managing earthquake debris: environmental issues, health impacts, and risk reduction measures.Environments 2023; 10(11):192.

[11]	Gleeson M.Earthquake early warning systems.Nat Rev Earth Environ 2020; 1(2):84.

[12]	Chung AI.The development of earthquake early warning methods.Nat Rev Earth Environ 2020; 1(7):331.

[13]	Cremen G, Galasso C, Zuccolo E.Investigating the potential effectiveness of earthquake early warning across Europe.Nat Commun 2022; 13(1):639.

[14]	Allen RM, Stogaitis M.Global growth of earthquake early warning.Science 2022; 375(6582):717-718.

[15]	United Nations Office for Disaster Risk Reduction (UNDRR).Sendai framework for disaster risk reduction 2015–2030. Report. Boon: United Nations Office for Disaster Risk Reduction; 2022.

[16]	Galasso C, Opabola EA.The 2023 Kahramanmaraş earthquake sequence: finding a path to a more resilient, sustainable, and equitable society.Commun Eng 2024; 3(1):24.

[17]	Kajitani Y, Chang SE, Tatano H.Economic impacts of the 2011 Tohoku-Oki earthquake and tsunami.Earthq Spectra 2013; 29(Suppl 1):457-478.

[18]	Gundogdu S.Turkey’s poor earthquake waste management.Science 2023; 380(6643):353.

[19]	Copenhagen Consensus Center.Natural disasters [Internet].Copenhagen: Copenhagen Consensus Center; [cited 2025 Feb 18]. Available from: https://copenhagenconsensus.com/copenhagen-consensus-iii/natural-disasters-0.

[20]	Kodera Y.An earthquake early warning method based on Huygens principle: robust ground motion prediction using various localized distance‐attenuation models.J Geophys Res Solid Earth 2019; 124(12):12981-12996.

[21]	Hoshiba M.Real‐time prediction of ground motion by Kirchhoff–Fresnel boundary integral equation method: extended front detection method for earthquake early warning.J Geophys Res Solid Earth 2013; 118(3):1038-1050.

[22]	Grapenthin R, Johanson IA, Allen RM.Operational real‐time GPS‐enhanced earthquake early warning.J Geophys Res Solid Earth 2014; 119(10):7944-7965.

[23]	Chung AI, Henson I, Allen RM.Optimizing earthquake early warning performance: ElarmS‐3.Seismol Res Lett 2019; 90(2A):727-743.

[24]	Suzuki W, Aoi S, Sekiguchi H, Kunugi T.Rupture process of the 2011 Tohoku-Oki mega-thrust earthquake (M9.0) inverted from strong-motion data.Geophys Res Lett 2011; 38(7):2011GL049136.

[25]	Zhang J, Zhang H, Chen E, Zheng Y, Kuang W, Zhang X.Real-time earthquake monitoring using a search engine method.Nat Commun 2014; 5(1):5664.

[26]	Kanamori H, Hauksson E, Heaton T.Real-time seismology and earthquake hazard mitigation.Nature 1997; 390(6659):461-464.

[27]	Kuang W, Yuan C, Zhang J.Real-time determination of earthquake focal mechanism via deep learning.Nat Commun 2021; 12(1):1432.

[28]	Cremen G, Galasso C.Earthquake early warning: recent advances and perspectives.Earth Sci Rev 2020; 205:103184.

[29]	Iaccarino AG, Picozzi M, Bindi D, Spallarossa D.On-site earthquake early warning: predictive models for acceleration response spectra considering site effects.Bull Seismol Soc Am 2020; 110:1289-1304.

[30]	Bhardwaj R, Sharma ML, Kumar A.Multi-parameter algorithm for earthquake early warning.Geomatics Nat Hazards Risk 2016; 7(4):1242-1264.

[31]	Caruso A, Colombelli S, Elia L, Picozzi M, Zollo A.An on-site alert level early warning system for Italy.J Geophys Res Solid Earth 2017; 122(3):2106-2118.

[32]	Lin CCJ, Lin PY, Chang TM, Lin TK.Development of on-site earthquake early warning system for Taiwan.S. D’Amico (Ed.), Earthquake research and analysis—new frontiers in seismology, IntechOpen Ltd., London 2012; 329-358.

[33]	Meier M, Kodera Y, Böse M, Chung A, Hoshiba M, Cochran E, et al.How often can earthquake early warning systems alert sites with high‐intensity ground motion?.J Geophys Res Solid Earth 2020; 125(2):1-17.

[34]	Münchmeyer J, Bindi D, Leser U, Tilmann F.The transformer earthquake alerting model: a new versatile approach to earthquake early warning.Geophys J Int 2021; 225(1):646-656.

[35]	Wu S, Beck JL, Heaton TH.ePAD: earthquake probability-based automated decision-making framework for earthquake early warning.Comput-Aided Civ Infrastruct Eng 2013; 28(10):737-752.

[36]	Fayaz J, Galasso C.A deep neural network framework for real‐time on‐site estimation of acceleration response spectra of seismic ground motions.Comput Aided Civ Infrastruct Eng 2023; 38(1):87-103.

[37]	Tajima F, Hayashida T.Earthquake early warning: what does “seconds before a strong hit” mean?.Prog Earth Planet Sci 2018; 5(1):63.

[38]	Wald DJ.Practical limitations of earthquake early warning.Earthq Spectra 2020; 36(3):1412-1447.

[39]	Minson SE, Meier MA, Baltay AS, Hanks TC, Cochran ES.The limits of earthquake early warning: timeliness of ground motion estimates.Sci Adv 2018; 4(3):eaaq0504.

[40]	Minson SE, Baltay AS, Cochran ES, Hanks TC, Page MT, McBride SK, et al.The limits of earthquake early warning accuracy and best alerting strategy.Sci Rep 2019; 9(1):2478.

[41]	Mittal H, Yang BM, Wu YM.Progress on the earthquake early warning and shakemaps system using low-cost sensors in Taiwan.Geosci Lett 2022; 9(1):42.

[42]	Chiang YJ, Chin TL, Chen DY.Neural network-based strong motion prediction for on-site earthquake early warning.Sensors 2022; 22(3):704.

[43]	Campbell KW, Bozorgnia Y.NGA-West2 ground motion model for the average horizontal components of PGA, PGV, and 5% damped linear acceleration response spectra.Earthq Spectra 2014; 30(3):1087-1115.

[44]	Fayaz J, Xiang Y, Zareian F.Generalized ground motion prediction model using hybrid recurrent neural network.Earthq Eng Struct Dyn 2021; 50(6):1539-1561.

[45]	Abrahamson NA, Silva WJ, Kamai R.Summary of the ASK14 ground motion relation for active crustal regions.Earthq Spectra 2014; 30(3):1025-1055.

[46]	McBride SK, Bostrom A, Sutton J, de Groot RM, Baltay AS, Terbush B, et al.Developing post-alert messaging for ShakeAlert, the earthquake early warning system for the West Coast of the United States of America.Int J Disaster Risk Reduct 2020; 50:101713.

[47]	Baker JW.Efficient analytical fragility function fitting using dynamic structural analysis.Earthq Spectra 2015; 31(1):579-599.

[48]	Baker JW, Allin CC.A vector-valued ground motion intensity measure consisting of spectral acceleration and epsilon.Earthq Eng Struct Dyn 2005; 34(10):1193-1217.

[49]	Jalayer F, Beck JL, Zareian F.Analyzing the sufficiency of alternative scalar and vector intensity measures of ground shaking based on information theory.J Eng Mech 2012; 138(3):307-316.

[50]	Bojórquez E, Iervolino I.Spectral shape proxies and nonlinear structural response.Soil Dyn Earthq Eng 2011; 31(7):996-1008.

[51]	Buratti N, Stafford PJ, Bommer JJ.Earthquake accelerogram selection and scaling procedures for estimating the distribution of drift response.J Struct Eng 2011; 137(3):345-357.

[52]	Tothong P, Luco N.Probabilistic seismic demand analysis using advanced ground motion intensity measures.Earthq Eng Struct Dyn 2007; 36(13):1837-1860.

[53]	Bazzurro P, Cornell CA, Shome N, Carballo JE.Three proposals for characterizing MDOF nonlinear seismic response.J Struct Eng 1998; 124(11):1281-1289.

[54]	Fayaz J, Zareian F.An efficient algorithm to simulate site-based ground motions that match a target spectrum.Earthq Eng Struct Dyn 2021; 50(13):3532-3549.

[55]	Fayaz J, Azar S, Dabaghi M, Zareian F.Methodology for validation of simulated ground motions for seismic response assessment: application to cybershake source-based ground motions.Bull Seismol Soc Am 2021; 111(1):226-241.

[56]	Fayaz J, Azar S, Dabaghi M, Zareian F.An efficient algorithm to simulate hazard-targeted site-based synthetic ground motions.Earthq Spectra 2021; 37(2):876-902.

[57]	Vamvatsikos D, Cornell CA.Incremental dynamic analysis.Earthq Eng Struct Dyn 2002; 31(3):491-514.

[58]	Kohrangi M, Vamvatsikos D, Bazzurro P.Site dependence and record selection schemes for building fragility and regional loss assessment.Earthq Eng Struct Dyn 2017; 46(10):1625-1643.

[59]	Kingma DP, Welling M.An introduction to variational autoencoders.Found Trends Mach Learn 2019; 12(4):307-392.

[60]	Rosenblatt F.The perceptron: a perceiving and recognizing automaton. Report. Amherst: University of Massachusetts Amherst; 1957.

[61]	Williams CKI, Rasmussen CE.Gaussian Processes for Regression. In: Proceedings of the 9th International Conference on Neural Information Processing Systems; 1995 Nov 27–Dec 2; Denver, CO, USA. Cambridge: MIT Press; 1995. p. 514–20.

[62]	Holzer TL, Savage JC.Global earthquake fatalities and population.Earthq Spectra 2013; 29(1):155-175.

[63]	National Research Institute for Earth Science and Disaster Resilience.Strong-motion seismograph networks (NIED K-NET, KiK-Net) [Internet].Tokyo: National Research Institute for Earth Science and Disaster Resilience; [cited 2025 Feb 18]. Available from: https://www.kyoshin.bosai.go.jp/.

[64]	Hsu TY, Wu RT, Chang KC.Two novel approaches to reduce false alarm due to non-earthquake events for on-site earthquake early warning system.Comput Aided Civ Infrastruct Eng 2016; 31(7):535-549.

[65]	Kalkan E.An automatic p-phase arrival-time picker.Bull Seismol Soc Am 2016; 106(3):971-986.

[66]	Akazawa T.A technique for automatic detection of onset time of p- and s-phases in strong motion records.In: Proceedings of the 15th World Conference on Earthquake Engineering (15WCEE); 2012 Sep 24–28; Lisbon, Portugal. Tokyo: International Association for Earthquake Engineering (IAEE); 2004.

[67]	Sleeman Eck R, van T .Robust automatic p-phase picking: an on-line implementation in the analysis of broadband seismogram recordings.Phys Earth Planet Inter 1999; 113(1–4):265-275.

[68]	Kullback S, Leibler RA.On information and sufficiency.Ann Math Stat 1951; 22(1):79-86.

[69]	Bergstra J, Bengio Y.Random search for hyper-parameter optimization.J Mach Learn Res 2012; 13:281-305.

[70]	Campbell KW, Bozorgnia Y.Ground motion models for the horizontal components of Arias intensity (AI) and cumulative absolute velocity (CAV) using the NGA-West2 database.Earthq Spectra 2019; 35(3):1289-1310.

[71]	Sand MAıkkaya, Akkar S.Cumulative absolute velocity, Arias intensity and significant duration predictive models from a pan-European strong-motion dataset.Bull Earthq Eng 2017; 15(5):1881-1898.

[72]	Wu YM, Kanamori H, Allen RM, Hauksson E.Determination of earthquake early warning parameters, τ_c and P_d, for Southern California.Geophys J Int 2007; 170(2):711-717.

[73]	Fayaz J, Medalla M, Zareian F.Sensitivity of the response of box-girder seat-type bridges to the duration of ground motions arising from crustal and subduction earthquakes.Eng Struct 2020; 219:110845.

[74]	Fayaz J, Medalla M, Torres-Rodas P, Galasso C.A recurrent-neural-network-based generalized ground-motion model for the Chilean subduction seismic environment.Struct Saf 2023; 100:102282.

[75]	Cortes C.Support-vector networks.Mach Learn 1995; 20:273-297.

[76]	Chen T, Guestrin C.XGBoost: a scalable tree boosting system.In: Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining; 2016 Aug 13–17; San Francisco, CA, USA. New York City: Association for Computing Machinery (ACM); 2016. p. 785–94.

[77]	Demidenko E.Mixed models: theory and applications with R. (2nd ed.), Wiley Publisher, Hoboken (2013)

[78]	Neal RM.Bayesian learning for neural networks. Springer New York, New York City (1996)

[79]	Kramer SL.Geotechnical earthquake engineering. Prentice Hall, Upper Saddle River (1996)

[80]	Kroese DP, Taimre T, Botev ZI.Handbook of Monte Carlo methods. John Wiley and Sons, New York City (2011)-->

[81]	Vats D.Flegal JM, Jones GL.Multivariate output analysis for Markov chain Monte Carlo. 2015. arXiv: 1512.07713.

[82]	Gelman A, Rubin D.Inference from iterative simulation using multiple sequences.Stat Sci 1992; 7(4):457-511.

[83]	Wang M, Wang G.Damping‐dependent correlations between horizontal–horizontal, horizontal–vertical, and vertical–vertical pairs of spectral accelerations.Earthq Eng Struct Dyn 2024; 53(7):2423-2444.

[84]	Lundberg SM, Lee SI.A unified approach to interpreting model predictions. In: Proceedings of the 31st International Conference on Neural Information Processing Systems; 2017 Dec 4–9; Long Beach, CA, USA. Red Hook: Curran Associates Inc.; 2017. p. 4768–77.