# Towards robust statistical damage localization via model-based sensitivity clustering

Saeid Allahdadian, Michael Döhler, Carlos Ventura, Laurent Mevel

► To cite this version:

Saeid Allahdadian, Michael Döhler, Carlos Ventura, Laurent Mevel. Towards robust statistical damage localization via model-based sensitivity clustering. Mechanical Systems and Signal Processing, 2019, 134, pp.106341. 10.1016/j.ymssp.2019.106341. hal-02293021

HAL Id: hal-02293021

https://inria.hal.science/hal-02293021

Submitted on 20 Sep 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Towards robust statistical damage localization via model-based sensitivity clustering

Saeid Allahdadianª, Michael Döhlerᵇ,*, Carlos Venturaª, Laurent Mevelᵇ

ªThe University of British Columbia, Department of Civil Engineering, Vancouver, B.C., V6T 1Z4, Canada

ᵇUniv. Rennes, Inria, Ifsttar COSYS/SII, Campus de Beaulieu, 35042 Rennes, France

## Abstract

Damage diagnosis is a fundamental task for structural health monitoring (SHM). With the statistical sensitivity-based damage localization approach, a residual vector is computed from vibration measurements in the reference and the damaged state. The residual is analyzed statistically in hypothesis tests with respect to change directions defined by the sensitivities of the structural parameters associated to elements of a finite element (FE) model of the investigated structure. If the test for a parameter reacts, then the respective element of the structure is indicated as damaged. This approach offers a very generic and theoretically sound framework to analyze parametric changes in systems, and takes into account the intrinsic statistical uncertainty related to measurement data. Depending on the definition of the residual and of the parameterization, the approach offers a simple computation of the test statistics directly from the measurement data in the damaged system, without the need of system identification. Since an FE model is used, the approach is applicable on arbitrary structures, while no model updating is required and therefore the requirements on the FE model accuracy are less strict. While the theoretical framework has been developed previously, it lacked robustness so far for an application on real structures. The purpose of this paper is the development of this framework into a working damage localization method that is applicable on real data from complex structures. To achieve this goal, robust hypothesis tests are used, the sensitivity computation of the residual is revisited for more precision thanks to reduced modal truncation errors, and an adequate clustering approach is proposed for the case of a high-dimensional FE parameterization for complex structures. Furthermore, several robustness properties of the method are proven. Finally, an application of this framework is shown for the first time on experimental data for damage localization, namely in an ambient vibration test of a 3D steel frame at the University of British Columbia.

**Keywords:** Damage localization, Structural vibration monitoring, Ambient excitation, Hypothesis testing, Subspace methods

## 1. Introduction

The diagnosis of damages is a fundamental task for structural health monitoring (SHM) [1]. With the advent of new technologies, instrumentation of structures is becoming widespread. In particular, vibration-based techniques have been actively developed in the last decades [2–4], for example for the monitoring of bridges, buildings or offshore structures [5]. Physical changes in the structure due to damage induce changes in the modal characteristics of the structure, which can be monitored through output-only vibration measurements.

Damage diagnosis is usually divided into five sub-tasks of increasing difficulty [1]: damage detection (level 1), damage localization (level 2), identification of the damage type (level 3), quantification of the damage extent (level 4) and prediction of the remaining service life (level 5). Methods for damage detection are the most developed since they can operate purely data-based and do not require a physical model of the monitored structure. Examples are algorithms from pattern classification and statistical process control, e.g. [6–8].

Methods for damage localization are more sophisticated since some link between the measurement data and the physical properties of the structure is required in order to localize the physical change in the structure due to damage.

<footnote>*Corresponding author; E-mail address: michael.doehler@inria.fr</footnote>
Such a link is often given by a finite element (FE) model or by directly assuming specific structural types. Data-driven damage localization methods [3, 4] are usually designed for specific structural types, like beams, plates or rotating machinery, and often require dense sensor grids. However, they are not easily generalizable to arbitrary structural types. Using an FE model usually leads to more universal damage localization approaches which are well-adapted for complex structures. E.g., in model-based methods the parameters of an FE model of the healthy structure are updated based on measurements from the damaged system, and damage is localized in the regions of the model where parameters have changed [9–11]. A shortcoming of model-based approaches is often poor conditioning due to the huge FE parameter dimension in comparison to relatively few modal parameters that can be extracted from data, which requires user interaction by an experienced engineer in the updating process [9, 10]. In these approaches, the damage extent quantification is not decoupled from the damage localization, which also contributes to poor conditioning due to the demand of both information at the same time.

Alternative damage localization approaches, which also use an FE model of the investigated structure, intend to achieve damage localization (level 2) without addressing the more difficult quantification task (level 4) [12]. Therefore, by addressing a simpler task, their requirements on the accuracy of the FE model are usually less strict than for updating approaches. Usually, these approaches are strongly based on data-driven features from measurements of the reference and damaged states, as well as on some information from an FE model. Damage indicators are defined with respect to the elements of the FE model without updating it, combining properties of data-driven and model-based approaches. For example, damage is localized by interrogating changes in the flexibility or in the transfer matrix of a structure by applying certain load vectors (computed from the data) to an FE model in [13–16]. In [17], damage is located by analyzing stress fields on the boundaries of closed regions of the structure that are computed from measurement data and the FE model.

Also belonging to this class of combined data-driven and model-based approaches, the statistical subspace-based damage localization approach introduced in [18, 19] computes a residual vector from the measurements that is based on subspace properties of a Hankel matrix containing output covariances. The residual is analyzed statistically in hypothesis tests with respect to change directions defined by sensitivities of the structural parameters associated to elements of the FE model. If the test for a parameter reacts, then the respective element of the structure is indicated as damaged, so damage is localized in that element. Indeed, the underlying asymptotic local approach to change detection [20] on which this work is built, combined with adequate hypothesis tests [21], offers a very generic and theoretically sound framework to analyze parametric changes in systems. Depending on the definition of the residual and of the parameterization, the approach offers a simple computation of the test statistics directly from the measurement data in the damaged system, without the need of system identification, for a decision of which system parameters have changed and thus where the damage is localized. Therefore, the questions of damage localization and quantification are decoupled, and the FE updating problem is avoided.

While the underlying framework has been thoroughly developed for statistical subspace-based damage detection [18, 22–25] with successful applications on structures in the lab and in the field [26–29], the respective damage localization approach introduced in [18, 19] has only been applied in simulations so far despite its appealing theoretical properties and its generality. In this paper, we revisit the approach with the purpose of developing it into a working method applicable on real ambient vibration data from complex structures. To achieve this goal, the framework for the hypothesis tests of [21] is taken into account, and developments for the appropriate sensitivity computation and clustering of high-dimensional FE parameters are required. Besides laying out the details of the framework and proving several robustness properties of the method, the principal contributions of this paper to achieve a working damage localization method are the following with respect to previous work [18, 19]:

* Choice of the hypothesis tests: While previously so-called sensitivity or direct tests where used [18, 19] that only take into account the sensitivity of the currently tested parameter, we propose to use minmax tests, where the sensitivities of all parameters are taken into account each time. We show that it leads to more robust localization results by avoiding false alarms in healthy elements.

* Sensitivity computation: The sensitivity of the residual with respect to the structural parameters is a key element for the damage localization approach. It is computed in two steps: first, the sensitivity of the residual with respect to the modal parameters is evaluated, and second the sensitivity of the modal parameters with respect to the structural parameters is computed. Previously, the first part has been computed using the identified modal parameters from measurements in the reference state [18, 19]. Instead, we propose to compute this sensitivity
using the modal parameters from the FE model. We show that in spite of possible model errors, the model-based approach may reduce modal truncation errors of the sensitivity significantly, which leads to more robust localization results.

* **Clustering:** Depending on the physical parameterization and the number of parameters, there may be parameters that are *close* in the sense that a change therein may have nearly the same effect on the residual. In this case, a preliminary clustering step is required. Previously, *k*-means clustering [30] was proposed [18, 19], which however depends strongly on a random initialization and does not necessarily cluster close elements satisfactorily. We propose a hierarchical clustering strategy, namely complete-linkage clustering [30], and show that it is well adapted to our damage localization approach.

With these developments, our localization framework is well adapted to meet the requirements of real applications. This is shown by its first application to experimental data in this paper, namely in an ambient vibration test of the Yellow Frame [31], which is a 3D steel frame at the University of British Columbia.

This paper is organized as follows. In Section 2 the basic vibration models and parameterizations are introduced. In Section 3, the subspace residual previously used in [18, 19] and other residuals [25, 32] are presented that are suitable for damage localization in our framework. In Section 4, the details of the damage localization method are developed and illustrated on a numerical example. Its application to experimental data of the Yellow Frame is shown in Section 5, before concluding the paper in Section 6.

## 2. Vibration model and parameters

The considered damage localization method connects vibration measurements to structural system parameters to evaluate changes in those parameters statistically. In this section, the underlying vibration models and parameterizations are recalled.

### 2.1. Vibration model

The behavior of mechanical structures subject to unknown ambient excitation is assumed to be described by the differential equation

$$ M\ddot{z}(t) + C\dot{z}(t) + Kz(t) = v(t) \quad (1) $$

where $t$ denotes continuous time, $M, C, K \in \mathbb{R}^{m \times m}$ are mass, damping, and stiffness matrices, respectively, the vector $z(t) \in \mathbb{R}^m$ collects the displacements of the $m$ degrees of freedom (DOF) of the structure, and $v(t)$ is the external unmeasured force (random disturbance). Observing the displacements, velocities or accelerations of system (1) yields the measurement vector

$$ y(t) = L_d z(t) + L_v \dot{z}(t) + L_a \ddot{z}(t) + e(t), \quad (2) $$

where the selection matrices $L_d, L_v, L_a \in \{0, 1\}^{r \times m}$ indicate the observed displacements, velocities or accelerations at the respective DOFs, and $e(t)$ denotes the measurement noise.

Sampling at discrete time instants $t = k\tau$ (with sampling rate $1/\tau$), system (1)-(2) can also be described by a discrete-time state space system model [33]

$$ \begin{cases} x_{k+1} = A x_k + w_k \\ y_k = C x_k + v_k \end{cases} \quad (3) $$

with the state vector $x_k = [z(k\tau)^T \quad \dot{z}(k\tau)^T]^T \in \mathbb{R}^n$, $n = 2m$, the measured output vector $y_k = y(k\tau) \in \mathbb{R}^r$ and the system matrices

$$ A = \exp\left(\left[-\mathcal{M}^{-1}\mathcal{K} \quad -\mathcal{M}^{-1}C\right]\tau\right) \in \mathbb{R}^{n \times n}, \quad C = \left[L_d - L_a\mathcal{M}^{-1}\mathcal{K} \quad L_v - L_a\mathcal{M}^{-1}C\right] \in \mathbb{R}^{r \times n}. \quad (4) $$

The state noise $w_k$ is related to the unmeasured force $v$, and the output noise $v_k$ is related to measurement noise $e$ as well as to the unmeasured force $v$ in case of acceleration measurements. Both noise terms are unmeasured and assumed to be zero-mean and white, while in practice also colored noise may be feasible [34].
## 2.2. System parameterization for damage localization

Damage in the system is related to a change in the physical properties of the structural elements. A basic step in our damage localization methodology is the definition of a damage-sensitive system parameter vector. Let such a parameter vector be called $\theta \in \mathbb{R}^l$, where $l$ is the number of parameters, and let its value in the (undamaged) reference state be $\theta_0$. The system parameterization is user-defined and adapted to the localization problem at hand. It should be defined such that damage in any structural element (or in groups of structural elements) of the system corresponds to a change in one entry of parameter vector $\theta_0$. The problem of damage localization is then to find out which entries of the parameter vector are responsible for the change in the system based on measurements from the healthy and from the damaged system. The parameter vector consists of, e.g., stiffness, cross sections or mass parameters of the different structural elements, or any other physical parameter that shall be monitored and whose change affects the modal parameters of the system. The parameter vector is usually defined in connection with an FE model that approximates the monitored structure described by (1), where $\mathcal{K}, \mathcal{M}$ and/or $\mathcal{C}$ are functions of $\theta$. Note that damage localization requires the parameter vector to be explicitly linked to the physical properties of the different structural elements, while damage detection can be performed more simply by detecting changes in the modal characteristics of the system [18, 19, 21].

# 3. Residual vectors – damage features

Any damage diagnosis method requires the extraction of damage-sensitive features from the measurement data of the monitored system [1]. With our damage localization approach, this feature vector is defined in a way that it is approximately Gaussian distributed with zero mean in the reference state and non-zero mean in the damaged state, hence the designation residual vector. In this section, residual vectors are presented that fit into this framework and that are suitable for the subsequent damage localization approach.

## 3.1. Subspace residual

In covariance-driven subspace identification a block Hankel matrix $\mathcal{H}$ is built from the output covariances $R_i = E(y_{k+i}y_k^T) \in \mathbb{R}^{r \times r}$, whose column space contains the information on the system matrices $(A, C)$ and thus the modal characteristics [35–37]. This property is used to define a residual vector from the measurements without doing system identification [18, 21–25] as follows. The left null space matrix $S(\theta_0)$ of $\mathcal{H}$ in the reference state is computed, which satisfies $S(\theta_0)^T \mathcal{H} = 0$ if and only if $\mathcal{H}$ is obtained from the system in the reference state, i.e. $\theta = \theta_0$. To check this property using measured data $\mathcal{Y}_N = \{y_1, y_2, \dots, y_N\}$ of the system in the current state under unknown parameter $\theta$, a consistent estimate $\hat{\mathcal{H}}$ is obtained from the output covariance estimates $\hat{R}_i = \frac{1}{N} \sum_{k=1}^N y_{k+i}y_k^T$, and the residual vector

$$ \zeta_1(\theta_0, \mathcal{Y}_N) = \sqrt{N} \operatorname{vec}(S(\theta_0)^T \hat{\mathcal{H}}) \quad (5) $$

is defined, where $\operatorname{vec}(\cdot)$ denotes the column stacking vectorization operator. Note that an estimate of $S(\theta_0)$ is easily obtained from a singular value decomposition (SVD) of $\hat{\mathcal{H}}$ in the reference state.

## 3.2. Subspace residual robust to excitation changes

Changes in the ambient excitation properties do not affect the dynamic properties of the monitored structure, but they may have an influence on the Hankel matrix $\hat{\mathcal{H}}$ and thus on (5). The robust subspace residual from [25] is based on the same subspace properties as in the previous section, but it is not affected by changes in the ambient excitation properties. There, the Hankel matrix estimate $\hat{\mathcal{H}}$ computed on the test data is replaced by the left singular vector matrix $\hat{U}_1 \in \mathbb{R}^{(p+1)r \times n}$ from an SVD of $\hat{\mathcal{H}}$, leading to definition

$$ \zeta_2(\theta_0, \mathcal{Y}_N) = \sqrt{N} \operatorname{vec}(S(\theta_0)^T \hat{U}_1). \quad (6) $$
### 3.3. Transfer matrix-based residual

A residual that is not built on subspace properties, but on the transfer matrix difference between reference and damaged states, was introduced in [32]. It requires the estimation of the system matrices $C$ and $A_c$ of the continuous-time state space system corresponding to (3) on the test data, following [14, 16]. Define

$$R(s) = C(sI - A_c)^{-1} \begin{bmatrix} CA_c \\ C \end{bmatrix}^{\dagger} \begin{bmatrix} I \\ 0 \end{bmatrix},$$

where $s$ is a Laplace parameter in the complex plane. This matrix is computed in the reference state, denoted by $R^{\theta_0}(s)$, and let its estimate from the test data $\mathcal{Y}_N$ be $\hat{R}(s)$. The respective residual vector is

$$\zeta_3(\theta_0, \mathcal{Y}_N) = \sqrt{N} \mathrm{vec}(\hat{R}(s) - R^{\theta_0}(s))_{\mathrm{re}}, \quad (7)$$

where $(\cdot)_{\mathrm{re}}$ denotes the stacking of the real and the imaginary part of the vector.

### 3.4. Further residuals

While a detailed statistical analysis of the residuals from the previous sections has been made in the respective works, many other residuals have been used in the literature without such an analysis, but which would also be compatible with the following damage localization approach after slight adaptations. To name a few, this concerns e.g. residuals built on a null-space based comparison of data Hankel matrices [38], on the difference of output covariance Hankel matrices [39] or directly on the modal parameter differences [40, 41].

## 4. Statistical sensitivity-based damage localization method

The residuals defined in the previous section are computed from the measurement data, and are thus random variables. It can be shown that they are asymptotically Gaussian distributed, and their asymptotic mean changes when damage appears. Linking this mean to the system parameterization by the respective sensitivity, and taking into account the residual covariance, damage is localized by using statistical tests for changes in each of the system parameters (see also Section 2.2).

In this section, the basics of these statistical tests are first laid out in a generic Gaussian framework, where all aforementioned residuals fit in. Second, the residual sensitivity, which plays an essential role in these tests, is elaborated in detail. Two different practical ways of its computation are developed. Third, a tailor-made clustering approach is developed based on the residual sensitivity and covariance. This ensures the robustness of the damage localization approach in particular when the number of system parameters is high in comparison to the number of sensors.

**Illustrative example.** In the course of these developments, a simple numerical example will be used for illustration purposes, namely a shear wall model in two configurations containing six or 15 elements shown in Figure 1. The stiffness and the mass are equal for all levels. The damping ratio is set to 2% in all modes. Vibration data at three sensors (see Figure 1) is simulated from Gaussian white noise excitation in all story levels. Gaussian measurement noise with 5% magnitude of the signals is added. Damage is simulated by reducing stiffness in one or some of the walls. In all cases, the subspace residual from Section 3.1 is applied.

### 4.1. Hypotheses definition and Gaussian residual vector

The residual vector is analyzed statistically for changes in the parameter vector in hypothesis tests. It is analyzed if $\theta = \theta_0$ or if $\theta \neq \theta_0$ for each entry of the parameter vector, and it needs to be decided which entries in $\theta_0$ have actually changed for damage localization. For mathematical convenience, the underlying hypotheses are formulated based on the asymptotic local approach for change detection [20] as

$$\mathbf{H}_0: \theta = \theta_0 \quad (\text{reference system}), \\ \mathbf{H}_1: \theta = \theta_0 + \delta / \sqrt{N} \quad (\text{damaged system}), \qquad (8)$$

where $\delta \in \mathbb{R}^l$ is an unknown but fixed change vector. With this statistical framework, very small changes in the system parameter $\theta_0$ can be detected if $N$ is large enough. Moreover, this framework allows to characterize the asymptotic
<figure><img src="image_1.png" /><figcaption>Figure 1: Shear wall model in two configurations with sensor positions, composed of (a) six elements, and (b) 15 elements.</figcaption></figure>

probability distributions of the residuals defined in Section 3. Indeed, it can be shown that all those residuals are asymptotically Gaussian distributed for $N \to \infty$ with [21, 22]

$$ \zeta(\theta_0, \mathbf{y}_N) \xrightarrow{d} \begin{cases} \mathcal{N}(0, \Sigma) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\mathcal{J}\delta, \Sigma) & \text{under } \mathbf{H}_1, \end{cases} \quad (9) $$

where $\mathcal{J}$ and $\Sigma$ are the asymptotic sensitivity and covariance of the respective residuals. Besides the residuals defined in Section 3, any other residual can be treated with the subsequent localization approach if it satisfies property (9). The computation of the residual sensitivity and covariance depends on the particular residual and general guidelines are given in Section 4.3 and 4.4. Thanks to property (9), the computed residual can be assumed to be approximately Gaussian distributed when the number of samples $N$ is large enough, having zero mean in the reference state and non-zero mean in the damaged state, and the same covariance in both states. Note that the factor $\sqrt{N}$ in the residuals in Section 3 is essential for convergence to a Gaussian distributed random variable. In this Gaussian framework, classical statistical hypothesis tests can be applied for damage diagnosis as detailed in the following section.

## 4.2. Statistical tests

For damage localization, it has to be decided which entries in the parameter vector $\theta_0$ have changed, i.e. which entries in the change vector $\delta = [\delta_1 \delta_2 \dots \delta_l]^T$ are non-zero. This is done by testing the entries $\delta_i$, $i = 1, \dots, l$ corresponding to each parameter separately, and damage localization in element $i$ corresponds to testing $\delta_i = 0$ against $\delta_i \neq 0$. For each $\delta_i$ that is tested, define also the change vector corresponding to the non-tested elements as

$$ \delta_{\bar{i}} = [\delta_1 \dots \delta_{i-1} \delta_{i+1} \dots \delta_l]^T. $$

Let the respective columns of the sensitivity matrix $\mathcal{J} = [\mathcal{J}_1 \mathcal{J}_2 \dots \mathcal{J}_l]$ be selected accordingly as $\mathcal{J}_i$ and

$$ \mathcal{J}_{\bar{i}} = [\mathcal{J}_1 \dots \mathcal{J}_{i-1} \mathcal{J}_{i+1} \dots \mathcal{J}_l], \quad (10) $$

and let the respective parts of the Fisher information matrix $F = \mathcal{J}^T \Sigma^{-1} \mathcal{J}$ be rearranged as

$$ \begin{bmatrix} F_{i,i} & F_{i,\bar{i}} \\ F_{\bar{i},i} & F_{\bar{i},\bar{i}} \end{bmatrix} = \begin{bmatrix} \mathcal{J}_i^T \Sigma^{-1} \mathcal{J}_i & \mathcal{J}_i^T \Sigma^{-1} \mathcal{J}_{\bar{i}} \\ \mathcal{J}_{\bar{i}}^T \Sigma^{-1} \mathcal{J}_i & \mathcal{J}_{\bar{i}}^T \Sigma^{-1} \mathcal{J}_{\bar{i}} \end{bmatrix}. \quad (11) $$
In the following subsections, two tests for parameter subsets are recalled from [21, 42].

### 4.2.1. Direct test

The simplest possibility for testing $\delta_i = 0$ against $\delta_i \neq 0$ is to assume no change in the remaining parameters, i.e. $\delta_{\bar{i}} = 0$. Then, the generalized likelihood ratio (GLR) test amounts to the test statistic

$$t_{\text{dir}}^{i} = \zeta^{T} \Sigma^{-1} \mathcal{J}_{i} (\mathcal{J}_{i}^{T} \Sigma^{-1} \mathcal{J}_{i})^{-1} \mathcal{J}_{i}^{T} \Sigma^{-1} \zeta, \quad (12)$$

This test is called sensitivity test or direct test. The test variable $t_{\text{dir}}^i$ is $\chi^2$-distributed with one degree of freedom and non-centrality parameter $F_{i,i} \delta_i^2$ if $\delta_{\bar{i}} = 0$ is actually true. If the assumption $\delta_{\bar{i}} = 0$ does not hold, the non-centrality parameter of $t_{\text{dir}}^i$ is $(F_{i,i}^{1/2} \delta_i + F_{i,i}^{-1/2} F_{i,\bar{i}} \delta_{\bar{i}})^2$ [21].

With this test, damage is located in the elements whose respective test values $t_{\text{dir}}^i$ are the highest for $i = 1, \dots, l$. This test was used in previous related works on damage localization [18, 19]. However, the underlying assumption $\delta_{\bar{i}} = 0$ is not true when more than one element is damaged, or when testing undamaged elements while others are damaged. In this case, the non-centrality parameter of the respective tests depends on the relation of the tested parameter to the remaining parameters in $F_{i,\bar{i}}$ and of the change $\delta_{\bar{i}}$. Thus, the results from this test may be unreliable in these situations.

### 4.2.2. Minmax test

Instead of assuming $\delta_{\bar{i}} = 0$, the minmax test is designed such that it reacts only to a change in the tested element, while being blind to changes in the other elements [21, 42]. Assume that $\mathcal{J}$ is of full column rank. First, define the partial residuals $\zeta_i \stackrel{\text{def}}{=} \mathcal{J}_i^T \Sigma^{-1} \zeta$ and $\zeta_{\bar{i}} \stackrel{\text{def}}{=} \mathcal{J}_{\bar{i}}^T \Sigma^{-1} \zeta$, the robust residual $\zeta_i^* \stackrel{\text{def}}{=} \zeta_i - F_{i,\bar{i}} F_{\bar{i},\bar{i}}^{-1} \zeta_{\bar{i}}$ and $F_i^* \stackrel{\text{def}}{=} F_{i,i} - F_{i,\bar{i}} F_{\bar{i},\bar{i}}^{-1} F_{\bar{i},i}$. Then, the mean of the robust residual $\zeta_i^*$ is only sensitive to changes $\delta_i$, but blind to $\delta_{\bar{i}}$, and the corresponding GLR test statistic for $\delta_i = 0$ against $\delta_i \neq 0$ turns out to be

$$t_{\text{mm}}^{i} = \zeta_{i}^{*T} F_{i}^{*-1} \zeta_{i}^{*}, \quad (13)$$

which is called minmax test. The test variable $t_{\text{mm}}$ is $\chi^2$-distributed with one degree of freedom and non-centrality parameter $\delta_i^T F_i^* \delta_i$, independently of $\delta_{\bar{i}}$. For a decision, the test variable should be compared to a threshold. Damage is located in the elements $i$ with the highest values of the test variable $t_{\text{mm}}^i$. A numerically efficient computation of the test variable is described in detail in [21].

### 4.2.3. Comparison of direct and minmax test

If there is no damage in the tested element, i.e. $\delta_i = 0$, the non-centrality parameter of the minmax test is always zero, while this is not necessarily the case in the direct test. This may lead to false positives in the direct test. On the other side, if there is actually damage in the tested element, i.e. $\delta_i \neq 0$, but also damage in other elements, it may happen that the direct test actually reacts less (or even not at all) in the case $F_{i,\bar{i}} \delta_{\bar{i}} < 0$, leading to a false negative result, while the minmax test reacts. Thus, the minmax test should be preferred to the direct test, while in related previous works on damage localization [18, 19] only the direct test was considered.

*Illustrative example.* The direct and the minmax test are applied to the shear wall example containing six elements in Figure 2, where damage is introduced in elements 4 and 6 by a 5% stiffness reduction. While the test reacts for the damaged elements in the direct test (Figure 2(a)), also the test for the undamaged element 5 reacts strongly, even more than for the actual damaged element 4, causing a false positive. In the minmax test (Figure 2(b)), only the tests for the damaged elements react, as expected.

### 4.3. Sensitivity computation

As seen in the previous section, the sensitivity matrix $\mathcal{J}$ plays an essential role in the statistical tests for the damage localization. It appears in the asymptotic property (9) of the residual, where it results from its Taylor expansion, being its derivative with respect to system parameter $\theta$ and evaluated in $\theta_0$,

$$\mathcal{J} = \lim_{N \to \infty} \frac{\partial}{\partial \theta} \mathbf{E}_{\theta} \left. \frac{1}{\sqrt{N}} \zeta(\theta_0, \mathcal{Y}_N) \right|_{\theta=\theta_0}$$
<figure><img src="image_2.png" /><figcaption>Figure 2: Damage localization tests of shear wall (six elements) with 5% damage in elements 4 and 6.</figcaption></figure>

where $E_θ$ denotes the expectation when the data $Y_N$ is recorded under system under parameter $θ$. Note that the sensitivity matrix is the derivative of the residual where the normalization factor $\sqrt{N}$ is removed, which ensures its (almost sure) convergence. The normalization factor does not disappear of course, and is part of $δ = \sqrt{N}(θ - θ_0)$ in the product $Jδ$ in (9) instead.

This sensitivity matrix is obtained in two steps using the chain rule. First, since all the proposed residuals are linked to the modal parameters of the system, the analytical residual derivative with respect to the modal parameters can be directly obtained based on the definition of the residual as

$$ \mathcal{J}_{\text{modal}} = \lim_{N \to \infty} \frac{\partial}{\partial \eta} \mathbf{E}_{\eta} \left. \frac{1}{\sqrt{N}} \zeta(\theta_0, \mathcal{Y}_N) \right|_{\eta=\eta_0}, \quad (14) $$

where $η$ is a vector containing the modal parameters of the system and $η_0$ is its value in the reference state. Second, the derivative of the modal parameters with respect to the system parameter $θ$ is obtained with the help of an FE model of the structure as

$$ \mathcal{J}_{\text{FE}} = \left. \frac{\partial \eta}{\partial \theta} \right|_{\theta=\theta_0}. $$

Finally, the desired sensitivity matrix is the product of both,

$$ \mathcal{J} = \mathcal{J}_{\text{modal}} \mathcal{J}_{\text{FE}}. $$

In practice, its estimate is obtained by using consistent estimates from measurement data recorded in the reference state and quantities obtained from the FE model, as detailed in the following.

Note that the accuracy of the statistical tests for damage localization depends in particular on the accuracy of this sensitivity matrix up to some scaling constants (see Appendix C), while accuracy of the FE model itself is not an explicit requirement.

### 4.3.1. Sensitivity $\mathcal{J}_{\text{modal}}$ of residual with respect to modal parameters

This part of the sensitivity computation is entirely dependent on the chosen residual and is usually obtained analytically. The modal parameter vector $η$ is user-defined for the problem at hand. Its components can e.g. be the frequencies and real-valued mode shapes for compatibility with a standard FE model under classical damping, or the whole set of modal parameters including frequencies, damping ratios and complex-valued mode shapes. Instead of the natural frequencies, the circular frequencies could be taken, or directly the eigenvalues of the continuous-time system (1). For a convenient computation of the derivatives, $η$ should be defined as a real-valued vector, i.e. instead of containing complex-conjugated pairs of eigenvalues or eigenvectors it should contain their real and imaginary parts, respectively. The only requirements on $η$ are, first, that it should be a complete parameterization, i.e. the modal parameters of all the modes in the frequency range of the measured data of system (3) are included to avoid modal truncation errors, and second, that the derivative of the chosen modal parameter vector $η$ with respect to the desired system parameters will be available in the second part of the sensitivity computation based on the FE model.
To give an example of $\eta$, we assume its composition of the natural frequencies and the pairs of conjugate complex mode shapes, or equivalently their real and imaginary parts of one of each pair, i.e.

$$ \eta = [f_1 \dots f_{m^*} \mathfrak{R}(\varphi_1)^T \dots \mathfrak{R}(\varphi_{m^*})^T \mathfrak{J}(\varphi_1)^T \dots \mathfrak{J}(\varphi_{m^*})^T]^T, \quad (15) $$

where $m^*$ is the number of modes in the frequency range until $1/(2\tau)$. The damping ratios are omitted in this example definition of $\eta$ for simplicity. An example of the sensitivity computation for the subspace residual (5) with respect to $\eta$ as defined in (15) is given in Section 4.3.3.

Besides other estimates related to the residual computation in the reference state, the actual value $\eta_0$ of the modal parameter vector in the reference state is required in the computation of this part of the sensitivity matrix. Since there might be some differences between the modal parameters that are obtained either from the measurement data or from the FE model, the following choices can be made:

**Data-based Computation:** Vector $\eta_0 = \eta_0^{\text{data}}$ is estimated from measurement data in the reference state, as previously proposed in [18, 19]. This requires a precise modal analysis from the reference data. In addition, the identified modes from the data need to be matched with the modes from the FE model, and the mode shapes need to be appropriately scaled to be consistent with the second part $\mathcal{J}_{\text{FE}}$ of the sensitivity computation related to the FE model.

**Model-based Computation:** Vector $\eta_0 = \eta_0^{\text{FE}}$ is obtained from the FE model. We propose this computation as an alternative to the data-based computation. In this way, $\eta_0$ is directly consistent with the second part $\mathcal{J}_{\text{FE}}$ of the sensitivity computation, which is related to the FE model, and no mode matching or scaling are required. Moreover, weakly excited modes that cannot be accurately identified from measurement data can be taken into account in the sensitivity computation, reducing modal truncation.

*Illustrative example.* In Figures 3 and 4, the influence of both computations on the localization performance is illustrated with the shear wall example using the minmax test. First, assuming the same (correct) FE model for data simulation and sensitivity computation, only the first 3 identified modes from the data were used for the computation of $\mathcal{J}_{\text{modal}}$ in Figure 3(a), while 4 modes from the model were used in Figure 3(b). In the former case, only one of two damages can be localized, while in the latter case both damages can be localized correctly. This illustrates that for cases with only few identified modes from the measurements, it might be advantageous to use the model-based sensitivity computation where more modes are available, reducing modal truncation.

In a second setting in Figure 4, both sensitivity computations are compared assuming an incorrect FE model in the sensitivity computation, since in practice FE models hardly describe a structure perfectly. Such a perturbed model of the shear wall was obtained by changing each wall stiffness randomly by up to ±20%. Then, the data-based sensitivity computation indeed improves the situation, which can be seen in Figure 4(a) where 4 identified modes were used, compared to Figure 4(b) where 4 modes from the model were used. In the former case, both damages are localized correctly, while in the latter case only one of two damages can be localized due to the model errors. However, using 6 modes from the model in Figure 4(c), both damages can be localized correctly despite the model errors, illustrating that the proposed model-based computation of $\mathcal{J}_{\text{modal}}$ may indeed remedy modal truncation errors even in the presence of model errors.

<figure><img src="image_3.png" /><figcaption>Figure 3: Damage localization tests of shear wall (six elements) with 5% damage in elements 4 and 6, where sensitivity is computed with the correct FE model.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4: Damage localization tests of shear wall (six elements) with 5% damage in elements 4 and 6, where sensitivity is computed with the perturbed FE model.</figcaption></figure>

### 4.3.2. Sensitivity $J_{FE}$ of modal parameters with respect to system parameter

The derivative of the modal parameters with respect to the system parameter is based on the FE model of the structure. Its computation can be performed analytically or by using mathematical differential approximation methods such as the finite difference approach. In the following these two approaches are described.

#### 4.3.2.1. Analytical computation

The sensitivity of the modal parameters is developed based on the eigenvalue equation of the mechanical model (1)

$$ (M\mu_j^2 + C\mu_j + K)\psi_j = 0, \qquad (16) $$

for each eigenvalue and eigenvector pair $(\mu_j, \psi_j)$, $j = 1, \dots, 2m$ with in general complex-valued $\psi_j \in \mathbb{C}^m$. Note that they are related to the eigenvalues $\lambda_j$ and mode shapes $\varphi_j \in \mathbb{C}^r$ of the monitored system (3), which satisfy

$$ A\phi_j = \lambda_j\phi_j, \quad \varphi_j = C\phi_j, \qquad (17) $$

through

$$ \lambda_j = e^{\mu_j \tau}, \quad \varphi_j = L_j \psi_j \quad \text{where} \quad L_j = L_d + L_v \mu_j + L_a \mu_j^2, \qquad (18) $$

see also (4). Natural frequencies and damping ratios are

$$ f_j = \frac{|\mu_j|}{2\pi}, \quad \xi_j = -\frac{\Re(\mu_j)}{|\mu_j|}, \qquad (19) $$

respectively.

The i-th column of $J_{FE}$ contains the derivative of modal parameter vector $\eta$ with respect to system parameter $\theta^i$ in the parameter vector $\theta = [\theta^1 \ \theta^2 \ \dots \ \theta^l]^T$. It is based on deriving the eigenvalues and eigenvectors $(\mu_j, \psi_j)$ with respect to $\theta^i$ based on (16), supposing that the mass, damping and/or stiffness matrices are functions of $\theta$. Then, $\eta$ is linked to $(\mu_j, \psi_j)$ using the relations (17)–(19). The details of these derivations are presented in Appendix A. Finally,
the modal parameter sensitivity $J_{FE}$ for $\eta$ as defined in (15), as an example, can be composed from these derivatives as

$$ J_{FE} = \begin{bmatrix} \frac{\partial f_1}{\partial \theta^1} & \frac{\partial f_1}{\partial \theta^2} & \dots & \frac{\partial f_1}{\partial \theta^l} \\ \vdots & \vdots & & \vdots \\ \frac{\partial f_{m^*}}{\partial \theta^1} & \frac{\partial f_{m^*}}{\partial \theta^2} & \dots & \frac{\partial f_{m^*}}{\partial \theta^l} \\ \mathbf{R}(\frac{\partial \varphi_1}{\partial \theta^1}) & \mathbf{R}(\frac{\partial \varphi_1}{\partial \theta^2}) & \dots & \mathbf{R}(\frac{\partial \varphi_1}{\partial \theta^l}) \\ \vdots & \vdots & & \vdots \\ \mathbf{R}(\frac{\partial \varphi_{m^*}}{\partial \theta^1}) & \mathbf{R}(\frac{\partial \varphi_{m^*}}{\partial \theta^2}) & \dots & \mathbf{R}(\frac{\partial \varphi_{m^*}}{\partial \theta^l}) \\ \mathfrak{J}(\frac{\partial \varphi_1}{\partial \theta^1}) & \mathfrak{J}(\frac{\partial \varphi_1}{\partial \theta^2}) & \dots & \mathfrak{J}(\frac{\partial \varphi_1}{\partial \theta^l}) \\ \vdots & \vdots & & \vdots \\ \mathfrak{J}(\frac{\partial \varphi_{m^*}}{\partial \theta^1}) & \mathfrak{J}(\frac{\partial \varphi_{m^*}}{\partial \theta^2}) & \dots & \mathfrak{J}(\frac{\partial \varphi_{m^*}}{\partial \theta^l}) \end{bmatrix} \quad (20) $$

### 4.3.2.2. Finite difference method

In evaluating the analytical sensitivities in the previous section, the mass, damping and stiffness matrices of the structure are needed to be evaluated and then their sensitivities to the parameters be calculated. This might be a problem when using commercial software in modeling big structures or complex FE models, where no access to the full mass and stiffness matrices may be available.

Instead, the finite difference method can be used to compute the required sensitivities. To obtain the i-th column of $J_{FE}$, the parameter $\theta_0^i$ is slightly perturbed in the FE model to $\theta_0^i + \Delta\theta_0^i$, while the other components of $\theta_0$ remain the same. Then, the modal parameters $\eta(\theta_0^i + \Delta\theta_0^i)$ are evaluated after the perturbation and

$$ \left. \frac{\partial \eta}{\partial \theta^i} \right|_{\theta=\theta_0} \approx \frac{\eta(\theta_0^i + \Delta\theta_0^i) - \eta(\theta_0)}{\Delta\theta_0^i}. $$

The perturbation $\Delta\theta_0^i$ should be chosen small enough to ensure a small approximation error, but not too small considering machine precision and computational errors.

### 4.3.3. Example: sensitivity of subspace-based residual

The sensitivity of the subspace residual (5) has been developed in [18, 19] and is summarized in this section. Furthermore, important properties regarding its uniqueness and its independence from mode shape scaling are revisited.

From definition (14), the residual sensitivity with respect to modal parameter vector $\eta$ writes as

$$ J_{\text{modal}} = \left. \frac{\partial \text{vec}(S(\theta_0)^T \mathcal{H}(\eta))}{\partial \eta} \right|_{\eta=\eta_0}. $$

Assuming $\eta$ as defined in (15), this sensitivity is computed analytically in three steps with the chain rule. First, the derivative of $\text{vec}(S(\theta_0)^T \mathcal{H})$ is obtained with respect to the set of eigenvalues $\lambda_j$ and mode shapes $\varphi_j$ of the discrete-time system (see (17)-(18)), denoted as $J_{(\lambda, \varphi)}^\zeta$. It is obtained based on the decomposition $\mathcal{H} = \mathcal{O}\mathcal{C}$ into observability and controllability matrix, where $\mathcal{O}$ can be written in the modal basis directly in terms of $\lambda_j$ and $\varphi_j$. Second, the derivative $J_{(\mu, \varphi)}^{(\lambda, \varphi)}$ of the $\lambda_j$'s with respect to the eigenvalues $\mu_j$ of the continuous-time system is obtained, based on (18). Third, the derivative $J_{\eta}^{(\mu, \varphi)}$ with respect to the natural frequencies based on (19) concludes the computation of $J_{\text{modal}} = J_{(\lambda, \varphi)}^\zeta J_{(\mu, \varphi)}^{(\lambda, \varphi)} J_{\eta}^{(\mu, \varphi)}$. The details of this computation are developed in Appendix B.1.

In the final step, the sensitivity $J_{FE}$ of $\eta$ with respect to the system parameter vector $\theta$ is obtained from the FE model as described in Section 4.3.2, leading to the residual sensitivity $J = J_{\text{modal}}J_{FE}$ that is required for the localization approach.

Note that the mode shapes are always identified or computed up to some particular scaling that may be chosen by the user or that may be arbitrary, but which is coherent between $J_{\text{modal}}$ and $J_{FE}$. Furthermore, the mode shape sensitivity with respect to the physical parameters in $J_{FE}$ is not unique (see Appendix A). Nevertheless, the subspace residual sensitivity $J$ is well-defined and unique:
* Mode shape scaling factors are carried through the computation of $J_{\text{modal}}$ and cancel out with the respective scaling factors in $J_{\text{FE}}$, even if the factors are dependent on the model parameters. This is proved in Appendix B.2.

* The non-unique part of the mode shape sensitivity in $J_{\text{FE}}$ lies in the null space of $J_{\text{modal}}$ and therefore cancels out. This is proved in Appendix B.3.

Thanks to these properties, the chosen mode shape scaling and the particular solution of the mode shape sensitivity do not play a role in damage localization with the subspace residual.

### 4.4. Covariance computation

The residual covariance $\Sigma$ is estimated for the chosen residual by a sample covariance computed from data in the reference state, and for some residuals also using data in the tested, possibly damaged, state. The computation may be direct in the reference state, where the residual is evaluated on several datasets (or on several blocks of a long dataset) $Y_{N_b}^k$, $k = 1, \dots, n_b$, each of length $N_b$ [18, 25]. This applies for example to the subspace residual from Section 3.1. Then the sample covariance writes as

$$ \widehat{\Sigma} = \frac{1}{n_b - 1} \sum_{k=1}^{n_b} \zeta(\theta_0, Y_{N_b}^k) \zeta(\theta_0, Y_{N_b}^k)^T. $$

For other residuals, the computation may be indirect, where first a sample covariance related to the data is obtained, e.g. of the system's output covariances, and then propagated in a sensitivity-based approach to the residual covariance [25, 32].

### 4.5. Clustering of parameters

In practice, the number of physical parameters in parameter vector $\theta$ may be large, which is usually based on an FE model. In other cases, a damage in different but close elements of the structure may have a similar effect on the residual. Then, the information provided by the sensors is not sufficient to be able to distinguish if a damage is due to a change in one or another close physical parameter. In this case, the direct test will react for all close parameters, and the minmax test becomes unpredictable since the assumption of $J$ being full column rank is violated in this case. This is illustrated in Figure 5 for the shear wall from Figure 1(b) containing 15 elements but only three sensors. In the example, elements 3 and 5 are damaged. Indeed, besides the damaged elements also many other elements react in the direct test in Figure 5(a), while the minmax test shows a strong reaction only for element 3 in Figure 5(b), missing the damage in element 5.

The closeness of parameters stems from the modal behavior of the structural elements which in turn is related to their geometrical and physical closeness and modal direction in the considered mode shapes. This closeness is reflected in the sensitivity matrix $J$: if changes in some parameters have (nearly) the same effect on the residual, then their respective columns in $J$ must be (nearly) identical. Note that only the direction defined by the columns of $J$ is relevant, but not their magnitude, since the magnitude is cancelled out in the localization tests as shown in Appendix

<figure><img src="image_5.png" /><figcaption>Figure 5: Damage localization tests of shear wall (15 elements) with 10% damage in elements 3 and 5 without using clustering.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 6: Schematic illustration of closeness of sensitivity vectors in three dimensions.</figcaption></figure>

C. Hence, the angle between the columns is the relevant parameter for measuring closeness. Figure 6 illustrates the possible closeness of some parameters in three dimensions.

To be consistent with test statistics of the direct and minmax test, the vectors for clustering are the columns of $\mathcal{J}$ pre-multiplied with $\Sigma^{-1/2}$, which is the square root inverse of $\Sigma$ yielding $\Sigma^{-1} = (\Sigma^{-1/2})^T \Sigma^{-1/2}$. Furthermore, they are normalized to unit norm to remove any effect of their scaling on the clustering approach, since only their angles are relevant. This leads to the normalized vector set $\tilde{\mathcal{J}}_i$, $i = 1, \dots, l$, with

$$ \tilde{\mathcal{J}}_i \stackrel{\text{def}}{=} \frac{1}{\|\Sigma^{-1/2} \mathcal{J}_i\|} \Sigma^{-1/2} \mathcal{J}_i $$

for clustering. The distance between two such vectors is defined as

$$ d_{ij} = 1 - |\tilde{\mathcal{J}}_i^T \tilde{\mathcal{J}}_j|, \qquad (21) $$

which varies between zero when both vectors define the same change direction, and one when they are orthogonal. Closeness of two parameters $\theta^i$ and $\theta^j$ is then reflected in a small value of $d_{ij}$. In the following, two approaches of clustering close parameters are described. The first approach is vector quantization that has been proposed previously for damage localization in [18, 25]. We propose the use of hierarchical clustering in the second approach.

### 4.5.1. Vector quantization: k-means

k-means clustering is a vector quantization approach frequently used in signal processing, image processing and machine learning [30]. In this algorithm, the number of clusters $n_c$ is defined by the user. Then, $n_c$ vectors are selected randomly from all vectors as the centers of the clusters. Subsequently, the other vectors in space are categorized to each of these clusters based on their minimum distance to the centers. Iteratively, the mean of each cluster is calculated and each vector in space is re-associated to the group with the closest center. The iteration converges when no point is re-associated to other clusters.

This algorithm is highly dependent on the number of clusters $n_c$ and the random starting points. It is not guaranteed to converge, and it can converge to local minima. Therefore, different starting points can result in different classifications. Moreover, the number of clusters of structural elements are unknown a priori and the resulting localization tests, in particular the minmax test, are dependent on it. A criterion on closeness of vectors in the resulting clusters is only implicit, and it is possible that some close vectors are categorized in different clusters. Finally, attention needs to be paid on the orientation of the vectors, since a multiplication of a vector with $-1$ does not affect our distance in (21) nor the results of the localization tests, but it affects strongly the outcome of the k-means clustering approach.
### 4.5.2. Hierarchical clustering: complete-linkage

Complete-linkage clustering is an agglomerative hierarchical clustering method [30]. In the first step, each element defines a cluster. In an iterative process, the two clusters with the shortest distance are combined into a larger cluster, where the distance between two clusters is defined as the distance of the cluster elements that are farthest to each other. The result of the clustering can be illustrated in a dendrogram, which shows the sequential cluster fusion and the cluster distance at which each fusion took place. With this classification, a threshold $d_{max}$ needs to be selected, which defines the maximum distance of vectors in each cluster. The clusters in the dendrogram that are combined below this threshold are our final clusters. For each cluster, it is ensured that the pairwise scalar product between the associated vectors has the same sign, otherwise a vector is multiplied by -1. Finally, their mean is computed to obtain the cluster centers.

With this clustering approach, it is not necessary to determine the number of clusters beforehand. Instead, clusters are formed up to a chosen maximum distance between their elements, where our distance measure is defined in (21).

### 4.5.3. Application of clusters in the tests

Define the normalized residual as $\tilde{\zeta} \stackrel{\text{def}}{=} \Sigma^{-1/2}\zeta$, and let the cluster centers be $c_k$, $k = 1, \dots, n_c$ and $k(i)$ the cluster index of the cluster containing element $i$. With this notation, the columns in the sensitivity matrix are replaced by the cluster centers in the tests as follows.

**Direct test:** While the direct test (12) can be applied directly without clustering, it can also be applied to the cluster center $c_{k(i)}$ when testing parameter $i$, which writes

$$ \tilde{t}_{\text{dir}}^i = \frac{\tilde{\zeta}^T c_{k(i)} c_{k(i)}^T \tilde{\zeta}}{c_{k(i)}^T c_{k(i)}}. $$

In this way, the same test value is assigned to all elements in the same cluster.

**Minmax test:** In the computation of the minmax test (13), the sensitivities of the non-tested elements $\tilde{\mathcal{J}}_i$ are replaced by the cluster centers $C_i = [c_1 \dots c_{k(i)-1} \ c_{k(i)+1} \dots c_{n_c}]$ of the clusters that do not contain element $i$. Then, the test for element $i$ writes

$$ \tilde{t}_{\text{mm}}^i = \tilde{\zeta}_i^* \tilde{F}_i^{-1} \tilde{\zeta}_i^*, $$

where the partial residuals are $\tilde{\zeta}_i = \tilde{\mathcal{J}}_i^T \tilde{\zeta}, \tilde{\zeta}_i = C_i^T \tilde{\zeta}, \tilde{F} = [\tilde{\mathcal{J}}_i \ C_i]^T [\tilde{\mathcal{J}}_i \ C_i]$ is the cluster Fisher information with respective partitions analogous to (11), $\tilde{\zeta}_i^* = \tilde{\zeta}_i - \tilde{F}_{i,\tilde{i}} \tilde{F}_{\tilde{i},\tilde{i}}^{-1} \tilde{\zeta}_i$ is the robust residual and $\tilde{F}_i^* = \tilde{F}_{i,i} - \tilde{F}_{i,\tilde{i}} \tilde{F}_{\tilde{i},\tilde{i}}^{-1} \tilde{F}_{\tilde{i},i}$. In this way, test values can be computed for all elements in a cluster, while being insensitive to changes in the other clusters. Keeping in mind that all elements in a cluster will yield similar test values, this still can give a slightly higher localization resolution compared to computing test values only for the cluster centers. In the latter case, $\tilde{\mathcal{J}}_i$ is replaced by $c_{k(i)}$ in the computation of $\tilde{t}_{\text{mm}}^i$.

Note that due to the invariance of the tests to the scaling of the sensitivity matrix (see Appendix C), the different scaling between the cluster centers and the columns of the sensitivity matrix does not affect the tests. Hence, the replacement of the columns of $\Sigma^{-1/2}\tilde{\mathcal{J}}$ by the cluster centers is coherent in the tests above.

<figure><img src="image_7.png" /><figcaption>Figure 7: Clustering for shear wall (15 elements).</figcaption></figure>
*Illustrative example.* Using the same data as in Figure 5 for localization on the 15 element shear wall, the elements were clustered with the *k*-means and the hierarchical clustering approach as described above. The resulting clusters are depicted in Figure 7. For the comparison of both approaches, a threshold of 0.3 was chosen in the dendrogram of the hierarchical approach in Figure 7(a), which leads to 7 clusters, and 7 clusters were chosen in the *k*-means approach, see Figure 7(b). The resulting minmax localization tests are shown in Figure 8, where the hierarchical approach leads to a clear localization of the damage in elements 3 and 5, while the undamaged elements react more strongly when using the *k*-means approach for clustering.

<figure><img src="image_8.png" /><figcaption>Figure 8: Damage localization tests of shear wall (15 elements) with 10% damage in elements 3 and 5 with clustering.</figcaption></figure>

**5. Case study: Yellow Frame**

The Yellow frame is a modular 4 story, scaled (1/3) steel frame reestablished during this research at the University of British Columbia (UBC), shown in Figure 9. Previously, this structure has been used as an IASC/ASCE benchmark structure since 2002 [31]. Several damage scenarios are designed and tested by removal of braces of the structure.

The structure is 3.6 m high and is composed of 2 spans in each direction with a total length of 2.5 m. Each floor of the structure is carrying dead loads applied to the structure by using 4 steel plates distributed on each level. For

<figure><img src="image_9.png" /><figcaption>Figure 9: Photo of the Yellow Frame structure (south-east corner).</figcaption></figure>
the lateral stability, four pairs of threaded steel rods with a diameter of 12 mm are used as braces on each side of the structure in each floor. In each floor, three sensors are located at the north, south and west side of the structure, amounting to altogether 12 sensors. In Figure 10, a plan of each level with the mass plates and sensor locations is shown, as well as the numbering of the braces. In our study, the following damage scenarios are considered in detail:

* Scenario 1: removal of braces 21 and 23

* Scenario 2: removal of braces 2, 4, 18 and 20

The goal is to identify the removed braces. Hence, the stiffnesses of each of the 32 braces are chosen as the physical parameter vector $\theta$. In order to compute the sensitivity $J_{FE}$ of the modal parameters with respect to the physical

<figure><img src="image_10.png" /><figcaption>Figure 10: Left: plan of each level with sensor locations, right: brace numbering.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 11: Finite element model of the Yellow frame.</figcaption></figure>
parameters, an FE model of the Yellow frame has been created using the Abaqus software (see Figure 11). In the model, the section properties of the elements of the Yellow frame are used in modeling the beams, braces and columns. The plates on each floor are modeled as four lumped masses in the four corners of each plate on the structure. Since the plates are bolted and connected with pretensioned rods to their surrounding beams (in friction), their contribution to the stiffness of the structure for ambient vibration is not negligible. Therefore the plates are also modeled as two parallel beams. The base of the structure is modeled as a fixed connection to the ground. The connections of beams and columns are also modeled as fixed connections and the braces are connected as moment free hinge connections to the structure. Because each group of two braces in each floor at each span is only under axial force, they are modeled as one element with cross section area equal to the total area of both braces. Localizing damage in one brace element indicates the possibility of damage in both of these parallel braces. The FE model is not updated based on the ambient vibration measurements, and is used as is for the sensitivity computation without additional calibration. The first four mode shapes computed from the model are depicted in Figure 12. Note that the model needs to be approximately representative of the dynamic behavior of the structure, although it is not required to match perfectly. It is only used for the sensitivity computation in the reference state, and it is not needed for updating or testing. The computation of $J_{FE}$ has been made using the finite difference method for the first ten modes of the FE model.

From ambient vibration measurements, the first six identified modes correspond to the modes of the FE model, confirming the validity of the model. Two further modes correspond to the ninth and tenth mode of the model, so altogether eight identified modes could be matched. This is in particular relevant when the *data-based* computation of the sensitivity $J_{modal}$ of the residual with respect to the modal parameters is used. The first four identified mode shapes are shown in Figure 13. Note that due to the symmetry of the structure some modes appear in pairs with frequencies that are distinct but close to each other. While the respective mode shapes of these modes are distinct in x and y directions in the FE model, the identified mode shapes from the data are mixed and need to be decoupled and scaled beforehand. Besides using the data-based computation of $J_{modal}$, the results will be compared to the proposed *model-based* computation of $J_{modal}$ in the following, where all ten modes from the model will be used and no mode matching is necessary. Clustering is done in both cases with the proposed hierarchical linkage clustering approach for the minmax tests in Figures 16 and 17, and results are compared to *k*-means clustering in Figure 19. In Figures 14 and 15, the dendrograms are shown for linkage clustering with the data-based and with the model-based sensitivity

<figure><img src="image_12.png" /><figcaption>Figure 12: First four mode shapes from the FE model.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: First four mode shapes identified from the data.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: Dendrogram depicting the hierarchical complete-linkage clustering with the data-based sensitivity computation.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: Dendrogram depicting the hierarchical complete-linkage clustering with the model-based sensitivity computation.</figcaption></figure>

computation, respectively, where a threshold of 0.15 was chosen for the closeness of elements. It can be seen that neighboring braces are always in the same cluster, and furthermore there are few clusters containing four braces.

In Figure 16, the direct and minmax test results are shown for the data-based and the model-based sensitivity computation for damage scenario 1, where braces 21 and 23 have been removed. Using the data-based sensitivity computation, the damage cannot be localized with the direct test in Figure 16(a), while the minmax test in Figure 16(b) shows the highest test values at the damaged braces. However, results are much clearer with the model-based sensitivity computation in Figures 16(c) and 16(d), where the direct test and minmax test show a similar performance and both tests localize the damage correctly. Since both braces 21 and 23 are in the same cluster, note that this damage case corresponds actually to *one* damaged element.

Damage scenario 2, where braces 2, 4, 18 and 20 are removed, corresponds to a truly multiple damage scenario since those elements are not in the same cluster. Using the data-based sensitivity computation, the damage cannot be localized neither with the direct test nor with the minmax test in Figures 17(a) and 17(b). Using the model-based sensitivity computation, the direct test in Figure 17(c) reacts well for the damaged braces, but also reacts for the undamaged braces 6, 8, 30 and 32 which makes a correct damage localization impossible. Finally, the minmax test in Figure 17(d) localizes the damage correctly.

In the minmax tests in Figures 16 and 17, the proposed hierarchical complete-linkage clustering was applied. To
<figure><img src="image_16.png" /><figcaption>Figure 16: Damage localization tests for scenario 1: removal of braces 21 and 23.</figcaption></figure>

compare to previously used *k*-means clustering, two datasets from damage scenarios 1 and 2 were considered. Here, only the model-based sensitivity computation was used which performed best in Figures 16 and 17. The clusters obtained from one *k*-means clustering instance are shown in Figure 18, where the number of clusters was set to 15, analogous to the hierarchical clustering approach in Figure 15. Note that due to the random initialization, the clusters are different when *k*-means is run another time. The localization results are shown in Figure 19. The localization of damage scenario 1 in Figure 19(a) is not successful, which may be caused by close elements being in different clusters, like neighboring damaged braces 21 and 23. Regarding damage scenario 2 in Figure 19(b), only the damage in braces 18 and 20 can be localized, while the tests for other braces (in particular 17 and 19) also react, and the damage in braces 2 and 4 cannot be localized.

From the results in Figures 16, 17 and 19 it can be seen that the minmax approach using the model-based sensitivity computation and the hierarchical complete-linkage clustering approach is the most robust approach for damage localization, illustrating the new developments of this paper.
<figure><img src="image_17.png" /><figcaption>Figure 17: Damage localization tests for scenario 2: removal of braces 2, 4, 18 and 20.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 18: Clusters obtained from k-means with the model-based sensitivity computation.</figcaption></figure>

Finally, this approach is applied to two further challenging damage scenarios. In Scenario 3, one brace of each pair is removed on all four sides of the frame in the ground floor, namely braces 1–4 and 17–20. As shown in Figure 20(a), the damage on three of the four sides is localized correctly, while the damage on the east side in braces 1 and 3 is
<figure><img src="image_19.png" /><figcaption>Figure 19: Minmax tests with *k*-means clustering, using model-based computation of *J*<sub>modal</sub>.</figcaption></figure>

not localized. This may be due to the absence of sensors on this side of the frame, but a proper analysis of damage detectability based on the Fisher information may give more insight in future work. In Scenario 4, braces 6, 8, 29 and 31 are removed, which includes a brace pair in the highest floor. The tests for all the damaged elements react, as shown in Figure 20(b), however the tests for some of the neighboring elements show also some reaction in this case, namely for elements 9, 30 and 32, which may be due to more difficult detectability in the highest floor.

It should be noted that these results are the first successful application of the considered generic damage localization framework to experimental data, which has become possible with the developments of this paper.

<figure><img src="image_20.png" /><figcaption>Figure 20: Robust damage localization (minmax tests with model-based computation of $J_{modal}$ and linkage clustering) on two further damage scenarios.</figcaption></figure>
## 6. Conclusions

In this paper, the theoretical framework for statistical sensitivity-based damage localization was developed into a working damage localization method that is robust for real applications. The essential steps are the use of minmax tests for detecting which structural elements are responsible for the change in the residual, the model-based computation of the residual sensitivity and the use of the hierarchical complete-linkage clustering approach. The resulting method offers a flexible and generic framework for damage localization, taking into account both physical model information and measurement based uncertainties. The method was successfully applied to experimental data from different damage cases of the Yellow Frame, a 3D steel frame, at the University of British Columbia.

## Appendix A. Analytical sensitivity of modal parameters with respect to physical parameters

Deriving (16) leads to the expression [19, 43]

$$ \frac{\partial \mu_j}{\partial \theta^i} (2\mu_j M + C) \psi_j + \left( \mu_j^2 \frac{\partial M}{\partial \theta^i} + \mu_j \frac{\partial C}{\partial \theta^i} + \frac{\partial K}{\partial \theta^i} \right) \psi_j + (M\mu_j^2 + C\mu_j + K) \frac{\partial \psi_j}{\partial \theta^i} = 0, $$

from where the equations for the eigenvalue and eigenvector sensitivities follow as

$$ \frac{\partial \mu_j}{\partial \theta^i} = - \frac{\psi_j^T \left( \mu_j^2 \frac{\partial M}{\partial \theta^i} + \mu_j \frac{\partial C}{\partial \theta^i} + \frac{\partial K}{\partial \theta^i} \right) \psi_j}{\psi_j^T (2\mu_j M + C) \psi_j}, \quad (A.1) $$

$$ (M\mu_j^2 + C\mu_j + K)\frac{\partial \psi_j}{\partial \theta^i} = -\frac{\partial \mu_j}{\partial \theta^i} (2\mu_j M + C)\psi_j - \left(\mu_j^2 \frac{\partial M}{\partial \theta^i} + \mu_j \frac{\partial C}{\partial \theta^i} + \frac{\partial K}{\partial \theta^i}\right)\psi_j. \quad (A.2) $$

Since $(M\mu_j^2 + C\mu_j + K)\psi_j = 0$, the solution for the eigenvector sensitivity is not unique. Thus, for any particular solution $\frac{\partial \psi_j}{\partial \theta^i}$ of (A.2), $\frac{\partial \psi_j}{\partial \theta^i} + \alpha \psi_j$ is also a solution with any $\alpha \in \mathbb{C}$. A particular solution is, e.g.,

$$ \frac{\partial \psi_j}{\partial \theta^i} = (M\mu_j^2 + C\mu_j + K)^{\dagger} \left( -\frac{\partial \mu_j}{\partial \theta^i} (2\mu_j M + C)\psi_j - \left( \mu_j^2 \frac{\partial M}{\partial \theta^i} + \mu_j \frac{\partial C}{\partial \theta^i} + \frac{\partial K}{\partial \theta^i} \right) \psi_j \right), \quad (A.3) $$

where $^\dagger$ denotes the pseudoinverse. In general, any solution can be chosen. In Appendix B.3 it is shown that the final sensitivity $\mathcal{J}$ of the subspace residual is unique, even though the above eigenvector sensitivity is not unique.

In Equations (A.1) and (A.3), the sensitivities of the eigenvalues and eigenvectors $(\mu_j, \psi_j)$ with respect to the system parameters are obtained. Since the sensitivity of $\eta$ with respect to $\theta$ is required, the relationship between $\eta$ and $(\mu_j, \psi_j)$ needs to be derived for all modes $j$ that are also present in the data. When $\eta$ is defined as in (15), the derivatives of the natural frequencies and mode shapes at the sensor coordinates are required. Then, based on (19) and (A.1) it follows

$$ \frac{\partial f_j}{\partial \theta^i} = \frac{\partial f_j}{\partial \mathcal{R}(\mu_j)} \mathcal{R} \left( \frac{\partial \mu_j}{\partial \theta^i} \right) + \frac{\partial f_j}{\partial \mathcal{I}(\mu_j)} \mathcal{I} \left( \frac{\partial \mu_j}{\partial \theta^i} \right) \\ = \frac{1}{2\pi|\mu_j|} \left( \mathcal{R}(\mu_j) \mathcal{R} \left( \frac{\partial \mu_j}{\partial \theta^i} \right) + \mathcal{I}(\mu_j) \mathcal{I} \left( \frac{\partial \mu_j}{\partial \theta^i} \right) \right) = \frac{1}{2\pi|\mu_j|} \mathcal{R} \left( \frac{\partial \mu_j}{\partial \theta^i} \right), $$

and from (18) and (A.3) it follows

$$ \frac{\partial \varphi_j}{\partial \theta^i} = \frac{\partial L_j}{\partial \theta^i} \psi_j + L_j \frac{\partial \psi_j}{\partial \theta^i}, \quad \text{where} \quad \frac{\partial L_j}{\partial \theta^i} = (L_v + 2\mu_j L_a) \frac{\partial \mu_j}{\partial \theta^i}. \quad (A.4) $$

Note that this general expression is based on relation (18) between the mode shapes $\varphi_j$ contained in $\eta$ and the eigenvectors $\psi_j$ from the FE model. It is in particular applicable when different sensor types are used at the same time, e.g. measuring displacements and accelerations. In practice, we may have the following simplifications regarding the mode shape scaling when only one sensor type is used, e.g. accelerometers:
* If the data-based computation is used for $\mathcal{J}_{\text{modal}}$, where the mode shapes $\varphi_j$ are obtained from the data and need to be rescaled to match the corresponding mode shapes $L_a\psi_j$ from the FE model through some scaling factor $s_j \in \mathbb{C}$ with $\varphi_j \approx s_j L_a \psi_j$, e.g. through $s_j = \varphi_j^H \varphi_j / (\varphi_j^H L_a \psi_j)$, then

$$ \frac{\partial \varphi_j}{\partial \theta^i} = \frac{\partial s_j}{\partial \theta^i} L_a \psi_j + s_j L_a \frac{\partial \psi_j}{\partial \theta^i}. $$

Note that the first term $\frac{\partial s_j}{\partial \theta^i}$ cancels out in the sensitivity of the subspace residual as shown in Appendix B.2, effectively leading to $\frac{\partial \varphi_j}{\partial \theta^i} = s_j L_a \frac{\partial \psi_j}{\partial \theta^i}$ in this case.

* If the model-based computation is used for $\mathcal{J}_{\text{modal}}$, where the mode shapes $\varphi_j$ are directly obtained from the FE model as $\varphi_j = L_a \psi_j$ without rescaling, then the components of the FE eigenvector sensitivity are simply selected at the measured DOFs as

$$ \frac{\partial \varphi_j}{\partial \theta^i} = L_a \frac{\partial \psi_j}{\partial \theta^i}. $$

# Appendix B. Subspace residual sensitivity analysis

## Appendix B.1. Residual sensitivity computation

For the first step of the derivative computation with respect to the eigenvalues $\lambda_j$ and mode shapes $\varphi_j$ of the discrete-time system (see (17)-(18)), the intermediate parameter vector $\tilde{\eta}$ is defined as follows. Since all modes appear in conjugated complex pairs in structural vibration analysis, the eigenvalues and mode shapes of system (3) can be separated into conjugated complex pairs, being $\lambda_j$ and $\varphi_j$ where $\mathfrak{J}(\lambda_j) > 0$, $j = 1, \dots, m^*$, and their conjugated complex counterparts being $\bar{\lambda}_j$ and $\bar{\varphi}_j$. Define $\Delta_c \stackrel{\text{def}}{=} \text{diag}(\lambda_1, \dots, \lambda_{m^*})$, $\Phi_c \stackrel{\text{def}}{=} [\varphi_1 \dots \varphi_{m^*}]$ and $\tilde{\eta}_c \stackrel{\text{def}}{=} [\lambda_1 \dots \lambda_{m^*} \varphi_1^T \dots \varphi_{m^*}^T]^T$, containing each one value of each conjugated complex pair. Then, the system is fully described by the complex-valued parameter vector $[\tilde{\eta}_c^T \tilde{\eta}_c^T]^T$, and we define the corresponding real-valued parameter vector $\tilde{\eta} \stackrel{\text{def}}{=} [\mathfrak{R}(\tilde{\eta}_c)^T \mathfrak{J}(\tilde{\eta}_c)^T]^T$. The goal of the first step in the derivative computation is to obtain

$$ \mathcal{J}_{(\lambda, \varphi)}^{\zeta} \stackrel{\text{def}}{=} \left. \frac{\partial \text{vec}(S(\theta_0)^T \mathcal{H}(\tilde{\eta}))}{\partial \tilde{\eta}} \right|_{\eta=\eta_0} \qquad (B.1) $$

This derivative is obtained through the decomposition $\mathcal{H} = \mathcal{O}\mathcal{C}$, where $\mathcal{O} = \mathcal{O}(\tilde{\eta})$ is defined in the modal basis in terms of the eigenvalues $\lambda_j$ and mode shapes $\varphi_j$ of the discrete-time system as follows. In the complex format, it is obtained from a similarity transform of the state space system, where $A = \text{diag}(\Delta_c, \Delta_c)$ and $C = [\Phi_c \quad \Phi_c]$. In the following, we use the equivalent real-valued format for ease of computation which is obtained after another simple similarity transform, yielding

$$ \mathcal{O}(\tilde{\eta}) \stackrel{\text{def}}{=} [\mathfrak{R}(\mathcal{O}_c(\tilde{\eta}_c)) \quad \mathfrak{J}(\mathcal{O}_c(\tilde{\eta}_c))] = \begin{bmatrix} \mathfrak{R}(\Phi_c) & \mathfrak{J}(\Phi_c) \\ \mathfrak{R}(\Phi_c \Delta_c) & \mathfrak{J}(\Phi_c \Delta_c) \\ \vdots & \vdots \\ \mathfrak{R}(\Phi_c \Delta_c^p) & \mathfrak{J}(\Phi_c \Delta_c^p) \end{bmatrix}, \quad \text{where } \mathcal{O}_c(\tilde{\eta}_c) \stackrel{\text{def}}{=} \begin{bmatrix} \Phi_c \\ \Phi_c \Delta_c \\ \vdots \\ \Phi_c \Delta_c^p \end{bmatrix}. \qquad (B.2) $$

Then, (B.1) is obtained in this first step as follows. It holds for each entry $i$ of $\tilde{\eta}$

$$ \left. \frac{\partial(S(\theta_0)^T \mathcal{O}(\tilde{\eta}) \mathcal{C}(\tilde{\eta}))}{\partial \tilde{\eta}_i} \right|_{\eta=\eta_0} = S(\theta_0)^T \left. \frac{\partial \mathcal{O}(\tilde{\eta})}{\partial \tilde{\eta}_i} \right|_{\eta=\eta_0} C(\eta_0) + S(\theta_0)^T \mathcal{O}(\eta_0) \left. \frac{\partial C(\tilde{\eta})}{\partial \tilde{\eta}_i} \right|_{\eta=\eta_0}. \qquad (B.3) $$

Since $S(\theta_0)^T \mathcal{O}(\eta_0) = 0$, the second term vanishes, noting that the physical and modal parameter vectors $\theta_0$ and $\eta_0$ describe the same system state, respectively. Vectorizing the expression, taking the derivative with respect to the entire vector $\tilde{\eta}$ and replacing $C(\eta_0) = \mathcal{O}(\eta_0)^{\dagger} \mathcal{H}(\eta_0)$ leads to

$$ \mathcal{J}_{(\lambda, \varphi)}^{\zeta} = (\mathcal{O}(\eta_0)^{\dagger} \mathcal{H}(\eta_0) \otimes S(\theta_0))^{\mathrm{T}} \mathcal{O}'(\eta_0), \quad \text{where } \mathcal{O}'(\eta_0) \stackrel{\text{def}}{=} \left. \frac{\partial \text{vec}(\mathcal{O}(\tilde{\eta}))}{\partial \tilde{\eta}} \right|_{\eta=\eta_0}, \qquad (B.4) $$
where ⊗ denotes the Kronecker product. The derivative $O'(\eta_0)$ is calculated analytically based on the complex-valued derivative of $O_c(\tilde{\eta}_c)$ [18, 22, 25]. Define

$$ \Lambda_j \stackrel{\text{def}}{=} [1 \ \lambda_j \ \lambda_j^2 \ \dots \ \lambda_j^p]^T, \quad \Lambda'_j \stackrel{\text{def}}{=} [0 \ 1 \ 2\lambda_j \ \dots \ p\lambda_j^{p-1}]^T, $$

then

$$ O'_c(\eta_0) \stackrel{\text{def}}{=} \left. \frac{\partial \operatorname{vec}(O_c(\tilde{\eta}_c))}{\partial \tilde{\eta}_c} \right|_{\eta=\eta_0} = \left[ \begin{array}{cccc} \Lambda'_1 \otimes \varphi_1 & 0 & \Lambda_1 \otimes I_r & 0 \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \Lambda'_{m^*} \otimes \varphi_{m^*} & 0 & \Lambda_{m^*} \otimes I_r \end{array} \right]_{\eta=\eta_0} \quad (B.5) $$

and the derivative of the real-valued parametric observability matrix $O'(\eta_0)$ in (B.4) yields

$$ O'(\eta_0) = \begin{bmatrix} \mathfrak{R}(O'_c(\eta_0)) & -\mathfrak{I}(O'_c(\eta_0)) \\ \mathfrak{I}(O'_c(\eta_0)) & \mathfrak{R}(O'_c(\eta_0)) \end{bmatrix}. \quad (B.6) $$

In the second step, the derivatives of $\lambda_j$ with respect to $\mu_j$ and $f_j$ are required. From (18) it follows

$$ \begin{bmatrix} \frac{\partial \mathfrak{R}(\lambda_j)}{\partial \mathfrak{R}(\mu_j)} & \frac{\partial \mathfrak{R}(\lambda_j)}{\partial \mathfrak{I}(\mu_j)} \\ \frac{\partial \mathfrak{I}(\lambda_j)}{\partial \mathfrak{R}(\mu_j)} & \frac{\partial \mathfrak{I}(\lambda_j)}{\partial \mathfrak{I}(\mu_j)} \end{bmatrix} = \tau \begin{bmatrix} \mathfrak{R}(\lambda_j) & -\mathfrak{I}(\lambda_j) \\ \mathfrak{I}(\lambda_j) & \mathfrak{R}(\lambda_j) \end{bmatrix}, $$

and thus

$$ J_{(\mu, \varphi)}^{(\lambda, \varphi)} = \begin{bmatrix} \begin{array}{c} \tau \mathfrak{R}(\lambda_1) \\ \vdots \\ \tau \mathfrak{R}(\lambda_{m^*}) \end{array} & \begin{array}{c} -\mathfrak{I}(\lambda_1) \\ \vdots \\ -\mathfrak{I}(\lambda_{m^*}) \end{array} \\ \begin{array}{c} \tau \mathfrak{I}(\lambda_1) \\ \vdots \\ \tau \mathfrak{I}(\lambda_{m^*}) \end{array} & \begin{array}{c} I_{m^*r} \\ \mathfrak{R}(\lambda_1) \\ \vdots \\ \mathfrak{R}(\lambda_{m^*}) \end{array} \\ \begin{array}{c} 0_{m^*r} \end{array} & \begin{array}{c} 0_{m^*r} \\ I_{m^*r} \end{array} \end{bmatrix}_{\eta=\eta_0} \quad (B.7) $$

where $I_a$ and $0_a$ are the identity and zero matrix of size $a \times a$, respectively. Following from the relation $\mu_j = -2\pi f_j \xi_j + 2\pi f_j \sqrt{1 - \xi_j^2}$, the derivatives with respect to the frequencies yield

$$ \frac{\partial \mathfrak{R}(\mu_j)}{\partial f_j} = -2\pi \xi_j, \quad \frac{\partial \mathfrak{I}(\mu_j)}{\partial f_j} = 2\pi \sqrt{1 - \xi_j^2}, $$

and thus

$$ J_{\eta}^{(\mu, \varphi)} = \begin{bmatrix} -2\pi\xi_1 & & & & \\ & \ddots & & & \\ & & -2\pi\xi_{m^*} & & \\ & & & I_{m^*r} & 0_{m^*r} \\ & & & & \\ 2\pi\sqrt{1-\xi_1^2} & & & & \\ & \ddots & & & \\ & & 2\pi\sqrt{1-\xi_{m^*}^2} & & \\ & & & 0_{m^*r} & I_{m^*r} \end{bmatrix}_{\eta=\eta_0} \quad (B.8) $$

Finally, an estimate of $J_{\text{modal}}$ is obtained by plugging into (B.4) the estimate of $S(\theta_0)$ that is used for the residual computation, and an estimate $\tilde{H}$ of $H(\eta_0)$ from the reference state. The matrices $O(\eta_0)$, $O'(\eta_0)$, $J_{(\mu, \varphi)}^{(\lambda, \varphi)}$ and $J_{\eta}^{(\mu, \varphi)}$ are computed using the modal parameters and the respective eigenvalues of the discrete-time and continuous-time systems, either based on $\eta_0 = \eta_0^{\text{data}}$ obtained from modal analysis of the reference data when using the data-based computation, or based on $\eta_0 = \eta_0^{\text{FE}}$ obtained from the FE model when using the model-based computation.
Appendix B.2. Independence from mode shape scaling

First it is shown how the mode shape scaling affects the sensitivity matrix. Suppose the set of mode shapes $\varphi_j$, $j = 1, \dots, m^*$, be given in the parameter vector $\eta$, and let the differently scaled mode shapes be $\varphi_j^d = d_j \varphi_j$ with deterministic scaling factors $d_j \in \mathbb{C}$ in respective parameter vector $\eta^d$, and $D \stackrel{\text{def}}{=} \text{diag}(d_1, \dots, d_{m^*})$. Comparing to (B.4), the derivative containing the $\varphi_j^d$'s is evaluated in the following. It holds $\Phi_c^d = \Phi_c D$, thus $O_c(\tilde{\eta}_c^d) = O_c(\tilde{\eta}_c)D$ and

$$O(\tilde{\eta}^d) = O(\tilde{\eta}) D_{\text{Re}}^T, \quad D_{\text{Re}} \stackrel{\text{def}}{=} \begin{bmatrix} \mathfrak{R}(D) & -\mathfrak{J}(D) \\ \mathfrak{J}(D) & \mathfrak{R}(D) \end{bmatrix} \qquad (B.9)$$

in (B.2). Then, since $D$ is invertible, it holds for the Kronecker product in (B.4)

$$\begin{align*}
(O(\eta_0^d)^\dagger H(\eta_0^d) \otimes S(\theta_0))^T &= H(\eta_0)^T (O(\eta_0^d)^\dagger)^T \otimes S(\theta_0)^T \\
&= H(\eta_0)^T (D_{\text{Re}}^{-T} O(\eta_0)^\dagger)^T \otimes S(\theta_0)^T \\
&= (O(\eta_0)^\dagger H(\eta_0) \otimes S(\theta_0))^T (D_{\text{Re}}^{-1} \otimes I_{(p+1)r}) \tag{B.10}
\end{align*}$$

In (B.5) it holds

$$O'_{c}(\eta_{0}^{d}) = \begin{bmatrix} \Lambda'_{1} \otimes d_{1}\varphi_{1} & 0 & d_{1}(\Lambda_{1} \otimes I_{r})\frac{1}{d_{1}} & 0 \\ 0 & \cdots & 0 & d_{m^{*}}(\Lambda_{m^{*}} \otimes I_{r})\frac{1}{d_{m^{*}}} \\ 0 & \Lambda'_{m^{*}} \otimes d_{m^{*}}\varphi_{m^{*}} & 0 & 0 \\ 0 & \cdots & 0 & \cdots \\ 0 & \cdots & D^{-1} \otimes I_{r} & 0 \end{bmatrix} \\ = (D \otimes I_{(p+1)r}) O'_{c}(\eta_{0}) E, \quad E \stackrel{\text{def}}{=} \begin{bmatrix} I_{m^{*}} & 0 \\ 0 & D^{-1} \otimes I_{r} \end{bmatrix}, \quad E_{\text{Re}} \stackrel{\text{def}}{=} \begin{bmatrix} \mathfrak{R}(E) & -\mathfrak{J}(E) \\ \mathfrak{J}(E) & \mathfrak{R}(E) \end{bmatrix},$$

and, following (B.6),

$$O'(\eta_0^d) = (D_{\text{Re}} \otimes I_{(p+1)r}) O'(\eta_0) E_{\text{Re}}. \qquad (B.11)$$

Finally, since the mode shape part is not affected when calculating the derivatives with respect to the frequencies, it follows from (B.4), (B.7), (B.8), (B.10) and (B.11)

$$\mathcal{J}_{\text{modal}}^{d} = \mathcal{J}_{\text{modal}} F, \quad F \stackrel{\text{def}}{=} \begin{bmatrix} I_{m^{*}} & D_{\text{Re}}^{-1} \otimes I_{r} \end{bmatrix}.$$

Since $\varphi_j^d = d_j \varphi_j$, the derivative of $\varphi_j^d$ with respect to the physical parameters is premultiplied by $d_j$ compared to the derivative of $\varphi_j$, leading to

$$\mathcal{J}_{\text{FE}}^{d} = F^{-1} \mathcal{J}_{\text{FE}}$$

and finally $\mathcal{J}^d = \mathcal{J}_{\text{modal}}^d \mathcal{J}_{\text{FE}}^d = \mathcal{J}_{\text{modal}} \mathcal{J}_{\text{FE}} = \mathcal{J}$, independently of the mode shape scaling.

Suppose now that the scaling factors depend on the parameter vector, e.g. when derived from the FE model. Then, a change in the scaling factors cannot simply be replaced in the derivative $\mathcal{J}$ as it was done through relation (B.9), but the derivative of this relation needs to be taken into account. Similar to (B.3), the residual derivative with respect to one parameter component writes as

$$\left. \frac{\partial(S(\theta_0)^T O(\tilde{\eta}^d) D_{\text{Re}}(\tilde{\eta}^d)^T C(\tilde{\eta}^d))}{\partial \tilde{\eta}_i^d} \right|_{\eta=\eta_0} = S(\theta_0)^T \left. \frac{\partial O(\tilde{\eta}^d)}{\partial \tilde{\eta}_i^d} \right|_{\eta=\eta_0} D_{\text{Re}}(\eta_0^d)^T C(\eta_0^d) + S(\theta_0)^T O(\eta_0^d) \left. \frac{\partial(D_{\text{Re}}(\tilde{\eta}^d)^T C(\tilde{\eta}^d))}{\partial \tilde{\eta}_i^d} \right|_{\eta=\eta_0},$$

and since the second term vanishes thanks to $S(\theta_0)^T O(\eta_0^d) = 0$, also the derivative of possible parameter dependent scaling factors vanishes. Hence, $\mathcal{J}$ is indeed independent of the chosen mode shape scaling, whether it is dependent on the physical parameters or not.
Appendix B.3. Independence from non-uniqueness of mode shape derivative

In Appendix A it was shown that the eigenvector sensitivity $\frac{\partial \psi_j}{\partial \theta^i}$ is not unique and any $\frac{\partial \psi_j}{\partial \theta^i} + \alpha_{ij} \psi_j$ with $\alpha_{ij} \in \mathbb{C}$ is valid in the computation of $\mathcal{J}_{FE}$ for all modes $j = 1, \dots, m^*$ and parameters $i = 1, \dots, l$. Nevertheless, it is shown in the following that the residual sensitivity $\mathcal{J} = \mathcal{J}_{modal}\mathcal{J}_{FE}$ is independent from this non-uniqueness, by showing that the contribution of $\alpha_{ij} \psi_j$ in $\mathcal{J}_{FE}$ lies in the null space of $\mathcal{J}_{modal}$ for any $i$ and $j$. The mode shape sensitivity follows from the eigenvector sensitivity $\frac{\partial \psi_j}{\partial \theta^i} + \alpha_j \psi_j$ from (A.4) and (18) as

$$ \frac{\partial L_j}{\partial \theta^i} \psi_j + L_j \left( \frac{\partial \psi_j}{\partial \theta^i} + \alpha_{ij} \psi_j \right) = \frac{\partial \varphi_j}{\partial \theta^i} + \alpha_{ij} \varphi_j, \qquad (B.12) $$

and denote the sensitivity matrices $\mathcal{J}_{FE}$ and $\mathcal{J}_{FE}^{\alpha}$ when replacing $\frac{\partial \varphi_j}{\partial \theta^i}$ by $\frac{\partial \varphi_j}{\partial \theta^i} + \alpha_{ij} \varphi_j$ in (20), respectively. Analyze the difference $d = \mathcal{J}_{modal}\mathcal{J}_{FE} - \mathcal{J}_{modal}\mathcal{J}_{FE}^{\alpha}$. In (B.4) it holds

$$ \begin{aligned} \mathcal{J}_{(\lambda, \varphi)}^{\zeta} &= (\mathcal{O}(\eta_0)^{\dagger} \mathcal{H}(\eta_0) \otimes I_s)^T (I_n \otimes S(\theta_0)^T) \mathcal{O}'(\eta_0) \\ &= (\mathcal{O}(\eta_0)^{\dagger} \mathcal{H}(\eta_0) \otimes I_s)^T \left[ \begin{array}{cc} (I_{m^*} \otimes S(\theta_0)^T) \mathfrak{R}(\mathcal{O}'_c(\eta_0)) & -(I_{m^*} \otimes S(\theta_0)^T) \mathfrak{F}(\mathcal{O}'_c(\eta_0)) \\ (I_{m^*} \otimes S(\theta_0)^T) \mathfrak{F}(\mathcal{O}'_c(\eta_0)) & (I_{m^*} \otimes S(\theta_0)^T) \mathfrak{R}(\mathcal{O}'_c(\eta_0)) \end{array} \right], \end{aligned} \qquad (B.13) $$

where $s = (p+1)r-n$. Since the mode shape part in the derivative $\mathcal{J}_{modal} = \mathcal{J}_{(\lambda, \varphi)}^{\zeta} \mathcal{J}_{(\mu, \varphi)}^{(\lambda, \varphi)} \mathcal{J}_{\eta}^{(\mu, \varphi)}$ is not affected in the last two derivatives, the only non-zero blocks in $d$ are the derivative with respect to the real and imaginary parts of mode shape $\varphi_j$ in $\mathcal{J}_{(\lambda, \varphi)}^{\zeta}$ multiplied by the real and imaginary part of $\alpha_{ij} \varphi_j$. Plugging in (B.13), the block of $(I_{m^*} \otimes S(\theta_0)^T) \mathcal{O}'_c(\eta_0)$ for mode shape $\varphi_j$ gets multiplied by $\alpha_{ij} \varphi_j$, which is the only non-zero block in this part of the product and writes in the respective complex-valued way as

$$ S(\theta_0)^T (\Lambda_j \otimes I_r) \alpha_{ij} \varphi_j = S(\theta_0)^T \begin{bmatrix} I_r \\ \lambda_j I_r \\ \vdots \\ \lambda_j^p I_r \end{bmatrix} \alpha_{ij} \varphi_j = \alpha_{ij} S(\theta_0)^T \begin{bmatrix} \varphi_j \\ \lambda_j \varphi_j \\ \vdots \\ \lambda_j^p \varphi_j \end{bmatrix} = 0. $$

Note that the vector in the last term is the $j$-th column of $\mathcal{O}_c(\eta_0)$, hence the product is zero since $S(\theta_0)^T \mathcal{O}_c(\eta_0) = 0$ (cf. [24, Proof of Theorem 7(a)]). Thus, $d=0$ and $\mathcal{J} = \mathcal{J}_{modal}\mathcal{J}_{FE} = \mathcal{J}_{modal}\mathcal{J}_{FE}^{\alpha}$, independently of the particular solution of the mode shape sensitivity.

Appendix C. Invariance of the tests to column scaling of the sensitivity matrix

If the $i$-th column of the sensitivity matrix $\mathcal{J}_i$ is scaled with some factor $\alpha \in \mathbb{R}$, then the direct test (12) writes

$$ \begin{aligned} \tilde{t}_{\text{dir}}^{\text{def}} &= \zeta^T \Sigma^{-1} \alpha \mathcal{J}_i (\alpha \mathcal{J}_i^T \Sigma^{-1} \alpha \mathcal{J}_i)^{-1} \alpha \mathcal{J}_i^T \Sigma^{-1} \zeta \\ &= \zeta^T \Sigma^{-1} \mathcal{J}_i (\mathcal{J}_i^T \Sigma^{-1} \mathcal{J}_i)^{-1} \mathcal{J}_i^T \Sigma^{-1} \zeta = t_{\text{dir}}^i, \end{aligned} $$

hence a scaling of the columns of the sensitivity matrix does not affect the direct test.

The minmax test computation involves the entire sensitivity matrix. Suppose column $\mathcal{J}_i$ corresponding to the tested element is scaled with factor $\alpha$, and the columns of $\mathcal{J}_i$ are scaled as $\mathcal{J}_i D$ with a diagonal invertible matrix $D$. Then, the Fisher information (11) writes with the scaled sensitivities

$$ \begin{bmatrix} \tilde{F}_{i,i} & \tilde{F}_{i,\tilde{i}} \\ \tilde{F}_{\tilde{i},i} & \tilde{F}_{\tilde{i},\tilde{i}} \end{bmatrix} = \begin{bmatrix} \alpha^2 F_{i,i} & \alpha F_{i,\tilde{i}} D \\ \alpha D F_{\tilde{i},i} & D F_{\tilde{i},\tilde{i}} D \end{bmatrix}. $$

The partial residuals defined in Section 4.2.2 write $\tilde{\zeta}_i = \alpha \mathcal{J}_i^T \Sigma^{-1} \zeta = \alpha \zeta_i$, $\tilde{\zeta}_i = D \mathcal{J}_i^T \Sigma^{-1} \zeta = D \zeta_{\tilde{i}}$ and the robust residual writes $\tilde{\zeta}_i^* = \tilde{\zeta}_i - \tilde{F}_{i,\tilde{i}} \tilde{F}_{\tilde{i},\tilde{i}}^{-1} \tilde{\zeta}_{\tilde{i}} = \alpha \zeta_i - \alpha F_{i,\tilde{i}} D (D F_{\tilde{i},\tilde{i}} D)^{-1} D \zeta_{\tilde{i}} = \alpha \zeta_i^*$. The projected Fisher information in Section 4.2.2 writes $\tilde{F}_i^* = \tilde{F}_{i,i} - \tilde{F}_{i,\tilde{i}} \tilde{F}_{\tilde{i},\tilde{i}}^{-1} \tilde{F}_{\tilde{i},i} = \alpha^2 F_i^*$, and finally the minmax test (13) writes

$$ \begin{aligned} \tilde{t}_{\text{mm}}^i &= \tilde{\zeta}_i^{*T} \tilde{F}_i^{*-1} \tilde{\zeta}_i^* \\ &= (\alpha \tilde{\zeta}_i^*)^T (\alpha^2 F_i^*)^{-1} \alpha \tilde{\zeta}_i^* = t_{\text{mm}}^i, \end{aligned} $$

hence a scaling of the columns of the sensitivity matrix does not affect the minmax test either.
References

[1] C. Farrar, K. Worden, An introduction to structural health monitoring, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365 (1851) (2007) 303–315.

[2] C. Farrar, S. Doebling, D. Nix, Vibration-based structural damage identification, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Science 359 (1778) (2001) 131–149.

[3] E. Carden, P. Fanning, Vibration based condition monitoring: a review, Structural Health Monitoring 3 (4) (2004) 355–377.

[4] W. Fan, P. Qiao, Vibration-based damage identification methods: a review and comparative study, Structural Health Monitoring 10 (1) (2011) 83–111.

[5] J. Brownjohn, Structural health monitoring of civil infrastructure, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 365 (1851) (2007) 589–622.

[6] K. Worden, G. Manson, N. Fieller, Damage detection using outlier analysis, Journal of Sound and Vibration 229 (3) (2000) 647–667.

[7] A. Deraemaeker, E. Reynders, G. De Roeck, J. Kullaa, Vibration-based structural health monitoring using output-only measurements under changing environment, Mechanical Systems and Signal Processing 22 (1) (2008) 34–56.

[8] G. Comanducci, F. Magalhães, F. Ubertini, Á. Cunha, On vibration-based damage detection by multivariate statistical techniques: Application to a long-span arch bridge, Structural Health Monitoring 15 (5) (2016) 505–524.

[9] J. M. W. Brownjohn, P.-Q. Xia, H. Hao, Y. Xia, Civil structure condition assessment by FE model updating: methodology and case studies, Finite Elements in Analysis and Design 37 (10) (2001) 761–775.

[10] M. Friswell, Damage identification using inverse methods, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365 (1851) (2007) 393–410.

[11] E. Simoen, G. De Roeck, G. Lombaert, Dealing with uncertainty in model updating for damage assessment: A review, Mechanical Systems and Signal Processing 56 (2015) 123–149.

[12] M. Limongelli, E. Chatzi, M. Döhler, G. Lombaert, E. Reynders, Towards extraction of vibration-based damage indicators, in: Proc. 8th European Workshop on Structural Health Monitoring, Bilbao, Spain, 2016.

[13] D. Bernal, Load vectors for damage localization, Journal of Engineering Mechanics 128 (1) (2002) 7–14.

[14] D. Bernal, Load vectors for damage location in systems identified from operational loads, Journal of Engineering Mechanics 136 (1) (2010) 31–39.

[15] L. Marin, M. Döhler, D. Bernal, L. Mevel, Robust statistical damage localization with stochastic load vectors, Structural Control and Health Monitoring 22 (3) (2015) 557–573.

[16] M. Bhuyan, M. Döhler, Y. Lecieux, L. Mevel, F. Schoefs, Statistical damage localization with stochastic load vectors using multiple mode sets, Structural Health Monitoring 16 (5) (2017) 518–535.

[17] D. Bernal, M. D. Ulriksen, Subspace exclusion zones for damage localization, Mechanical Systems and Signal Processing 114 (2019) 120–127.

[18] M. Basseville, L. Mevel, M. Goursat, Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios, Journal of Sound and Vibration 275 (3) (2004) 769–794.

[19] E. Balmès, M. Basseville, L. Mevel, H. Nasser, W. Zhou, Statistical model-based damage localization: a combined subspace-based and substructuring approach, Structural Control and Health Monitoring 15 (6) (2008) 857–875.

[20] A. Benveniste, M. Basseville, G. Moustakides, The asymptotic local approach to change detection and model validation, IEEE Transactions on Automatic Control 32 (7) (1987) 583–592.

[21] M. Döhler, L. Mevel, Q. Zhang, Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis, Annual Reviews in Control 42 (2016) 244–256.

[22] M. Basseville, M. Abdelghani, A. Benveniste, Subspace-based fault detection algorithms for vibration monitoring, Automatica 36 (1) (2000) 101–109.

[23] E. Balmès, M. Basseville, F. Bourquin, L. Mevel, H. Nasser, F. Treyssède, Merging sensor data from multiple temperature scenarios for vibration-based monitoring of civil structures, Structural Health Monitoring 7 (2) (2008) 129–142.

[24] M. Döhler, L. Mevel, Subspace-based fault detection robust to changes in the noise covariances, Automatica 49 (9) (2013) 2734–2743.

[25] M. Döhler, L. Mevel, F. Hille, Subspace-based damage detection under changes in the ambient excitation statistics, Mechanical Systems and Signal Processing 45 (1) (2014) 207–224.

[26] M. Döhler, F. Hille, L. Mevel, W. Rücker, Structural health monitoring with statistical methods during progressive damage test of S101 Bridge, Engineering Structures 69 (2014) 183–193.

[27] C. Ventura, P. Andersen, L. Mevel, M. Döhler, Structural Health Monitoring of the Pitt River Bridge in British Columbia, Canada, in: Proc. 6th World Conference on Structural Control and Monitoring, Barcelona, Spain, 2014.

[28] S. Gres, P. Andersen, R. J. Johansen, M. D. Ulriksen, L. Damkilde, A comparison of damage detection methods applied to civil engineering structures, in: International Conference on Experimental Vibration Analysis for Civil Engineering Structures, Springer, 2017, pp. 306–316.

[29] M. Döhler, F. Hille, L. Mevel, Vibration-based monitoring of civil structures with subspace-based damage detection, in: E. Ottaviano, A. Pelliccio, V. Gattulli (Eds.), Mechatronics for Cultural Heritage and Civil Engineering, Springer International Publishing, Cham, 2018, pp. 307–326.

[30] R. O. Duda, P. E. Hart, D. G. Stork, Pattern Classification, John Wiley & Sons, 2012.

[31] C. E. Ventura, J.-F. Lord, M. Turek, A. M. Sereci, D. Radulescu, C. Radulescu, Experimental studies and remote monitoring of IASC-ASCE benchmark test frame, in: Proc. 21st International Modal Analysis Conference (IMAC), 2003.

[32] M. Bhuyan, E. Viefhues, M. Döhler, Y. Lecieux, L. Mevel, F. Hille, F. Schoefs, Output-only subspace and transfer matrix-based damage localization and quantification, in: Proc. 7th International Operational Modal Analysis Conference, Ingolstadt, Germany, 2017.

[33] J.-N. Juang, Applied system identification, Prentice Hall, Englewood Cliffs, NJ, USA, 1994.

[34] A. Benveniste, L. Mevel, Nonstationary consistency of subspace methods, IEEE Transactions on Automatic Control AC-52 (6) (2007) 974–984.
[35] P. Van Overschee, B. De Moor, Subspace Identification for Linear Systems: Theory, Implementation, Applications, Kluwer, 1996.

[36] B. Peeters, G. De Roeck, Reference-based stochastic subspace identification for output-only modal analysis, Mechanical Systems and Signal Processing 13 (6) (1999) 855–878.

[37] M. Döhler, L. Mevel, Fast multi-order computation of system matrices in subspace-based system identification, Control Engineering Practice 20 (9) (2012) 882–894.

[38] A.-M. Yan, J.-C. Golinval, Null subspace-based damage detection of structures using vibration measurements, Mechanical Systems and Signal Processing 20 (3) (2006) 611–626.

[39] S. Gres, M. D. Ulriksen, M. Döhler, R. J. Johansen, P. Andersen, L. Damkilde, S. A. Nielsen, Statistical methods for damage detection applied to civil structures, in: Proc. 10th International Conference on Structural Dynamics (EURODYN), Rome, Italy, 2017.

[40] E. Parloo, P. Guillaume, M. Van Overmeire, Damage assessment using mode shape sensitivities, Mechanical Systems and Signal Processing 17 (3) (2003) 499–518.

[41] E. M. Hernandez, Identification of isolated structural damage from incomplete spectrum information using 11-norm minimization, Mechanical Systems and Signal Processing 46 (1) (2014) 59–69.

[42] M. Basseville, Information criteria for residual generation and fault detection and isolation, Automatica 33 (5) (1997) 783–803.

[43] W. Heylen, S. Lammens, P. Sas, Modal Analysis Theory and Testing, Katholieke Universiteit Leuven, Belgium, 1998.