A Comparative Study of Damage-Sensitive Features for Rapid Data-Driven Seismic Structural Health Monitoring

Abstract

Rapid post-earthquake damage assessment forms a critical element of resilience, ensuring a prompt and functional recovery of the built environment. Monitoring-based approaches have the potential to significantly improve upon current visual inspection-based condition assessment that is slow and potentially subjective. The large variety of sensing solutions that has become available at affordable cost in recent years allows the engineering community to envision permanent-monitoring applications even in conventional low-to-mid-rise buildings. When combined with adequate structural health monitoring (SHM) techniques, sensor data recorded during earthquakes have the potential to provide automated near-real-time identification of earthquake damage. Near-real time building assessment relies on the tracking of damage-sensitive features (DSFs) that can be directly and rapidly derived from dynamic monitoring data and scaled with damage. We here offer a comprehensive review of such damage-sensitive features in an effort to formally assess the capacity of such data-driven indicators to detect, localize and quantify the presence of nonlinearity in seismic-induced structural response. We employ both a parametric analysis on a simulated model and real data from shake-table tests to investigate the strengths and limitations of purely data-driven approaches, which typically involve a comparison against a healthy reference state. We present an array of damage-sensitive features which are found to be robust with respect to noise, to reliably detect and scale with nonlinearity, and to carry potential to localize the occurrence of nonlinear behavior in conventional structures undergoing earthquakes.

1. Introduction

Despite modern design codes, earthquakes continue to threaten the structural integrity of the built environment and impact the functionality of communities. The lacking seismic design of a salient potion of the existing building stock and the emergence of performance-based design methods, which focus on life safety of building occupants while tolerating structural damage under rare earthquakes [1,2], motivate the search for robust tools enabling the integrity assessment of buildings after major earthquakes. Currently, timely post-earthquake damage assessment and recovery actions are hindered by lengthy expert-conducted visual inspections [3,4], which—despite being increasingly standardized—also suffer from possible subjectivity [5]. However, the tagging of buildings in terms of safety for occupancy forms a much-needed task [6], especially facing the risk of inevitable aftershocks, during which damage may accumulate in mainshock-damaged buildings [7].

A direct quantification of sustained building damage would require the tracking of experienced forces and displacements on critical components, a task which is challenging, if not infeasible, especially when considering the assessment of entire building portfolios. An indirect assessment is, however, made feasible via recent advances in sensor development, which offer dynamic, or vibration-based, sensing hardware at low cost, combined with improved acquisition frameworks, for instance involving Internet of Things [8]. It is however the treatment of these collected raw data, which is key in diagnosing structural condition under operating or extreme conditions [9]. To this end, structural health monitoring (SHM) offers a framework, which relies on the principle of continuous tracking of structural response, in order to detect, localize, and quantify damage, or even attempt to prognose remaining useful life [10]. SHM-based information has the potential to significantly improve upon current post-earthquake evaluation of buildings [11]. Features extracted from vibration measurements have been used by Zhao et al. [12] to track the current state and predict the performance of roller bearings and by Martakis et al. [13] for automated sensor diagnosis. Leveraging the correlation between the stiffness of a structure and the vibration modes, many applications of SHM for the built environment rely on changes in dynamic properties to extract indicators of damage [14,15]. In particular, the link between modal properties, which can be extracted from ambient vibrations (AVs) using operational modal analysis (OMA), and damage has been extensively investigated [16,17,18,19,20,21,22,23,24,25,26,27,28,29]. As a result, despite the limitations of such global indicators of the structural state [30,31], dynamic properties derived from OMA, such as changes in fundamental frequencies, have been introduced into post-earthquake assessment frameworks [32,33,34,35,36,37,38,39]. However, the robustness of such a single feature to encode the structural state has not been assessed.

Modal properties constitute an essential instance of damage-sensitive features (DSFs), i.e., metrics that contain useful information to characterize the condition deterioration of a structure. Such DSFs, often more sophisticated than modal properties, can be extracted from dynamic measurements in near-real time. Such DSFs can subsequently be used to construct damage indices [9], whose probabilistic treatment allows building classifiers for damage detection and localization [40,41,42,43]. Owing to the lack of labeled data, such classifiers are typically calibrated on physics-based models before being applied to real-world data [44]. Shan et al. [45] have introduced a data-driven damage index that incorporates peak and accumulated components and subsequently formulated a link between the damage index and extent of damage, which has been derived from model simulations and successfully applied to real structural elements (concrete columns). Martakis et al. [46] have used domain adaptation techniques to improve the performance of model-based classifiers on real data. However, DSFs were exclusively derived from ambient data, which only contains information about residual stiffness degradation. Moving to the level of prognosis, data-driven indicators can further serve as a starting point for the updating of physics-based structural models [47,48], which can be used to forecast residual capacity after the strike of an earthquake event [49]. Several DSFs have been proposed in the past. Noh et al. [50,51] and Balafas and Kiremidjian [52] proposed wavelet decomposition [53] as a starting point to define DSFs; Johnson and Adams [54] and Luo et al. [55] introduced transmissibility-based DSFs to detect damage; changes in mode shapes and their curvatures were successfully introduced to detect and localize damage [56,57,58,59]; and Kaya and Safak [60] proposed the use of floor spectra for damage localization. However, a thorough comparison of DSFs and their application to highly transient signals, such as the structural response to seismic ground motions, is still lacking.

Many vibration-based damage detection approaches have been applied and tested on high-rise buildings or bridges [61,62,63,64,65,66,67,68]. While slender structures allow for the identification of higher modes, which are often found to be most sensitive to damage and its location [56], this is often impossible for low-to-mid rise buildings, which comprise the most common building typology and may be characterized by high lateral stiffness, such as is the case for un-reinforced masonry buildings, which present the most vulnerable elements within the European building stock. While in many applications, DSFs track changes in modal properties that are derived from AVs, the potential nonlinearity that occurs during seismic events may be undetectable at low amplitudes of shaking due to the closing of the cracks or the compensating effect of multiple damage sources [69]. Indeed, buildings have been found to slowly return to the pre-seismic natural frequencies [70,71], which is a phenomenon that could mask structural weakening. Permanent monitoring installations allow for the continuous tracking of structural response and allow for detecting nonlinearity; this is often separate to residual damage but allows for quantification of the level of cyclic forcing and deformation sustained by a structure, which relates to damaging effects. In this work, we make such a distinction and seek for features, which are not limited to residual damage but allow for the racking of nonlinearity, which is particularly relevant for structures experiencing plasticity and energy dissipation effects.

This paper evaluates the performance of acceleration-based DSFs with respect to their capacity to detect, localize, and scale with the quantity of nonlinearity. The specific contributions can be outlined as follows:

• We offer a systematic and thorough comparison on the capacity of DSFs to detect, localize, and quantify damage using highly transient vibration response data recorded during earthquakes.
• We compare the capacity of such established DSFs to distinguish between temporary nonlinearity and residual stiffness degradation.
• We offer a data-driven evaluation in terms of absolute and incremental damage occurrence during seismic sequences.

We restrict this study to purely data-driven schemes; hence, physics-based simulation models are not used, and exact damage identification is not sought. The robustness of DSFs against noise and transient signals is evaluated. The formulation of the DSFs that are evaluated is provided in Section 2.1. First, the performance of the DSFs on the earthquake-induced response of a simplified simulated spring-mass model is evaluated in Section 3.1, while Section 3.2 presents a comparative assessment on experimental shake-table datasets obtained from the testing of a mixed reinforced concrete unreinforced masonry shear building.

2. Materials and Methods

This paper offers a thorough comparison of acceleration-based DSFs in terms of their efficacy in the detection, quantification, and localization of instances of temporary nonlinearity and residual damage for building structures undergoing earthquake actions. DSFs are features that encode indirect information on the building state in the absence of direct measurements of displacement and reaction forces. The performance of a number of established DSFs, that are overviewed in Section 2.1, is compared in Section 3 in terms of their capacity to detect onset of damage, avoid false alarms, scale with damage, and allow for quantification of damage. Nonlinearity and damage are defined based on engineering demand parameters (EDPs), which are typically used in earthquake engineering to predict the onset and severity of damage and traditional damage indices (DI). The performance of DSFs is tested with the help of EDPs/DIs for a simulated toy model (Section 3.1) and a real building structure tested on a shake table (Section 3.2); however, neither EDPs nor DIs are easily measurable for regional applications of rapid post-earthquake damage assessment.

2.1. Damage-Sensitive Features

DSFs are derived from vibration data and, when defined for a sensor pair that covers the entire structure, can be defined as metrics that characterize the condition (state) of a structure and allow for differentiation at succinct damage levels (see Figure 1). In the context of vibration-based monitoring, which forms the most realistic instance of seismic monitoring, DFSs can be extracted from measured acceleration time-series, possibly in near-real time, which is a feature that conveys their particular potential for rapid loss assessment [11]. As aforementioned, several vibration-based damage detection approaches have been applied and tested on slender or flexible structures, for instance high-rise buildings or bridges [61,62,63,66,67], for which higher vibration modes can be extracted from measurement data. The dynamic response of typical low-to-mid-rise buildings, on the other hand, is often governed by a single dominant mode in each direction, undermining the use of some of the most common existing approaches to damage detection. In addition, while some buildings such as steel or modern reinforced concrete structures offer significant ductility, others, such as masonry buildings, exhibit more sudden and possibly reversible stiffness changes (for instance, rocking motion of slender masonry walls), which may cause a rather brittle behavior and thus requires the use of DSFs that can pick up small variations in the dynamic behavior in short time windows. Therefore, the capacity of DSFs to pick up nonlinearity and/or residual damage will be evaluated in Section 3. In addition, when DSFs are derived for multiple sensor pairs that define substructures of the monitored building, information about the damage location can be inferred (see Figure 1).

In what follows, we present DSFs that are constructed to reveal such variations in the response, mostly corresponding to signs of nonlinear structural responses, which can be associated with damage. The overviewed DSFs are formulated on the basis of five primary quantities, as briefly reviewed in the following sections prior to comparatively assessing their potential in Section 3:

• The first class contains DSFs, whose primary ingredient is the characteristic quantity of transmissibility, which operates in the frequency domain (see Section 2.1.1). Thus, stiffness changes, which may be provoked by nonlinear behavior, are revealed by transmissibility-based DSFs through changes in the (amplitude of the) frequency response function of an input–output signal pair.
• The second class of DSFs derives from correlation and coherence as primary quantity (see Section 2.1.2). These DSFs track the similarity in behavior of two time-series in the time domain (correlation) and frequency domain (coherence). The onset of nonlinear behavior typically appears as loss of correlation/coherence for specific frequency bandwidths.
• Energy-based DSFs (see Section 2.1.3) rely on the time-frequency representation offered by the wavelet decomposition in order to detect changes in energy distribution in the frequency domain and potentially reveal shifts linked to stiffness and damping.
• DSFs that approximate the relationship between input accelerations and output displacements directly (see Section 2.1.4) are used as a representation of the stiffness that links the restoring force to displacement.
• In addition, DSFs that derive from OMA (see Section 2.1.5) can be used to indicate permanent (residual) damage, which alters modal properties, such as natural frequencies and mode shapes.

2.1.1. Transmissibility-Based Features

Transmissibility is defined as the amplitude of the frequency-response function of an input–output pair of dynamic sensors. Typically, a transmissibility value that exceeds 1 indicates amplification and occurs when the input frequency lies in the vicinity of a characteristic frequency of the dynamic subsystem defined by the locations of the sensor pair. In buildings, peaks in transmissibility are linked to the stiffness and mass distribution of the structure above the reference (lowest) sensor. Thus, assuming that the mass does not change significantly in the short term, a change in the transmissibility indicates stiffness reduction, as encountered both during nonlinearity (which can be temporary or recoverable) as well as due to residual damage. In buildings, especially those featuring stiff floor diagrams such as reinforced concrete slabs, each floor can be considered a separate degree of freedom (DOF) and thus, transmissibility input–output pairs can be constructed for each pair of floors. As seismic excitation is ground-borne, the input DOF in this investigation is chosen as the lower floor of the pair.

The use of transmissibility for damage detection has been proposed before [55,72,73], as stiffness reduction in turn provokes a reduction in the frequency of fundamental modes of the substructure defined by the input–output sensor pair, for which transmissibility is calculated using Equation (1):

$$T_{io}\left({\mathit{\omega }}_{\mathit{m}}\right)=\frac{P_{oo}\left({\mathit{\omega }}_{\mathit{m}}\right)}{P_{ii}\left({\mathit{\omega }}_{\mathit{m}}\right)},$$

(1)

where $P_{kk}$ denotes the spectral density of an acceleration signal ${\ddot{x}}_k$, while the subscripts i, o refer to the input and output signal, respectively, and ${\mathit{\omega }}_{\mathit{m}}$ is the frequency band of interest. In the context of purely data-driven assessment, three DSF are formulated to track the change in transmissibility between a healthy reference state and the current (unknown) building state. The first DSF, $Tma{x}_{io,m}$, returns the frequency value that maximizes the transmissibility in a given bandwidth ${\mathit{\omega }}_{\mathit{m}}$:

$$Tma{x}_{io,m}=\underset{\mathit{\omega }\in {\mathit{\omega }}_{\mathit{m}}}{argmax}{T}_{io}\left({\mathit{\omega }}_{\mathit{m}}\right),$$

(2)

where $T_{io}$ is the transmissibility between the input signal (${\ddot{x}}_i$) and the output signal (${\ddot{x}}_o$), as defined in Equation (1). The DSF can then, for instance, be taken as the relative difference between the reference and the new state. While $Tma{x}_{io,m}$ only involves the peak in transmissibility, the frequency centroid, $Tcn{t}_{io,m}$, relates to information from the entire frequency bandwidth, ${\mathit{\omega }}_{\mathit{n}}$, as defined in Equation (3):

$$Tcn{t}_{io,m}=\frac{{\sum }_{i\in m}{\omega }_i{T}_{io}\left({\omega }_i\right)}{{\sum }_{i\in m}{T}_{io}\left({\omega }_i\right)}.$$

(3)

where ${\mathit{\omega }}_{\mathit{m}}$ and $T_{io}$ have already been defined. Again, the change in $Tcn{t}_{io,m}$ can be evaluated as the relative difference between a reference and a target (current) state.

A third transmissibility-based DSF is defined as the vector colinearity, which is used to quantify the similarity between the reference transmissibility and the one measured during a strong ground motion. The vector colinearity is linked to the modal assurance criterion (see Equation (4)), also known as cosine similarity, which is a common metric to assess the similarity between two modal shapes in OMA [74]:

$$TA{C}_{io,m}=1-\frac{{\left|T_{io}^r{\left({\mathit{\omega }}_{\mathit{m}}\right)}^T{T}_{io}^e\left({\mathit{\omega }}_{\mathit{m}}\right)\right|}^2}{\left|{T_{io}^r\left({\mathit{\omega }}_{\mathit{m}}\right)}^T{T}_{io}^r\left({\mathit{\omega }}_{\mathit{m}}\right)\right|\left|{T_{io}^e\left({\mathit{\omega }}_{\mathit{m}}\right)}^T{T}_{io}^e\left({\mathit{\omega }}_{\mathit{m}}\right)\right|},$$

(4)

where the superscript r refers to the baseline reference signal, when the structure is characterized by its original (linear) stiffness, and e denotes an earthquake response signal, with i and o once again denoting the input channel and output channels. In typical applications of seismic SHM, the reference signal is measured under ambient conditions, which can be considered to correspond to an excitation of flat frequency content (white noise). The values of $TA{C}_{io,m}$ range between 0, indicating perfect colinearity between two transmissibility vectors and thus, absence of nonlinearity or damage, and 1, when the two transmissibility vectors $T_{io}^r$ and $T_{io}^e$ are perpendicular. The normalized values of the $TA{C}_{io,m}$ facilitate an intuitive understanding of the corresponding DI.

All three DSFs, whose primary ingredient is transmissibility, depend on a frequency bandwidth (${\mathit{\omega }}_{\mathit{m}}$), which needs to be defined. The first two indicators, $Tma{x}_{io,m}$ and $Tcn{t}_{io,m}$, depend on frequencies that maximize the transmissibility and thus should be restricted to a smaller bandwidth around individual peaks of the transmissibility in order to avoid sensitivity to changes in the input signal. The colinearity metric that defines $TA{C}_{io,m}$ provides some robustness with respect to the relative height of multiple modes and thus can be defined across wider frequency bandwidths.

2.1.2. Coherence and Correlation

The magnitude-squared coherence is used to estimate the power transfer between the input and output of a system. For an ideal constant parameter single-input single-output linear system, the coherence will be equal to one. The magnitude-squared coherence is expressed as a function of the power spectral densities:

$$C_{io}\left(\omega \right)=\frac{{\left|P_{io}\left(\omega \right)\right|}^2}{P_{ii}\left(\omega \right)P_{oo}\left(\omega \right)},$$

(5)

where $C_{io}$ is the magnitude-squared coherence between the signal pair, ${\ddot{x}}_i$ and ${\ddot{x}}_o$ at the frequency $\omega$. $P_{io}$ denotes the cross-power spectral density between the same sensor pair, whereas $P_{ii}$ and $P_{oo}$ refer to the power spectral density of the input and output signal, respectively. In a similar manner to $TAC$ in Equation (4), a DSF to track the changes in the coherence can be formulated as the colinearity between a healthy reference state and a new state for a frequency range ${\mathit{\omega }}_{\mathit{n}}$.

While coherence as a DSF evaluates the correspondence between two signals in the frequency domain, their correlation can also be estimated in the time domain. The wavelet decomposition offers a time-frequency representation that enables tracking correlation for separate frequency bands. As shown by Goggins et al. [75], a decrease in the correlation of the wavelet coefficients outside the frequency range of a particular structural mode may indicate nonlinearity and residual damage. Thus, exploiting the time-frequency representation of the discrete wavelet decomposition, the correlation for a specific frequency bandwidth, ${\mathit{\omega }}_{\mathit{m}}$, can be derived using Equation (6):

$$WC{C}_{io,m}=1-\frac{{\psi }_m\left({\ddot{x}}_i\right){\psi }_m\left({\ddot{x}}_o\right)}{\sqrt{{\psi }_m{\left({\ddot{x}}_i\right)}^2}\sqrt{{\psi }_m{\left({\ddot{x}}_o\right)}^2}},$$

(6)

where ${\psi }_m$ is the wavelet coefficient corresponding to the pseudo-frequency (range) m and ${\ddot{x}}_i$/${\ddot{x}}_o$ are the accelerations measured at the input/output DOFs. To reduce sensitivity to noise, $WCC$ is calculated for a short moving time window, as described in Section 2.2. When the structure behaves linearly, an almost perfect correlation is observed outside the frequency values that correspond to the local structural modes, thus leading to a value of ($WC{C}_{io,m}=0$). When the frequency values of structural modes decrease, then the correlation in the bandwidth below the initial frequency in turn decreases. Thus, this DI is evaluated for a frequency bandwidth corresponding to ${\mathit{\omega }}_{\mathit{m}}=\left[0.5·Tma{x}_{io,m},0.85·Tma{x}_{io,m}\right]$, where $Tma{x}_{io,m}$ is the frequency that corresponds to a peak of the transmissibility function (see Equation (2)) and m designates the mode (or frequency bandwidth) of interest.

2.1.3. Energy-Based Indicators

The wavelet decomposition exposes changes in the energy distribution of a signal in the frequency domain over time, especially in case of non-stationary earthquake signals [76]. The wavelet transform decomposes a time history into a sum of wavelets that are obtained by dilating, scaling and time shifting a mother wavelet. The DSFs investigated in this paper have first been introduced by Noh et al. [50] as robust features.

A first DSF, termed $M{E}_{o,m}$, tracks the spread of the spectral energy as the ratio of the energy accumulated by wavelet coefficients in the region of the $m^{th}$ natural frequency and the total energy accumulated by the structure at the DOF o. To reduce the impact of sensor noise, the total energy, $E_o\left({\mathit{\psi }}_{\mathit{r}\mathit{e}\mathit{f}}\right)$, is calculated for discrete wavelet coefficients within a reduced frequency bandwidth, $\left[f_m/2,f_m\right]$, which is defined based on the frequency of a structural mode, $f_m$, as shown in Equation (7):

$$M{E}_{o,m}=\frac{E_o\left({\psi }_m\right)}{E_o\left({\mathit{\psi }}_{\mathit{r}\mathit{e}\mathit{f}}\right)},$$

(7)

where ${\psi }_m$ denotes the wavelet coefficients related to the (pseudo-)frequency of the $m^{th}$ mode, ${\mathit{\psi }}_{\mathit{r}\mathit{e}\mathit{f},\mathit{m}}$ refers to the reference frequency bandwidth (typically chosen as the range between half of the natural frequency of the $m^{th}$ mode, $0.5·f_m$, and the natural frequency, $f_m$); and E refers the the energy of the signal, which is derived as the absolute amplitude of the complex wavelet coefficients.

To account for the changing energy distribution in the frequency domain that characterizes the seismic ground motion and significantly violates white-noise (WN) assumptions, Equation (7) can be slightly rewritten in relative terms and expressed as a ratio with respect to the same quantity derived for the input signal ${\ddot{x}}_i$:

$$ME{R}_{io,m}=\frac{E_o\left({\psi }_m\right)·E_i\left({\mathit{\psi }}_{\mathit{r}\mathit{e}\mathit{f},\mathit{m}}\right)}{E_o\left({\mathit{\psi }}_{\mathit{r}\mathit{e}\mathit{f},\mathit{m}}\right)·E_i\left({\psi }_m\right)}.$$

(8)

Another feature that exploits the energy decay in the time domain in order to track the residual stiffness deterioration is $C{E}_{io,m}$ [51]. According to Equation (9), $C{E}_{io,m}$ is the time-domain centroid of the signal part beyond $95\%$ of the cumulative input energy. At this point, the measured response comprises mainly free vibrations, as the strong shaking part of the earthquake is past, and thus reflects the post-earthquake dynamic characteristics of the studied structure and can be considered independent from the input shaking. Lower frequencies are associated with longer cycles that shift and result in later time values of the time-domain centroid:

$$C{E}_{io,m}=\frac{{\sum }_{b=b_{95}}^K{E}_{f_m}\left(b\right)·b·t_s}{{\sum }_{b=b_{95}}^K{E}_{f_m}\left(b\right)}-t_{95,i},$$

(9)

where $t_{95},i$ designates the time to reach $95\%$ of the cumulative input energy; $E_{f_m}\left(b\right)$ is the output wavelet energy at (pseudo-)frequency $f_m$ at time sample b; $b_{95}$ is the time sample corresponding to $t_{95},i$, $t_s$ is the sampling period and K reflects the total number of samples that are measured for the earthquake signal.

An alternative formulation of a wavelet-based DSF in time domain is $T{S}_{io,m}$, which is expressed in terms of the difference between the timespans needed to accumulate from 5 to $95\%$ of the output and input signal, respectively [50]. This DSF is derived using Equation (10):

$$T{S}_{io,m}=\frac{\left(t_{95,o}\left({\psi }_m\right)-t_{05,o}\left({\psi }_m\right)\right)-\left(t_{95,i}\left({\psi }_m\right)-t_{05,i}\left({\psi }_m\right)\right)}{t_{95,i}\left({\psi }_m\right)-t_{05,i}\left({\psi }_m\right)},$$

(10)

where $t_{x,i}$ refers to the time required to accumulate $x\%$ of the energy of the input signal and $t_{x,o}$ is the timespan relating to the output signal. To reduce susceptibility to signal noise, the time for accumulation is calculated for the reconstructed signal of the wavelet coefficient ${\psi }_m$, corresponding to the $m^{th}$ mode.

It is mentioned that both $C{E}_{io,m}$ and $T{S}_{io,m}$ act exclusively on the input and output recordings during strong shaking without the need for reference measurements during a baseline (“healthy”) state.

2.1.4. Linear Stiffness Indicators

Two additional DSFs are defined in the time-domain and are computed as (i) an approximation of the structural stiffness and (ii) the prediction error resulting from a linear approximation. The stiffness proxy, $KPRX$, approximates the reaction force, which is defined as the sum of accelerations of all measured DOFs above the input DOF, i, and the displacement, obtained as the relative displacement at the output DOF, o, with respect to the input DOF, i, as formulated in Equation (11):

$$KPR{X}_{io,m}\left(t\right)=\frac{{\sum }_{dof=i}^o{\ddot{x}}_{dof}\left(t\right)}{\widehat{d_o}\left(t\right)-\widehat{d_i}\left(t\right)},$$

(11)

where ${\ddot{x}}_{dof}$ is the acceleration at a given DOF, and $\widehat{d_i}\left(t\right)$ and $\widehat{d_o}\left(t\right)$ are the displacement approximations at time t obtained through numerical integration at the input and output DOF, respectively. Changes in the stiffness proxy, $KPRX$, deliver a direct damage indicator and, when comparing against the $KPRX$ values from healthy reference data, this DI provides a direct measure of stiffness loss. In order to mitigate known limitations and uncertainties induced by numerical integration schemes, prior wavelet-based filtering [77] is performed to reduce noise effects and increase precision.

Another feature is the error induced by linear fitting of the output displacement estimate, $\widehat{d_o}\left(t\right)$, which is defined as $LFI{T}_{io,m}$ in Equation (12):

$$LFI{T}_{io,m}=\sqrt{\widehat{d_o}-{\ddot{x}}_i·KPR{X}_{io,m}},$$

(12)

where ${\ddot{x}}_i$ is the acceleration measured at the input DOF, i; $\widehat{d_o}$ is the displacement estimate at the output DOF, o, obtained from double integration; and $KPR{X}_{io,m}$ is the stiffness proxy, which is defined in Equation (11).

2.1.5. Features Originating from Operational Modal Analysis

OMA methods comprise a super-set of well-known output-only system identification methods (such as Frequency Domain Decomposition [78] and Stochastic Subspace Identification [79] schemes) that enable the extraction of modal properties, i.e., natural frequencies, mode shapes, and—to a lesser extent—damping, from ambient data. Given the non-stationary nature of earthquake excitation and the strong deviation from the WN assumption that defines output-only methods, OMA-based methods are not particularly fit for identifying instances of nonlinearity but can be used for identifying residual (permanent) damage on the basis of treatment of ambient data sets prior to and after an earthquake event.

In addition to two classical DSFs, namely the change in natural frequency $\Delta f_m$ of the m-th mode and the corresponding mode shape $MA{C}_m$, evaluated using the MAC criterion, two additional OMA-based DSFs are typically exploited: the change in mode-shape curvature, $\Delta MSC$, and the change in the phase angle, $PAAC$.

The mode-shape curvature (MSC), which is the second spatial derivative of the modal shapes, is estimated using the central difference approximation formula:

$$MS{C}_{k,m}=\frac{{\varphi }_{k-1,m}-2{\varphi }_{k,m}+{\varphi }_{k+1,m}}{h^2},$$

(13)

where ${\varphi }_{k,m}$ is the the k-th component of the mode-shape vector $\varphi$ of the m-th mode. The absolute change in the MSC, which provides a feature that is sensitive to residual damage and its location, is calculated with respect to the healthy (reference) state:

$$\Delta MS{C}_{k,m}=\frac{MS{C}_{k,m}^h-MS{C}_{k,m}^d}{MS{C}_{k,m}^h},$$

(14)

where h and d refer to the healthy reference and the potentially damaged post-earthquake building state.

Finally, the $PAAC$ is a DSF based on variations in the phase angle of an input–output sensor pair. The phase angle is estimated through Equation (15) on the basis of the the cross-power spectral density (CPSD) matrix, which can be obtained using Fourier analysis:

$$P{A}_{io}\left({\mathit{\omega }}_{\mathit{m}}\right)=arctan\left(\frac{imag\left(CPS{D}_{io}\left({\mathit{\omega }}_{\mathit{m}}\right)\right)}{real\left(CPS{D}_{io}\left({\mathit{\omega }}_{\mathit{m}}\right)\right)}\right),$$

(15)

where $CPS{D}_{io}$ is the CPSD evaluated for an input (i)–output (o) sensor pair; ${\mathit{\omega }}_{\mathit{m}}$ is the frequency bandwidth defining the m-th mode of the frequency-response function between i and o; and $imag$/$real$ refer to the imaginary/real part of the complex CPSD vector, respectively. The $PAA{C}_{io,m}$ reports the change in the $P{A}_{io}$ between a healthy reference state and a potentially damage post-earthquake state and is computed in a similar manner to the TAC, as defined in Equation (4).

2.1.6. Summary of DSFs

The DSFs reviewed in the previous sections are summarized in

Table 1. Summary of the damage-sensitive features (DSFs) reviewed in Section 2.1.

DSF	Domain	Signal	Nonlineariy	Residual Damage	Localization	Equation
$TAC$	f	AV & GM	✔	✔	(✔)	(4)
$Tmax$	f	AV & GM	✔	✔	(✔)	(2)
$Tcnt$	f	AV & GM	✔	✔	(✔)	(3)
$WLC$	tf	AV & GM	✔	✔	✔	(6)
$CohAC$	f	GM	✔			(4)
$CE$	tf	GM	✔		✔	(9)
$TS$	tf	GM	✔		✔	(10)
$ME$	tf	AV & GM	✔	✔		(7)
$MER$	tf	AV & GM	✔	✔	✔	(8)
$KPRX$	t	AV & GM	✔	✔	✔	(11)
$LFIT$	t	GM	✔		✔	(12)
$\Delta f_n$	f	AV		✔
$PAD$	f	AV		✔	✔	(15)
$MAC$	f	AV		✔
$\Delta MSC$	f	AV		✔	✔	(14)

, indicating whether their principle is acting in the time or frequency domain, whether they apply to the ground motion (GM) signal or to post-earthquake AV data, whether they evaluate nonlinearity, residual damage, or both, and whether they offer the possibility of localizing the damage (see Figure 1).

2.2. Extraction of Damage Indicators

The DSFs, presented in the previous Section 2.1, are either extracted from AVs or on the basis of the GM signal. In the context of continuous monitoring, the separation of the GM-induced response from AVs requires temporal delimitation of the continuously acquired GM signal. Therefore, in this work, a triggering condition based on the ratio between short-time average and long-time average (STA/LTA) of absolute acceleration values, recorded at the base of a building, is used [39,80], as indicated in the top subplot of Figure 2. Then, the GM signal is split into shorter time windows of a fixed pre-defined window length in order to track the evolution of nonlinearity and pick up temporary shifts in the evaluated DSFs (see bottom subplot of Figure 2).

As the majority of DSFs involve discrepancies between reference (typically “healthy”) and newly acquired data, it is essential to adopt appropriate metrics for comparison. Due to the uncertainties that underlie monitoring data, primarily due to noise contamination, probabilistic metrics offer a more complete representation of the observed discrepancies.

When comparing a single DSF value with a reference distribution, the Mahalonobis distance [81] can be adopted, which reflects a point-to-distribution metric, and is derived using Equation (16):

$$d_M=\sqrt{{\left(\mathbf{x}-\mathit{\mu }\right)}^T{\Sigma }^{-1}\left(\mathbf{x}-\mathit{\mu }\right)},$$

(16)

where $\mathbf{x}$ is the vector of DSFs characterizing a new state, $\mathit{\mu }$ is the vector containing the mean of the reference (“healthy”) distributions of the DSFs and $\Sigma$ is the corresponding covariance matrix. The traditional formulation of the Mahalonobis distance (see Equation (16)) combines multiple DSFs and adopts a Euclidian distance that is agnostic to the direction of a change. However, the adoption of a signed distance (which reflects the direction of a change) may be important for certain DSFs; for instance, in the case of $KPRX$ and $Tcnt$, an increase in stiffness is unlikely to result from nonlinearity. Therefore, a slightly differentiated formulation of $d_M$ is proposed, per the normalized distance for a single DSF, as presented in Equation (17)

$$d_{MS}\left(D{I}^e\right)=\frac{D{I}^e-\mathit{\mu }\left(D{I}^h\right)}{\sigma \left(D{I}^h\right)},$$

(17)

where $\mathit{\mu }$ and $\sigma$ reflect the mean and standard-deviation, respectively, and the superscripts e and h refer to the newly acquired and reference (“healthy”) data, respectively.

A second probabilistic distance metric can be defined on the basis of splitting the acquired records into shorter time-windows (as shown in Figure 2), or by considering multiple windows of AV. The Kullback–Leibler (KL) divergence [82,83,84] is defined as the difference between two probabilistic distributions; a reference (“healthy”) distribution, $f^h\left(x\right)$, and the distribution corresponding to newly acquired data, $f^e\left(x\right)$:

$$di{v}_{KL}\left(DSF\right)=\sum _x{f}_{DSF}^h\left(x\right)log\left(\frac{f_{DSF}^h\left(x\right)}{f_{DSF}^e\left(x\right)}\right)$$

(18)

where the superscripts h and e again denote a healthy reference state and a newly acquired unknown condition, respectively; $f_{DSF}\left(x\right)$ refers to a fitted distribution to the computed values of a specific DSF, on the basis of the individual time-windows. The fitting of a probability distribution allows to reduce sensitivity to noise and transient effects in short time-windows.

3. Results

3.1. Simulated Case Study

3.1.1. Model Simulations

In this section, we exploit a simulated toy example, namely a nonlinear spring-mass model, to offer a first comparison of the DSFs presented in Section 2.1 in terms of their efficacy in detecting the onset of nonlinearity and residual damage as well as localizing and quantifying nonlinearity/damage within elements of a building structure. The model comprises five DOFs, emulating a five-story shear building. The DOFs are connected by means of hysteretic Bouc-Wen springs [85], with the nonlinear time-history analyses carried out in the Opensees software [86]. The masses that are lumped into each of the DOFs are $\mathbf{M}=\left\{117.7,124.2,103.7,119.3,116.6\right\}$ tons, and the linear and nonlinear properties of the springs connecting the DOFs are summarized in

Table 2. Definition of the nonlinear link properties of the simulated mass-spring model.

Link	Stiffness (kN/mm)	Yield Force (kN)
0-1	270	500
1-2	250	433
2-3	230	350
3-4	210	300
4-5	190	225

. In addition, the following properties define the Bouc-Wen hysteretic model of the springs: the exponent regulating the transition between linear and nonlinear branch: $n=5$, stiffness degradation $\delta \eta =3.{610}^{-5}$, strength degradation $\delta \nu ={10}^{-3}$, and hardening of $2\%$. In addition, a Rayleigh-type damping (with a damping coefficient of $\zeta =3\%$ for the first two modes) is adopted. A schematic representation of the model and the behavior of the link 0-1 are shown in Figure 3.

Given the natural variability of GM signals, changes in spectral content and signal duration may undermine the capacity of DSFs to separate linear from nonlinear behavior. In an effort to cover a realistic range of GMs, 21 ground motions from historical earthquakes are scaled to multiple amplitudes for generating a structural response dataset with the adopted nonlinear model. An overview of the earthquake GMs is shown in Figure 4, with the details of the utilized historical earthquake records summarized in

Table A1. Earthquakes used in the simulated examples.

Earthquake	Year	Station	Fault Type	${\mathit{R}}_{\mathit{rup}}$	${\mathit{M}}_{\mathit{w}}$	Database
L’Aquila (IT)	2009	Gran Sasso (Assergi)	Normal	6.4	6.30	PEER
Montenegro (MG)	1979	Ulcinj—Hotel Albatros	Reverse	4.35	7.10	PEER
Gilroy (US)	2002	Palo Alto (Fire Stat)	Strike slip	72.49	4.90	PEER
Northridge (US)	1994	Bev. Hills 12520 Mulhol	Reverse	18.36	6.69	PEER
Imperial Valley (US)	1979	Delta	Strike slip	22.03	6.53	PEER
Christchurch (NZ)	2011	Riccarton High School	Rev. Obl.	9.44	6.20	PEER
Umbra Marche (IT)	1997	Nocera Umbra	Normal	8.92	6.00	PEER
Kozani (GR)	1995	Kozani	Normal	19.54	6.40	PEER
Kobe (JP)	1995	Fukushima	Strike slip	17.85	6.90	PEER
Chi-Chi (TW)	1999	CHY024	Rev. Obl.	9.62	7.62	PEER
Loma Prieta (US)	1989	Capitola Fire Station	Rev. Obl.	12.53	6.93	PEER
Friuli (IT)	1976	Buia	Reverse	11.03	5.91	PEER
Izmit (TK)	1999	AI_081_ IZN_KY	Strike slip	40.3	7.60	ESMD
Erzincan (TK)	1992	2402	Strike slip	0.9	6.60	PEER
Aigio (GR)	1995	AMIA	Normal	16.6	6.50	PEER
Adriatic Sea (ALB)	2019	DURR	Thrust	4.4	5.60	ESMD
Tirana (ALB)	1988		Strike slip	7	5.90	REXEL
Vallo di Nera (IT)	1980	NRC	Strike slip	10	5.00	ESMD
Amatrice (IT)	2017	MSC	Normal	8.2	5.60	ESMD
Colfionto (IT)	1997	CLC	Strike slip	1.2	4.30	ESMD
Lytle Creek (US)	1970	Cedar Springs Pump.	Rev. Obl.	22.94	5.33	ESMD

, which can be found in the Appendix. As the aim of the analysis is to assess the robustness of DSFs to identify earthquake-induced damage, the GMs are deliberately chosen to represent several faulting styles and source-to-site distances, yet their scaling from peak-ground-acceleration (PGA) values of 0.5 m/s${}^2$ to 13.5 m/s${}^2$ may not always reflect realistic earthquake GMs. The variability of the respective frequency spectra is shown in Figure 4 (top right), together with the scatter in two intensity measures characterising the earthquakes (Figure 4, bottom right): first, the duration during which $50\%$ of the GM energy is accumulated ($D_{25-75}$), and secondly, the ratio between the spectral acceleration at the fundamental frequency ($S_a\left(T_1\right)$) of the modeled structure and the PGA.

With the help of the nonlinear time-history model simulations, the DSFs, described in Section 2.1, are reviewed with respect to the capacity to detect and quantify nonlinearity/damage (Section 3.1.2), to locate damage (Section 3.1.3), and in terms of their sensitivity to noise (Section 3.1.5). Taking advantage of the simulated model, three EDPs (namely maximum transient roof drift ratio, hysteretic work, and residual stiffness drop), and one damage index (DI), namely the $D{I}_{Park\&Ang}$ are used for cross-validation. It should, however, be stressed that these EDPs cannot be directly measured or straightforwardly quantified at the required resolution for a reasonable cost, using current sensing solutions. Nevertheless, these EDPs are widely used in defining seismic damage sustained by a structure and its influence on residual performance.

3.1.2. Damage Detection and Quantification

The first level of damage identification consists in detecting the onset of damage (shift in stiffness) or nonlinearity (i.e., regions where the response lies outside the linear elastic range). Damage detection is first executed in a global sense. This implies that the DSFs presented in Section 2.1 are derived for the entire structure, with input acceleration measured at the base and output acceleration acquired at top storey (roof) level. This setup is most suitable for low-rise buildings, which are generally dominated by the fundamental vibration mode. If higher modes bear a more prominent role, for instance in high-rise buildings, the output (measured) acceleration can be assigned at an intermediate level, where higher modes contribute further to the deformation.

The potential of different DSFs for damage identification is evaluated using four performance indicators (PIs), as illustrated in Figure 5:

• The correlation between a DSF and an EDP, such as maximum transient roof displacement or the DI by Park&Ang [87], is defined using Equation (19):

  $$D{I}_{PA}=\frac{d_{max}-d_y}{d_u-d_y}+\beta \frac{E_h}{d_u{V}_y},$$

  (19)

  where $d_{max}$ is the absolute maximum transient displacement the structure undergoes during an earthquake, $d_u$ is the ultimate displacement capacity of the structure, $d_y$ is the idealized yield displacement of the structure, $E_h$ refers to the hysteretic energy that is dissipated by the structure, $V_y$ is the idealized yield force, and $\beta$ is a weighting factor that characterizes the influence of energy dissipation on the overall damage. Herein, it is assumed that $\beta =0.15$.
• The capacity of a DSF to distinguish between building responses of linear buildings ($D{I}_{PA}<0.025$) and buildings that sustained significant nonlinearity ($D{I}_{PA}>0.50$). This is quantified as the probability of exceeding the $D{I}_{PA}$ thresholds given a DSF value, $p(di⩾DI\mid DSF)$, which is computed using the so-called IM-based method, as described by Iervolino [88] and illustrated in Figure 5a,b. The efficacy in distinguishing between separate damage states is evaluated as the probability of exceeding a $D{I}_{PA}$ of $0.50$, $p(di⩾0.50\mid DSF)$ for the DSF value that produces a probability of $90\%$ to exceed $D{I}_{PA}=0.025$.
• The third PI evaluates the probability of erroneous damage detection. To this end, the number of GMs, for which the detection threshold (derived using reference AV data) is exceeded, is calculated. This procedure is described in Figure 5c.
• Finally, the fourth PI is related to the minimum value of EDP/$D{I}_{PA}$, for which the DSF detects damage (see Figure 5d). This PI is evaluated as the mean EDP/$D{I}_{PA}$ value, at which the DSF crosses the detection threshold. Similarly to the previous PI, this is defined using reference (“healthy”) AV data. In both cases, the detection threshold is set to the $DS{F}^h$ value corresponding to a cumulative probability of $0.95$ for the healthy reference distribution.

TAC values are derived for a specific frequency bandwidth, which influences the results, as described in Appendix Figure A1. The results of Figure 5 are obtained for the full bandwidth ($fb=4–25\mathrm{Hz}$) that is tightly defined around the first and second mode. In addition, the TAC is obtained by averaging the values of four-second windows and a time-bandwidth parameter of the Morse wavelet set to 100 (see Appendix B Figure A1).

For the robust detection of damage or tracking of the manifestation of nonlinearities during earthquake loading, the DSFs should scale with the amount of nonlinearity sustained by the structure and be invariant to the characteristics of the earthquake excitation. This forms a challenge as increasing shaking amplitudes naturally lead to higher amounts of nonlinearity, and an often-biased correlation to earthquake amplitudes may occur. For the tracking of nonlinear exposure, in particular, the proposed DSFs should be adept in discriminating between linear and nonlinear structural responses. To investigate such potential, performance indicators are analyzed for 21 GMs and scaled to increasing amplitudes that lead to changing structural conditions (from purely linear behavior to extensive damage). When simulating the “healthy” reference behavior, a pure WN excitation would lead to DSF values that are significantly differentiated with respect to realistic ambient or earthquake ground inputs, which typically comprise varying spectral and temporal components. In order to reflect realistic conditions, four 10 min AV datasets that have been measured at the base of an actual building in Switzerland during day-and-night time are used as input for simulation rather than artificial WN excitation.

Figure 6 illustrates the values of $Tcn{t}_{05,M1}$ with respect to four EDP values: hysteretic work (left), maximum transient roof displacement (center left), $D{I}_{PA}$ (center right), and the residual stiffness drop (right). It is observed that the correlation is very good and the scatter due to the GM signals remains low. However, the significant uncertainty in the ambient reference distribution ($Tcn{t}_{05,M1}^h$) undermines the successful detection of damage.

Reversely, the DSF based on the phase angle has a very low scatter in the ambient conditions and thus allows for a significantly better detection, even from lower levels of damage (see Figure 7). As the $PAAC$ is derived from post-earthquake AVs, the correlation with the residual stiffness drop is the highest (see Figure 7, right). In addition, when comparing Figure 6 and Figure 7, the variability of the $PAAC$ for linear model instances is much lower than for the $T_{cnt}$, highlighting the superior capacity of the $PAAC$ to identify absence of damage and the reduced influence of GM signals, given the $PAAC$ is derived from AVs.

The systematic evaluation of all DSFs based on the four defined PIs—derived with respect to $D{I}_{PA}$, which combines the influence of hysteretic work and maximum transient displacement—is shown in Figure 8. The evaluation of all DSFs allows ranking the DSFs with respect to their performance to detect and quantify structural damage. When focusing on the $T_{cnt}$, shown in detail in Figure 6, the good correlation with $D{I}_{PA}$ is reflected by a high ranking in the two first PIs, which evaluate damage quantification. In addition, the large variability in the healthy reference state avoids false positives (see PI3 in Figure 8) at the cost of very reduced capacity of early detection (see PI4 in Figure 8). While similar conclusions can be drawn for the quantification under use of the $PAAC$, the much lower variability of the healthy reference values increase the risk for erroneous damage detection yet significantly increase the capacity for early damage detection.

When comparing the capacity to quantify damage (PI1 and PI2 in Figure 8), the transmissibility- and the OMA-based DSFs yield superior results, while the DSFs that are based on the energy distribution have a poor performance. This is due to the inevitable sensitivity of the energy distribution with respect to the spectral content of the GM. The good performance of the OMA-based DSFs originates from the fact that AVs have reduced variability in comparison with DSFs that are derived from GMs that are transient in nature and thus characterized by significant variability in the frequency domain (as observed in Figure 4, top right). In terms of damage detection potential, as reflected by PI3 and PI4 in Figure 8), we observe a trade-off between probability of false detection (PI3) and early detection of damage. This underlines the need for combining multiple DSFs when constructing robust damage identification schemes.

As shown in Figure 8, DSFs that are derived from AVs have good overall performance (PAAC is the best DSF overall with respect to $D{I}_{PA}$) due to the reduced uncertainty from the GM signals. Transmissibility-based DSFs can be derived from both GMs and AV; Figure 9 contains a comparison of the respective values for all three DSFs that derive from changes in transmissibility: $T_{cnt}$, $T_{max}$, and $TAC$ on the right. While the scatter is low for $T_{cnt}$ (Figure 9 left) and the values are centered around the the diagonal, a noticeable scatter appears for $T_{max}$ (Figure 9 center). The largest difference is found for the $TAC$ (Figure 9 right), which is more sensitive when derived from GMs and picks up damage at earlier stages than when derived from post-earthquake AVs. Indeed, almost all values lie above the diagonal. This behavior derives from the fact that the temporary stiffness drops created by nonlinearity are more pronounced than permanent stiffness drops related to damage (see Figure 3 bottom right). The differentiation in DSFs derived from AVs versus those extracted during strong ground motions is further analyzed in Section 3.1.5.

3.1.3. Damage Localization

Beyond detection of the onset damage, the spatial information of the sensing network that is encoded within the various DSFs may also provide insights into the location of damage. In doing so, we resort to the use of input–output signal pairs, which are defined locally, in order to assess the onset of damage at the substructural level. In buildings, a reasonable level of localization may be achieved through floor-wise damage assessment, since the need for repair quantities of main structural elements, such as columns/walls, is regulated on the basis of the experienced inter-story drift ratio [89]. While some OMA-based features, namely the difference in natural frequencies ($\Delta f$) and mode shapes ($MAC$), do not contain information about the location of damage, the difference in mode-shape curvature ($\Delta MSC$) does contain information about the location, but rather than relying on local sensor pairs, it requires the systematic monitoring of multiple DOFs (in case of buildings, all floors should be measured for instance).

In order to assess the capacity of DSFs to achieve floor-wise damage localization, five alternate mass-spring models are created by limiting nonlinearity to a single link (storey) per model instance, while the other links remain linear elastic. The local evaluation of five DSFs for local input–output pairs is reported in Figure 10 for the Amatrice earthquake, and then, it is scaled to a PGA that leads to a hysteretic work of approximately 20 kNmm (analysis of all earthquake signals is done later on). Three models, differing in the location of the nonlinear link (left 01, center 23, and right 45), are compared (see Figure 3 for model definition). For ease of comparison, the DSFs are reported in terms of the modified Mahalonobis distance (see Equation (17)).

The first row of Figure 10 reports the evolution of the first examined DSF, namely the TAC. The transmissibility of a sensor pair characterizes the entire substructure from the input sensor (i) to the rooftop, irrespective of the output sensor (o). Therefore, while for the first damaged case (Link 01) the TAC correctly identifies damage as limited to the first link (top left of Figure 10), more sensor pairs tend to pick up damage when this is located above their input level, i. This is observed in the cases of damage in links 23 and 45, respectively. However, an absence of damage in links 34 and 45 is correctly identified in the case of damage in link 23.

The second row of Figure 10 plots the evolution of the DSF based on the linear correlation of wavelet coefficients, $WLC$. This DSF offers high localization precision, with distinct peaks occurring at the corresponding damage locations. The sensitivity to damage is the lowest at the ground floor, where linear correlation is influenced more by the ground motion than the structure, unlike in higher floors.

The DSF based on energy accumulation in the time-domain, $TS$, is reported in the third row of Figure 10. This further proves sensitive to the location of damage; however, spurious peaks appear when the damage is located above the first floor, where the energy accumulation is less influenced by the ground motion itself. In addition, when the damage occurs in the last link 4–5, it is not picked up.

Perhaps the most broadly established DSF for damage localization, in a more global sense, the change in mode-shape curvature, $\Delta MSC$ is formulated with respect to DOFs and not sensor pairs. This implies that this DSF does not directly pinpoint damaged links. However, when the damage appears at an extreme of the structure, either at the ground or at the top, the location is correctly identified, as shown in the fourth line of Figure 10. When damage appears at an intermediate level, the localization is less clear. Given that for unreinforced masonry buildings, which comprise the most vulnerable building type in central Europe and are a major contributor to seismic risk worldwide, the most common failure mechanisms are a soft story at the ground floor or the out-of-plane failure of walls in the highest floors, this indicator—albeit global in nature—can still be useful for rapid post-event assessment. It should be noted that the change in mode-shape curvature of the first mode (plotted in the fourth line of Figure 10) is most sensitive to damage in the lower links and is insufficient to achieve localization for all DOFs (see

Table A2. Theoretical values of OMA-based DSFs, when damage (stiffness reduction of 10%) is applied in single links of the spring-mass model, which is used in Section 3.1. The change in the first three natural frequencies ($\Delta f_n,n=1–3$) is reported as well as the mode-shape curvature, $\Delta MS{C}_{ij}$, of the first two modes ($j=1–2$) at the first four degrees of freedom ($i=1–4$).

Property	Ref.	Dam ${\mathit{K}}_{01}$	Dam ${\mathit{K}}_{12}$	Dam ${\mathit{K}}_{23}$	Dam ${\mathit{K}}_{34}$	Dam ${\mathit{K}}_{45}$
$f_1$	6.588	6.475	6.483	6.509	6.539	6.573
$\Delta f_1$		−0.017	−0.016	−0.012	−0.007	−0.002
$f_2$	18.248	17.962	18.190	18.175	17.928	17.991
$\Delta f_2$		−0.016	−0.003	−0.004	−0.018	−0.014
$f_3$	28.926	28.665	28.865	28.439	28.938	28.192
$\Delta f_3$		−0.009	−0.002	−0.017	0.000	−0.025
$\Delta MS{C}_{1,1}$		−231.227	224.273	3.818	2.364	0.773
$\Delta MS{C}_{2,1}$		−0.008	−1.108	1.069	0.045	0.012
$\Delta MS{C}_{3,1}$		0.026	0.021	−0.603	0.552	0.015
$\Delta MS{C}_{4,1}$		0.038	0.033	0.020	−0.254	0.162
$\Delta MS{C}_{1,2}$		−0.227	0.130	−0.006	0.052	0.060
$\Delta MS{C}_{2,2}$		0.037	−0.051	−0.070	0.036	0.048
$\Delta MS{C}_{3,2}$		0.070	0.006	0.128	−0.207	0.009
$\Delta MS{C}_{4,2}$		−0.234	−0.059	−0.043	−1.604	2.027

in Appendix C).

Finally, the stiffness proxy derived from accelerations ($KPRX$) is shown to successfully localize damage for the two lower locations of damage (link 01 and 23), as shown in the last row of Figure 10. However, the low difference in signals on the highest floor and the related large variance in the healthy reference case undermines successful localization of damage at the last floor.

As aforementioned, it is suggested that a robust localization scheme would require a combination of multiple DSFs, given the fact that individual DSFs perform well for a specific part of the structure (when relying on the first mode). In order to systematically evaluate the damage localization potential, the results for all 21 GM signals (see Figure 4) are reported in Figure 11. In a similar manner to confusion matrices for machine-learning applications, the true damage location is reported in each line, and the identified damage locations are reported in each column. Perfect localization corresponds to diagonal entries equaling 1 and zero off-diagonal entries. The confusion matrix of the $TAC$ (top right of Figure 11) confirms the previous findings: while the location of damage is correctly picked up, identifying the absence of damage is restricted to the links above the damaged link, and undamaged links cannot be discovered when their position is below the damaged one. Good localization performance is achieved by $WLC$ and $KPRX$ (Figure 11 top center and top right, respectively), despite some spurious detection for the $WLC$ and the incapacity of the $KPRX$ to detect damage at the last floor.

The $TS$ indicator yields satisfactory localization performance (Figure 11 bottom left); however, false positives are observed on the first floor, irrespective of the true damage location, and on the second floor when the real damage is at the third or fourth floor. In addition, damage at link 12 and 45 is not identified in all 21 cases. The $MER$, shown in Figure 11 (bottom center), performs well in picking up the damaged links and, in a similar manner to $TS$, to identify the links above the damaged one as being undamaged. However, as with many features that rely on the modal characteristics of the structure, spurious damage localization in the links (corresponding to floors in buildings) below the damaged one cannot be avoided. Finally, the change in the mode-shape curvature ($\Delta MSC$, shown in Figure 11 bottom right) correctly picks out the location of damage in all cases. Given the formulation of $\Delta MSC$ with respect to DOFs and not links, the specific damaged link cannot be precisely identified. In addition, changes in the MSC at a specific DOF often entail small changes in further locations, thus leading to erroneous damage localization.

In Figure 11, the localization of damage is treated as a Boolean variable (either presence or absence of damage). However, when referring to the value of the computed DSFs and the corresponding probabilistic distributions, further information on the extent of damage (quantification) can be extracted, as real-world structures may sustain damage in more than one floor. For instance, the Mahalonobis distance of the TAC assumes the highest value in vicinity of the damaged link (see first row of Figure 10). A more detailed analysis of the extent of damage per floor would require the use of a structural (physics-based) model, which however lies outside the scope of this paper.

While in Figure 10 and Figure 11, localization refers to cases with a single damage location, the correlation between DSFs evaluated locally for each link and the respective value of the $D{I}_{PA}$ calculated for each link is evaluated in Figure 12 for the original model (defined in Table 2 and Figure 3 in Section 3.1). In general, all DSFs yield a good correlation with damage location (correlation coefficient above 0.80), except for features $CE$ and $TS$, which are formulated with respect to the energy accumulation over time, and the changes in mode-shape curvature ($\Delta MSC$), given the inability to relate damage to a single link. Very good correlation is achieved for $WLC$, $KPRX$, and $MER$. This confirms that while $MER$ has a large scatter in global damage, due to the input energy of GMs, it yields good performance in relative assessment from the same GM, as required for localization.

3.1.4. Tracking DSFs in the Time Domain to Detect Nonlinearity

The detection of damage has been found to yield good results when deriving DSFs from AVs after a damaging earthquake, which has induced a residual shift in the structural stiffness. The performance of OMA-based DSFs thus strongly depends on the stiffness degradation that results from cyclic loading, as indicated by the high correlation between the $PAAC$ and the residual stiffness drop (Figure 7 right). Therefore, the capacity to detect damage is re-evaluated for two new models, for which the stiffness degradation parameter of the Bouc-Wen spring formulation is either reduced from $\delta \eta =3.{610}^{-5}$ to $\delta \eta ={10}^{-5}$ or increased to $\delta \eta ={10}^{-4}$, respectively.

When reducing the stiffness degradation effect, the OMA-based DSFs ($\Delta f_n$, $MAC$, $\Delta MSC$, and $PAAC$) fail to detect damage, given the very low residual effect of damage in this case. Given the reduced effect of nonlinearity on residual damage, DSFs are evaluated during the GM with short data windows of 4 s, with an overlap of $3.5$ s. Thus, temporary nonlinearities may be picked up. The results for one earthquake instance (Montenegro GM with PGA = 5 m/s${}^2$) are shown in Figure 13. The left column contains the simulated results with low stiffness degradation $\delta \eta ={10}^{-5}$, the center column contains the simulated results for high stiffness degradation $\delta \eta ={10}^{-4}$, and the right column contains the results for PGA = 1.5 m/s${}^2$, during which the structure remains in the linear elastic range. Each line contains the Mahalonobis distance (see Equation (17)), with respect to the healthy reference distribution, for five DSFs, which can be formulated for short time windows (see Table 1): $TAC$, $Tcnt$, $MER$, $KPRX$, and $WLC$.

The $TAC$, shown in the first row of Figure 13, clearly identifies the onset of nonlinearity for both models, with and without stiffness degradation, and correctly identifies the residual loss of stiffness, when this occurs (second column of Figure 13), or the absence of a residual effect, when the structure returns to the initial stiffness (first column in Figure 13). Such a capacity to detect temporary stiffness loss and residual stiffness loss is essential when separating failure modes. For instance, unreinforced masonry walls that fail in bending show rocking behavior with minor residual stiffness loss while showing residual stiffness loss when failing due to shear. In a similar manner to the $TAC$, the $Tcnt$ identifies temporary and residual stiffness loss (see Figure 13 second row). Tracking DSFs over time also depends on the window-size, as shorter windows are more noisy and longer time windows less sensitive to temporary nonlinearity (see Appendix B Figure A1). In this case, a window length of 4 s is chosen.

The $MER$, reported in the third row of Figure 13, shows potential to scale with instantaneous stiffness. However, the strong dependence of this DSF on the spectral content of the input motion undermines the detection of damage, as the threshold for detection is not crossed. A similar observation can be made for the $KPRX$ (fourth row of Figure 13), which scales better with the residual stiffness than with the instantaneous stiffness. Finally, the $WLC$, shown in the last row of Figure 13, is very sensitive to temporary stiffness reduction, and it also contains information about the residual stiffness drop. However, this high sensitivity makes it prone to false detection of damage for short time windows, as can be seen for the linear model (last column of Figure 13). The four other DSFs do not detect damage for the linear model.

As stated before, for failure mechanisms with limited stiffness deterioration, the capacity to pick up nonlinearities in short windows is key. The three variations of the toy model, defined by changing the stiffness degradation parameter of the Bouc-Wen spring, are used to assess the correlation of DSFs with respect to the Park and Ang index ($D{I}_{PA}$) (see Figure 14). Two DSFs that are derived from post-earthquake AVs, $\Delta f$ (Figure 14, top left) and $TAC$ (Figure 14, top center), are very sensitive to stiffness degradation: for the model with very weak stiffness deterioration ($\delta \eta ={10}^{-5}$), the residual damage is small, and thus, the nonlinearity is not detected. In addition, the high performance for OMA-based DSFs to detect and quantify nonlinearity, shown in Figure 8, would significantly drop when considering multiple models with changing hysteretic properties, as evidenced by the large scatter obtained for the change in natural frequency, $\Delta f_1$, in Figure 14 (top left).

When using DSFs that can be split into short moving time windows, residual damage and nonlinearity can be addressed separately (see Figure 13). To detect (peak) nonlinearity, the 90-percentile value of the DSF values derived for the windows defining the GM (see Figure 2 bottom) may be used. The 90-tile values of $TAC$ (Figure 14, top right) and $WLC$ (Figure 14, bottom left), scale well with the $D{I}_{PA}$, irrespective of the stiffness degradation parameter. This feature is important to detect nonlinearity that may weaken a structure but not result in residual changes.

Finally, two DSFs that pinpoint nonlinearity from the entire GM signal are evaluated: $CE$ (Figure 14, bottom center) and $CohAC$ (Figure 14, bottom right). While both features scale with $D{I}_{P}A$, the lack of differentiation between the three models reveals their dependence on nonlinearity alone and the related incapacity to detect residual damage. In addition, the CohAC fails to pick up nonlinearity for GM signals that correspond to short earthquakes, such as Northridge, Marche, and Amatrice (see Figure 4).

3.1.5. Sensitivity to Measurement Noise

Permanent sensor installation in multiple buildings requires the use of affordable sensors, which come with the drawback of reduced precision and higher noise floors. In Figure 15, the influence of changing noise levels on the variability of OMA-based DSFs and the related differentiation between linear and nonlinear model responses is assessed. As observed for all three DSFs ($\Delta MSC$, $PAAC$, and $TAC$), a signal-to-noise ratio of at least 5 is required to avoid large variability in the results. Given that the amplitude of AVs typically ranges between ${10}^{-3}$ and ${10}^{-1}{\mathrm{mm}/\mathrm{s}}^2$, this requires sensors with resolution around ${10}^{-4}$${\mathrm{mm}/\mathrm{s}}^2$. For DSFs derived from GMs, given the higher amplitude of shaking, lack of precision only initiates at noise levels of $5{\mathrm{cm}/\mathrm{s}}^2$, which is above the typical range of modern—even low-cost—sensors.

The DSFs evaluate structural behavior through an input–output sensor pair. For such a setup to be successful, good synchronization needs to be guaranteed between sensor nodes.

3.2. Shake-Table Data

The performance of DSFs for data-driven identification of nonlinearity/damage is further assessed on actual vibration-based monitoring data that has been recorded during shake-table tests performed by Beyer et al. [90] with a four-story mixed reinforced-concrete unreinforced-masonry building on the EUcentre shake table in Pavia. The building specimen has been tested with the Montenegro earthquake, and scaled to increasing amplitudes. The last test has led to failure (unstable building) and is therefore discarded, as SHM is not required for flagging collapsed buildings. Structural damage accumulated during the first eight earthquake, as reported in

Table 3. Description of the eight earthquake sequences and the evolution of damage observed by Beyer et al. [90].

Test #	PGA [g]	D.G.	Damage Description
1	0.05	1	Single hairline crack one URM wall of first floor.
2	0.1	1	Haircracks in one wall of first two floors and in the construction joint between wall and foundation
3	0.2	1	Cracks in all masonry walls of first two floors
4	0.3	2	Several diagonal cracks over the entire wall height of one wall of the 1st floor with negligible residual crack width. Many flexural cracks in the concrete slab of the first floor; masonry spandrels and concrete wall remained undamaged.
5	0.4	2	Same as previous.
6	0.6	2	Significant increase in damage to the structure. All masonry walls with diagonal cracks at all floors. First and second floor walls present residual cracks of 0.8 mm.
7	0.4	2	This test, with a smaller amplitude than the previous EQK6, was intended to simulate a possible aftershock but led only to very minor additional dam-age to the structure.
8	0.7	3	Structure was severely damaged. Damage in the masonry walls started concentrating in one diagonal crack. Diagonal cracks passed through bricks.

, and EMS-98 damage grades (DGs), ranging from 1 to 3, have been attributed based on visual inspection of the specimen. The building has been excited with uni-directional WN signals after each earthquake (except after Test # 7). More details about the test protocol and the building geometry can be found in [90].

The DSFs are derived from recorded accelerations, with one sensor per floor considered for the analysis. However, the base shear is also recorded by the table, and the top displacement has been measured [90]. While displacements and, above all, reaction forces are challenging to measure in a large population of existing buildings, these quantities are used to validate the capacity of DSFs to detect nonlinearity by computing the peak roof drift ratio and the hysteretic work dissipated by the structure. The residual damage is mostly characterized by the damage description and the EMS-98 DG (see Table 3).

The results for five DSFs, which have been found to perform well on the simulated case study, are represented in Figure 16. The DSFs are calculated with respect to two different references: the WN excitation corresponding to the healthy reference state, recorded prior to the first earthquake test, and the WN excitation performed prior to the respective earthquake (labeled $WN-1$ in Figure 16). Only for Test #7, no pre-earthquake WN data were available, and therefore, the equivalent data prior to test #6 have been used. While deriving DSFs with respect to the healthy reference state corresponds to tracking the absolute damage of the structure, using the WN before the earthquake provides an estimate of the additional damage that is sustained by the structure.

Generally, the values of the DSFs are in good agreement with the nonlinearity indicators (hysteretic work on the left of Figure 16) and peak roof drift ratio on the right of Figure 16) and the damage description (see Table 3). The relative TAC, computed with respect to the pre-earthquake WN) and shown in the first row of Figure 16, has the highest increase for tests #7 and #9, after which significantly increased structural damage states have been noted (see Table 3). The TAC value does not increase between tests #7 and #8, indicating that the aftershock did not further damage the structure. The capacity of data-driven indicators to accelerate post-earthquake assessment by identifying buildings that do not sustain additional damage could be a major contributor to accelerated post-earthquake recovery. In addition, the $TAC$, when calculated in an absolute manner, provides a clear separation between non-damaging earthquakes (#1 and #2) and damaging earthquakes as well as between slight and extensive damage. This is a necessary starting point for successful building tagging based on a traffic-light logic. In addition, when compared to visually observed DGs (EMS98 damage scale), the TAC provides a more refined information about damage damage, which could be of interest for subsequent re-pair design.

While a popular feature and producing good results with the simulated case study, the change in natural frequency ($\Delta f_1$, shown in the second row of Figure 16) is sensitive to damage but does not scale well with the observed severity of damage, as the relative decrease in this DSF is higher for the minor damage sustained in the first four tests than for the major structural damage sustained in tests #6 and #8. Given no WN tests have been conducted following the aftershock, this DSF cannot be evaluated for test #7.

The stiffness proxy, $KPRX$, is reported in the third row of Figure 16) and when calculated with respect to the pre-earthquake WN contains useful information to separate small damage ($DG1$) from more severe damage ($DG\ge 2$). However, the $KPRX$ does not scale well with the nonlinearity indicators and is not able to track the absolute damage with respect to the initial “healthy” building state, as the frequency range (m in Equation (11)) is not updated and thus, the $KPRX$ becomes noisy.

The change in the phase-angle of the cross power-spectral density, $PAAC$ (see fourth row of Figure 16), derived from post-earthquake WN excitation, is very sensitive to damage, yet it does not scale well with the two nonlinearity indicators (hysteretic work and peak roof drift ratio). While the relative formulation of this DSF is a good source to detect the appearance of new cracks and significant damage, it fails to detect the growth of cracks, such as those caused by earthquakes #4 and #5).

In a similar manner, the time-domain correlation of wavelet coefficients ($WLC$) provides a good separation between minor damage and significant structural damage (from test #6). In addition, the increment in damage after the sixth test was dominated by the peak roof drift ratio, while in the previous tests, the hysteretic work dominated the $PAindex$, as shown in Figure A2.

The application of DSFs on the shake-table data confirms the observation that they provide complementary information and thus, a combination of several DSFs would be optimal. For instance, the $TAC$ contains information about the overall damage state, while $PAAC$ provides useful information about increasing the extent of damage (appearance of new cracks). Furthermore, $WLC$ shows a good separation between minor damage and severe nonlinearity due to displacement beyond the elastic, mostly linear, range.

Given the observation that DSFs may track the absolute and the incremental damage sustained by structures, the $TAC$ and the $WLC$, which can be both calculated for shorter time windows (see Table 1), are evaluated in greater detail. The tracking of both DSFs with short time windows over the eight earthquake instances is reported in Figure 17: while the damage steadily increases during the first four instances, no increase is observed in the absolute $TAC$ values (Figure 17 top) during the fifth and the seventh signal. In addition, when comparing the absolute and the incremental formulation of both damage sensitive features, it is observed that the absolute value saturates after reaching a value around 0.8, which underlines the importance to track the damage with incremental values as well as the largest increment in DSF is found for the last—damaging—earthquake. In general, the TAC scales better with the amount of damage, which is due to the fact that it covers a wider frequency bandwidth. Finally, for many test instances, it can be observed that the large-amplitude shaking (at the beginning of the Montenegro earthquake, see Figure 4) leads to higher values of the DSF when compared with the final values of the DSF when the shaking has a lower amplitude. This underlines the capacity to separate the nonlinearity from the residual damage.

Formalizing the observations of Figure 17 into a quantitative assessment of the change using the KL divergence see (Equation (18)), the absolute damage can be tracked as well as the incremental damage with respect to the previous WN excitation and also the previous GM data, as shown in Figure 18. When using the previous GM as reference, earthquake signal #6, which has led to a significant increase in damage (see Table 3) clearly stands out. In addition, earthquakes #3 (appearance of cracks in the RC walls, thus significantly reducing the stiffness) and #8 (severe damage) correspond to significant increases in damage. The capacity to track incremental damage is crucial for safe, yet efficient, functionality of buildings during earthquake sequences.

3.3. Limitations

Both the toy model presented in Section 3.1 and the shake-table structure (Section 3.2) are rather ductile and thus, nonlinearity and damage forms without immediate failure. Therefore, the capacity of DSFs to diagnose brittle failures is not addressed in this paper. However, failure states do not necessary require monitoring data to be diagnosed. In addition, the application to very stiff and very flexible buildings is left for future work and thus, possible influences of signal propagation effects in high-rise buildings remain to be evaluated. In addition, the influence of strong redundancy in the load-bearing system remains to be investigated.

4. Conclusions

This work assesses data-driven damage-sensitive features (DSFs) in terms of their potential to identify earthquake-induced damage and manifested nonlinearity in the response of low-to-midrise residential buildings. Based on a thorough simulated analysis on an exemplary nonlinear hysteretic spring-mass system and a half-scale building specimen experimentally tested on a shake table, the following conclusions are drawn:

• Damage-sensitive features derived in the frequency domain, using either transmissibility or wavelet decomposition, are promising for data-driven nonlinearity detection. However, reference linear data are required to detect the onset of nonlinearity.
• Damage-sensitive features do not only indicate presence of damage and evolution of nonlinearity but further carry the potential to localize damage and quantify its extent.
• Damage-sensitive features allow to identify the evolution of the structural state and the accumulation of damage over the course of multiple earthquake instances, in addition to picking up the severity of damage that is incrementally introduced by a specific earthquake in a sequence.
• DSFs calculated for short time windows are necessary to track nonlinearity over time, which enables differentiating between reversible and irreversible stiffness drops.

The findings of this thorough comparative assessment of DSFs further show that an exclusive reliance on WN reference data may lead to false alarms, as the frequency content of earthquake signals may lead to DSFs that differ from the reference values even for linear responses. In addition, WN data recorded after an earthquake event can only serve for offering information on (possible) residual damage. Instead, the continuous monitoring and appropriate windowing (segmentation) of the vibration response, which is recorded during strong ground motions, is key to retrieving additional information on temporarily experienced nonlinearity. The latter, as well, may weaken the structure and consume part of its capacity to absorb hysteretic energy. Finally, the formulations of DSFs and their target domain differ and thus, robust damage diagnosis should involve a combination of multiple DSFs.

Future work will focus on the fusion of individual DSFs in a classification framework, which can more precisely localize damage and assess its influence on the residual seismic capacity. Moreover, the reliability of purely data-driven assessment will be compared against model-based approaches, which are required to diagnose the consequence of damage on the seismic performance.