<figure><img src="image_1.png" /></figure>

HAL
open science

# Subspace-based damage detection in engineering structures considering reference uncertainties and temperature effects

Eva Viefhues

► To cite this version:

Eva Viefhues. Subspace-based damage detection in engineering structures considering reference uncertainties and temperature effects. Dynamique, vibrations. Université de Rennes 1, 2021. English.
NNT: 2021REN1S116 . tel-04249624

HAL Id: tel-04249624

https://inria.hal.science/tel-04249624

Submitted on 19 Oct 2023

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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# THÈSE DE DOCTORAT DE

L'UNIVERSITÉ DE RENNES 1

ÉCOLE DOCTORALE N° 601

Mathématiques et Sciences et Technologies
de l'Information et de la Communication

Spécialité : signal, image, vision

Par

Eva VIEFHUES

## Subspace-based damage detection in engineering structures considering reference uncertainties and temperature effects

Thèse présentée et soutenue à Rennes, le 14/12/2021

Unité de recherche : Inria Rennes

## Rapporteurs avant soutenance :

Peter KRAEMER Professeur, University of Siegen, Allemagne
Julien WAEYTENS Chargé de recherche, Université Gustave Eiffel

## Composition du Jury :

<table><tr><td>Président :</td><td>Maria Pina LIMONGELLI</td><td>Professeur, Politecnico di Milano, Italie</td></tr><tr><td>Rapporteurs :</td><td>Peter KRAEMER</td><td>Professeur, University of Siegen, Allemagne</td></tr><tr><td></td><td>Julien WAEYTENS</td><td>Chargé de recherche, Université Gustave Eiffel</td></tr><tr><td>Examinateurs :</td><td>Michael DÖHLER</td><td>Chargé de recherche, Inria Rennes</td></tr><tr><td></td><td>Falk HILLE</td><td>Chercheur, BAM Federal Institute for Materials Research and Testing, Allemagne</td></tr><tr><td>Dir. de thèse :</td><td>Laurent MEVEL</td><td>Directeur de recherche, Inria Rennes</td></tr></table>
<figure><img src="image_3.png" /></figure>
ACKNOWLEDGEMENT

This thesis was written during my scientific work in the Division *Buildings and Structures* at the Federal Institute for Materials Research and Testing (BAM) and in very close cooperation with the research institute Inria Rennes (France). I would like to express my sincere thanks for the financial support from BAM - especially to my superiors Matthias Baeßler and Andreas Rogge, who offered me the opportunity to do research in this exciting field. Gratitude is also owed to the DAAD for the mobility grant within the Procope program making two research trips to Inria Rennes possible.

I am particularly grateful to my supervisor Laurent Mevel from the I4S team at Inria Rennes. In many fruitful discussions I could learn from his comprehensive knowledge in the field of statistics which enabled me to do research in this borderline area between statistics and civil engineering. To my co-supervisor Michael Döhler from the I4S team I would like to express my heartfelt thanks for his great support and mentoring. He not only inspired the research proposal, but also provided guidance through the challenges of the research process and was open for discussions at any time. The very frequent and reliable exchange with him has enabled me to take a critical yet constructive look at my results. With his wide knowledge in the research field, his ability to make complex things simple, and his openness and positivity it was always a great pleasure to work with him. Thanks also go to Julien Waeytens and Yann Lecieux who, as jury members of the CSID, encouraged me in my development as a PhD student.

Furthermore, I would like to thank my colleagues at BAM. Especially my supervisor Falk Hille helped me in many constructive discussions with his great experience to always look at the problems and results also from the user's point of view. I would like to thank Fred Ziegler for his support in implementing the experiments on the test bridge. Likewise, I am very thankful to Patrick Simon for putting up with extra work during his experiments in the climate chamber to make it possible to provide experimental data for my purposes. I thank my colleagues in the SHM group for the interesting technical and personal exchange. It was a great pleasure to be able to learn from their experiences, some of them of many years, and to get involved in some of their exciting projects. I also enjoyed the many friendly conversations with the other colleagues and especially with my roommate Borana
Kullolli. I would like to thank Simone Morawietz and Margrit Kaiser for their great organizational support.

I am very grateful to my reviewers, Peter Kraemer from the University of Siegen and Julien Waeytens from the University Gustave Eiffel, for their high interest in my work. Through their constructive comments also in the run-up to the defense and in the exciting discussions about the results, I was able to learn further from their great experience. I would like to thank them and Maria Pina Limongelli very much for taking the time to be the jury for my thesis defense. It is a great pleasure to have had the opportunity to meet them.

Finally, I would like to thank my family and friends for their support and words of encouragement. A special thanks is due to my parents who have always encouraged me in my life choices. Their personal commitment, especially in the last two years during the pandemic, made it possible for me to finish this work. Likewise, I would like to extend a very special thank you to my husband and children. Without their tireless support, patience, understanding and cheerfulness, I would never have been able to complete this work.
<figure><img src="image_4.png" /></figure>
# TABLE OF CONTENTS

<table><tr><td>Table of Contents</td><td>iii</td></tr><tr><td>List of figures</td><td>vii</td></tr><tr><td>List of tables</td><td>x</td></tr><tr><td>Notation</td><td>xiii</td></tr><tr><td>Introduction</td><td>1</td></tr><tr><td>Context of the work</td><td>1</td></tr><tr><td>Chosen method</td><td>2</td></tr><tr><td>Contribution to the research field</td><td>3</td></tr><tr><td>Outline of the work</td><td>5</td></tr><tr><td>I Preliminaries</td><td>7</td></tr><tr><td>1 State of the art</td><td>9</td></tr><tr><td>1.1 Introduction</td><td>9</td></tr><tr><td>1.2 Vibration-based damage diagnosis concepts</td><td>10</td></tr><tr><td>1.3 Feature evaluation and related uncertainties</td><td>16</td></tr><tr><td>1.4 Data-based damage detection under changing environmental conditions</td><td>20</td></tr><tr><td>1.4.1 Methods with relation to measured environmental parameters</td><td>21</td></tr><tr><td>1.4.2 Methods with unmeasured environmental conditions</td><td>22</td></tr><tr><td>2 Background theory</td><td>27</td></tr><tr><td>2.1 Introduction</td><td>27</td></tr><tr><td>2.2 Structural dynamics of output-only systems</td><td>28</td></tr><tr><td>2.3 System identification</td><td>29</td></tr><tr><td>2.3.1 Concept of stochastic subspace identification</td><td>29</td></tr><tr><td>2.3.2 Covariance-driven subspace identification</td><td>31</td></tr></table>
TABLE OF CONTENTS

<table><tr><td>2.3.3</td><td>Identification of the eigenstructure</td><td>31</td></tr><tr><td>2.4</td><td>Stochastic subspace-based damage detection</td><td>32</td></tr><tr><td>2.4.1</td><td>Subspace-based residual function</td><td>33</td></tr><tr><td>2.4.2</td><td>Statistical damage testing</td><td>34</td></tr><tr><td>2.4.3</td><td>Asymptotic residual covariance</td><td>36</td></tr><tr><td>2.4.4</td><td>Decision making from test results</td><td>37</td></tr><tr><td>2.5</td><td>Considerations for uncertainty computation</td><td>39</td></tr><tr><td>2.5.1</td><td>Mathematical remarks</td><td>39</td></tr><tr><td>2.5.2</td><td>Concept of uncertainty propagation</td><td>40</td></tr><tr><td>2.5.3</td><td>Relevant uncertainty computations from literature</td><td>41</td></tr><tr><td>2.6</td><td>Model interpolation of linear parameter varying systems</td><td>44</td></tr><tr><td>2.6.1</td><td>Concept</td><td>45</td></tr><tr><td>2.6.2</td><td>Model interpolation</td><td>45</td></tr><tr><td>2.7</td><td>Summary</td><td>46</td></tr><tr><td>II</td><td>Contribution to theory</td><td>47</td></tr><tr><td>3</td><td>Stochastic subspace-based damage detection with uncertain reference null space</td><td>49</td></tr><tr><td>3.1</td><td>Introduction</td><td>49</td></tr><tr><td>3.2</td><td>Damage detection with uncertain reference</td><td>50</td></tr><tr><td>3.2.1</td><td>Formalized problem statement</td><td>50</td></tr><tr><td>3.2.2</td><td>Asymptotic distribution of the residual</td><td>52</td></tr><tr><td>3.2.3</td><td>Asymptotic covariance of the residual</td><td>54</td></tr><tr><td>3.2.4</td><td>Test statistic</td><td>58</td></tr><tr><td>3.3</td><td>Proof of concept and discussion</td><td>59</td></tr><tr><td>3.3.1</td><td>Study setup</td><td>60</td></tr><tr><td>3.3.2</td><td>Statistical properties and test performance</td><td>61</td></tr><tr><td>3.3.3</td><td>Data length in the reference state</td><td>63</td></tr><tr><td>3.3.4</td><td>Data length of the testing data set</td><td>66</td></tr><tr><td>3.3.5</td><td>System parametrization</td><td>68</td></tr><tr><td>3.3.6</td><td>Reference null space computed from observability matrix</td><td>70</td></tr><tr><td>3.4</td><td>Summary</td><td>71</td></tr><tr><td>3.5</td><td>Dissemination</td><td>72</td></tr></table>
4 Stochastic subspace-based damaged detection under changing environmental conditions 75

<table><tr><td>4.1</td><td>Introduction</td><td>75</td></tr><tr><td>4.2</td><td>Damage detection robust to temperature effects</td><td>76</td></tr><tr><td>4.2.1</td><td>Formalized problem statement</td><td>76</td></tr><tr><td>4.2.2</td><td>Parameter adapted reference null space</td><td>79</td></tr><tr><td>4.2.3</td><td>Residual function and statistical evaluation</td><td>83</td></tr><tr><td>4.2.4</td><td>Residual covariance</td><td>84</td></tr><tr><td>4.2.5</td><td>Discussion and comparison to previous approaches</td><td>86</td></tr><tr><td>4.3</td><td>Proof of concept and discussion</td><td>87</td></tr><tr><td>4.3.1</td><td>Study setup</td><td>88</td></tr><tr><td>4.3.2</td><td>Damage detection under changing temperatures</td><td>89</td></tr><tr><td>4.3.3</td><td>Covariance computation from test data</td><td>93</td></tr><tr><td>4.3.4</td><td>Considerations for reference setup</td><td>94</td></tr><tr><td>4.3.5</td><td>Effects of changes in the excitation covariance</td><td>96</td></tr><tr><td>4.4</td><td>Summary</td><td>97</td></tr><tr><td>4.5</td><td>Dissemination</td><td>98</td></tr></table>

III Applications 101

5 Damage detection of a steel frame considering uncertainties in the reference matrix 103

<table><tr><td>5.1</td><td>Introduction</td><td>103</td></tr><tr><td>5.2</td><td>Experimental setup and dynamical properties</td><td>104</td></tr><tr><td>5.3</td><td>Damage detection with uncertain reference null space</td><td>106</td></tr><tr><td>5.4</td><td>Summary</td><td>109</td></tr><tr><td>5.5</td><td>Dissemination</td><td>109</td></tr></table>

6 Damage detection of a temperature affected laboratory beam structure 111

<table><tr><td>6.1</td><td>Introduction</td><td>111</td></tr><tr><td>6.2</td><td>Experimental setup and dynamical properties</td><td>112</td></tr><tr><td>6.3</td><td>Damage detection with temperature effects</td><td>113</td></tr><tr><td>6.4</td><td>Summary</td><td>119</td></tr><tr><td>6.5</td><td>Dissemination</td><td>119</td></tr></table>
# 7 Damage detection of a testing bridge exposed to ambient temperature variations 121

7.1 Introduction 121

7.2 Experimental setup and dynamical properties of the structure 122

7.3 Damage detection with temperature effects 125

7.3.1 General settings and temperature effect 125

7.3.2 Robust damage detection 127

7.3.3 Effect of reference setup 128

7.4 Summary 131

7.5 Dissemination 132

# Conclusion 133

Considered aspects 133

Summary, improvements and limits 133

Outlook 136

# Resume in French 139

# Bibliography 151
LIST OF FIGURES

<table><tr><td>1</td><td>Westend bridge and excitation device of EMPA [Roh+00b].</td><td>1</td></tr><tr><td>1.1</td><td>Damage detection procedure.</td><td>11</td></tr><tr><td>1.2</td><td>Classification concept.</td><td>12</td></tr><tr><td>2.1</td><td>χ<sup>2</sup>-distribution for different parameters (left) and definition of threshold (right).</td><td>38</td></tr><tr><td>2.2</td><td>Typical receiver operating characteristic curve.</td><td>39</td></tr><tr><td>3.1</td><td>Mass-spring chain with four sensors.</td><td>60</td></tr><tr><td>3.2</td><td>Theoretical χ<sup>2</sup>-distribution and histograms of the new test statistics and the conventional test with Σ<sub>1</sub> = 0, M = 50,000, N = 100,000.</td><td>61</td></tr><tr><td>3.3</td><td>Histogram in the reference and the damaged state of the conventional (left) and the new test (right), M = 50,000, N = 100,000, damage = 2%.</td><td>62</td></tr><tr><td>3.4</td><td>POD (left) and ROC diagram with damage of 2% (right) of the conventional and the new test, M = 50,000, N = 100,000.</td><td>63</td></tr><tr><td>3.5</td><td>Empirical mean of conventional and new tests for different data lengths M, N = 100,000.</td><td>64</td></tr><tr><td>3.6</td><td>Histogram in the reference and the damaged state of the conventional (left) and the new test (right) for different M, N = 100,000, damage = 2%.</td><td>65</td></tr><tr><td>3.7</td><td>POD of conventional and new test for different M, N = 100,000, damage = 1.4%.</td><td>65</td></tr><tr><td>3.8</td><td>POD (left) and ROC diagram with damage of 2% (right) of the conventional and the new test, M = 500,000, N = 100,000.</td><td>66</td></tr><tr><td>3.9</td><td>Histogram in the reference and the damaged state of the conventional (left) and the new test (right) for different N, M = 100,000, damage = 2%.</td><td>67</td></tr><tr><td>3.10</td><td>POD of conventional and new test for different N, M = 100,000, damage = 2%.</td><td>67</td></tr></table>
3.11 Empirical mean of conventional and new tests for different data lengths M for a system parametrized with respect to structural parameters, $N = 100,000$. 69

3.12 POD of conventional and new test for a system parametrized with respect to structural parameters for $M = 50,000$ (left) and $M = 500,000$ (right), $N = 100,000$. 69

3.13 POD of conventional and new test with non parametric test for $M = 50,000$ (left) and $M = 500,000$ (right), $N = 100,000$. 69

3.14 Empirical mean of conventional and new tests for different data lengths M when the reference null space is computed from observability matrix, $N = 100,000$. 70

3.15 POD of conventional and new test when the reference null space is computed from observability matrix for $M = 50,000$ (left) and $M = 500,000$ (right), $N = 100,000$. 71

4.1 Determination of weights at parameter working points for any $P$ from a Gaussian function. 80

4.2 Mass-spring chain and temperature dependent stiffness. 88

4.3 Test values at different testing temperatures computed from a single reference data set at 10°C for the whole temperature range (left) and with zoom on the temperatures close to the reference temperature (right). 89

4.4 Test values at different testing temperatures computed with the *average approach* [Bal+08a]. 90

4.5 Test values at different testing temperatures computed with the three strategies of the new method. 91

4.6 Test performance of the average method and the three strategies of the new method. Probability of detection (left) and ROC curves for 0.8 % stiffness reduction (right). 92

4.7 Test values at different testing temperatures when the uncertainties of the reference null space are neglected (left) and ROC curves for 0.8 % stiffness reduction computed with S(H) strategy. 93

4.8 Test values at different testing temperatures when the covariance is computed from the testing data (left) and ROC curves for 0.8 % stiffness reduction computed with S(H) strategy. 94

4.9 Test values at different testing temperatures for different reference setups. 96
4.10 Test values when reference and testing data sets have different excitation covariance, computed with S(H) interpolation (left) and S(O) interpolation (right) strategies. 97

5.1 Laboratory testing structure (left) and sensor positions (right). 104

5.2 Reversible damage detail (left) and screw configuration (right). 105

5.3 Test values of the conventional (left) and the new test (right) in the undamaged reference state computed with different reference data lengths *M*. 107

5.4 Test values of the conventional (left) and the new test (right) in the reference and at different damaged states. 107

5.5 Test values of the conventional (left) and the new test (right) in the reference and at different damaged states for different *N*. 108

6.1 Test setup of the concrete beam with locations of geophones, shaker, and loading for the damage introduction. 112

6.2 Concrete beam in the climate chamber. 113

6.3 1st to 4th eigenfrequencies at different system states and at different temperatures. The y-axis is scaled to about 5 % of each frequency. 114

6.4 Fixed reference null space and covariance matrix. Test values at different testing temperatures in the undamaged and damaged states. Reference data at 5°C. 116

6.5 Averaged reference null space and covariance matrix. Test values at different testing temperatures in the undamaged and damaged states. Reference data at -25°C, 5°C, and 40°C. 116

6.6 Interpolated reference null space and covariance matrix. Test values at different testing temperatures in the undamaged and damaged states. Reference data at -25°C, 5°C, and 40°C. 116

6.7 Effect of weighting function in interpolated reference null space and covariance matrix. Sharp weighting function (top) and wide weighting function (bottom). Test values at different testing temperatures in the undamaged and damaged states. Reference data at -25°C, 5°C, and 40°C. 118

7.1 Testing bridge. 122

7.2 Static system of the testing bridge. 123
<table>
  <tbody>
    <tr>
      <td>7.3</td>
      <td>Pre-stressing facility (left) and encased shaker (right).</td>
      <td>123</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Test setup testing bridge with geophones and shaker for excitation.</td>
      <td>124</td>
    </tr>
    <tr>
      <td>7.5</td>
      <td>Stabilization diagram at 15°C.</td>
      <td>125</td>
    </tr>
    <tr>
      <td>7.6</td>
      <td>First four modeshapes of the testing bridge.</td>
      <td>125</td>
    </tr>
    <tr>
      <td>7.7</td>
      <td>Temperature effect to the first four eigenfrequencies of the testing bridge.</td>
      <td>126</td>
    </tr>
    <tr>
      <td>7.8</td>
      <td>Test values with reference data from one fixed reference temperature.</td>
      <td>127</td>
    </tr>
    <tr>
      <td>7.9</td>
      <td>Test values of the 4 algorithms for different testing temperatures in the undamaged and the damaged case, computed for reference setup A.</td>
      <td>128</td>
    </tr>
    <tr>
      <td>7.10</td>
      <td>Test values of the 4 algorithms for different testing temperatures in the undamaged case, computed for reference setup B, C, D and E from top to bottom.</td>
      <td>129</td>
    </tr>
    <tr>
      <td>7.11</td>
      <td>Weights for different reference setups with weighting function W1 (left) and with optimized weighting functions (right).</td>
      <td>130</td>
    </tr>
    <tr>
      <td>7.12</td>
      <td>Effect of optimized weighting function to false positives for the interpolation approach in reference setups C, D, and E.</td>
      <td>131</td>
    </tr>
    <tr>
      <td>1</td>
      <td>Histogrammes du test dans les états de référence et endommagé, pour le test conventionnel (gauche) et le nouveau test (droite).</td>
      <td>144</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Probabilité de détection (gauche) and dépendance du test de la longueur de données M dans l'état de référence (droite).</td>
      <td>144</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Valeurs du test pour toutes les strategies développées, ainsi que pour l'approche existante [Bal+08a] aux températures de test variées et les références prises à 0, 5, 10, 15 and 20°C.</td>
      <td>145</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Valeurs du test conventionnel (gauche) et du nouveau test (droite) dans l'état de référence et plusieurs états endommagés.</td>
      <td>146</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Valeurs du test avec la stratégie d'interpolation des matrices Hankel (en haut), et avec l'approche existante [Bal+08a] (en bas). Les données de référence ont été utilisées à -25°C, 5°C, et 40°C.</td>
      <td>147</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Testing bridge.</td>
      <td>148</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Test values of the 4 algorithms for different testing temperatures in the undamaged and the damaged case, computed for reference setup A.</td>
      <td>148</td>
    </tr>
  </tbody>
</table>
LIST OF TABLES

<table><tr><td>4.1 Reference setups</td><td>95</td></tr><tr><td>5.1 Frequencies of the steel frame in Hz from [Hil18]</td><td>106</td></tr><tr><td>7.1 Frequencies and damping ratios at 15°C.</td><td>125</td></tr><tr><td>7.2 Reference setups for testing bridge</td><td>128</td></tr></table>
<figure><img src="image_5.png" /></figure>
List of variables

A state transition matrix

b bias

C observation matrix

C damping matrix

C controllability matrix

c ratio of N and M

D₁ diagonal matrix of n singular values of a matrix

D₂ diagonal matrix of noise singular values very close to zero of a matrix

E selection matrix

e selection vector

F selection matrix

f eigenfrequency

G cross-covariances between states and output

H Block Hankel matrix, data matrix in general for SSI

I identity matrix

J_{X1,X2} sensitivity of X1 with respect to X2

K stiffness matrix

k stiffness of elements in simulation

L selection matrices position sensor

M number of data samples in reference data

M mass matrix

m DOF of dynamical system

m masses of elements in simulation

N number of data samples in test data

N_b block length for covariance estimation

n model order

n_b number of blocks for covariance estimation
O observability matrix

P scheduling parameter

P permutation matrix

q displacement at DOF

R output covariances

r number of sensors

S reference null space

T temperature

t in the context of system representations: continuous time

t in the context of hypothesis testing: test statistics

U₁ column space of a matrix

U₂ left null space of a matrix

u number of available reference data sets / scheduling parameters P

uᵢ left singular vector of the Hankel matrix

V₁ row space of a matrix

V₂ null space of a matrix

v measurement noise

vᵢ right singular vector of the Hankel matrix

W invertible weighting matrix in SSI

w state noise

x system state

Y set of output data

y output data column vector at one time instance at all sensor locations

Z SSI algorithm depending factor matrix

γ change in the residual

Δ diagonal matrix of Λ

δ parameter change

ε perturbation magnitude in Taylor approximation

ζ subspace-based residual (see also explanation of $\tilde{\mathcal{D}}$)

θ parameter vector, contains e.g. the modal parameters

θ* parameter vector in the reference state

Λ vector of eigenvalues λ

λ in the context of system identification: eigenvalue of A in discrete state-space system
$\lambda$ in the context of hypothesis test: non-centrality parameter

$\mu$ eigenvalue of system in continuous time

$\nu$ in the context of equation of motion: excitation function

$\nu$ in the context of $\chi^2$-distribution: degree of freedom

$\xi$ damping ratio

$\rho$ weighting function for interpolation

$\Sigma_X$ covariance matrix of random variable $X$

$\sigma_i$ singular value of the Hankel matrix

$\sigma^2$ variance of the Gaussian weighting function

$\tau$ time increment

$\Phi$ matrix with eigenvectors $\phi$

$\phi$ eigenvector of $A$ in discrete state-space system

$\varphi$ mode shape

$\chi$ left eigenvector of $A$ in discrete state-space system

$\Psi$ matrix with mode shapes $\varphi$

$\omega$ angular frequency

List of indices

$a$ acceleration

$c$ separated into conjugate complex pairs

$d$ displacement

$i$ counting variable

$j$ counting variable

$k$ time index

$l$ counting variable

$p$ size parameter, e.g. $p+1$ is number of block rows in Hankel matrix $\mathcal{H}$

$q$ size parameter, e.g. number of block columns in Hankel matrix $\mathcal{H}$

$(s)$ counting variable

$s$ number of columns of $S$

$t$ number of rows of $S$

$v$ velocity
List of symbols

$T$ transposed

$\dagger$ pseudoinverse

$\Delta$ perturbation

$\otimes$ Kronecker product

$\mathcal{N}$ Normal distribution

$\text{cov}$ covariance

$\text{vec}$ stacking operator

$\mathbf{E}$ expectation

$\mathbf{Re}$ real part of a complex number or matrix

$\mathbf{Im}$ imaginary part of a complex number or matrix

$\hat{\Box}$ estimate of $\Box$

$\dot{\Box}$ derivation of $\Box$ with respect to time

$\Box'$ derivation of $\Box$ with respect to variable

$\bar{\Box}$ mean of $\Box$

$\Box'$ interpolated version of $\Box$

$\tilde{\Box}$ variable $\Box$ belongs to residual with consideration of uncertainty in the reference null space

List of acronyms

AR autoregressive

ARMA autoregressive moving average

ARX autoregressive model with exogenous inputs

CCP-VARX constant coefficient pooled vector autoregressive with exogenous variables

COMAC coordinate modal assurance criterion

DOF degrees of freedom

ERA eigensystem-realization-algorithm

EMPA Eidgenössische Materialprüfungs- und Forschungsanstalt

FE finite element

FP-VARX functionally pooled vector autoregressive with exogenous variables

FRF frequency response function

GMM Gaussian mixture model

GLR generalized likelihood ratio
<table><tr><td>IDDM</td><td>interpolated damage detection method</td></tr><tr><td>LPV</td><td>linear parameter varying system</td></tr><tr><td>LSWT</td><td>lag shifted whiteness test</td></tr><tr><td>LTI</td><td>linear time-invariant</td></tr><tr><td>MAC</td><td>modal assurance criterion</td></tr><tr><td>MAF</td><td>moving average filter</td></tr><tr><td>MCD</td><td>minimum covariance determinant</td></tr><tr><td>MD</td><td>Mahalanobis distance</td></tr><tr><td>MEWMA</td><td>multivariate exponentially weighted moving average</td></tr><tr><td>MMSE</td><td>minimum mean square error</td></tr><tr><td>MPC</td><td>modal phase collinearity</td></tr><tr><td>MVEE</td><td>maximum volume enclosing ellipsoid</td></tr><tr><td>PCA</td><td>principle component analysis</td></tr><tr><td>PCE</td><td>polynomial chaos expansion</td></tr><tr><td>PFA</td><td>probability of false alarms</td></tr><tr><td>POD</td><td>probability of detection</td></tr><tr><td>PSD</td><td>power spectral density</td></tr><tr><td>ROC</td><td>receiver operating characteristic</td></tr><tr><td>SDDLV</td><td>stochastic dynamic damage locating vector</td></tr><tr><td>SHM</td><td>structural health monitoring</td></tr><tr><td>SSDD</td><td>stochastic subspace-based damage detection</td></tr><tr><td>SSI</td><td>stochastic subspace identification</td></tr><tr><td>SVD</td><td>singular value decomposition</td></tr></table>
<figure><img src="image_6.png" /></figure>
Context of the work

Structural health monitoring (SHM) can be understood as the strategy to provide peri-
odically updated information on the condition or damage state of a mechanical structure.
Since the 1990s research on vibration-based methods for automated damage diagnosis of
engineering structures increased rapidly. Some effort was made in continuous monitor-
ing strategies, since many structures reached the end of their calculated life cycle and
construction methods from the 1960s with insufficient evaluation of their critical details
led to deteriorating conditions. One example is the long-term monitoring of the highway
Westend bridge in Berlin (Germany) in Fig. 1 that was started in 2000 by the German
Federal Institute for Materials Research and Testing (BAM) and is still ongoing [BHS16;
Hu+15].

The vibration-based methods are inspired by the developments in condition monitor-
ing from aerospace and offshore oil industries. The principle is based on the fact that
the dynamical behavior of a structure is affected by damage altering its stiffness, mass

<figure><img src="image_7.png" /><figcaption>Figure 1 – Westend bridge and excitation device of EMPA [Roh+00b].</figcaption></figure>
or damping properties, which is reflected in the measured vibration response of the sys-
tem. Following the Fundamental Axioms of Structural Health Monitoring introduced by
Worden et al. in [Wor+07], the structure condition is monitored by extracting damage
sensitive features from vibration measurements and comparing the feature value in a
(healthy) reference state to the value in the current testing state. Such automated SHM
procedures can complement and partly replace the obligatory visual inspections of engi-
neering structures by possibly continuous information, also at locations that are difficult
to reach, and thus help to guarantee high infrastructure integrity.

However, special conditions have to be considered for the application of vibration-
based damage diagnosis to engineering structures. On one hand, the monitoring setup
must consider rough environmental conditions and poor accessibility requiring new sensor
techniques [Leo06; SSS10]. It must be designed for a limited number of sensors that can be
placed on a large-scale and often complex structure, employing optimal sensor placement
and clustering approaches [Bal+08b; All+19; Men+20].

On the other hand the measured data must be processed with methods which are ade-
quate for engineering structures. Since the structure excitation usually is unmeasured and
always affected by e.g. wind, output-only methods are particularly interesting in this con-
text. They are based on the assumption of ambient excitation, and the dynamic properties
can be evaluated under normal operating conditions and without human interaction. It
has to be taken into account that the damage sensitive features are usually computed
from noisy data and thus are subject to statistical variability. Hence, the statistical prop-
erties of the feature need to be considered properly when comparing the testing data to
the reference data, in order to decide if a change is significant, i.e. if there is damage or
not. Moreover, when significant changes are detected it must be distinguished between
changes due to damage and changes due to variability in the environmental conditions
[Soh07], otherwise high false alarm rates may occur. These requirements to robust data
processing are in the scope of the present work and the developments are based on the
following methodology.

Methodology

Inspired by the covariance-driven stochastic subspace identification (SSI) the stochas-
tic subspace-based damage detection (SSDD) method was introduced in [BAB00]. It is
formulated for the output-only state-space representation of the dynamical system and is
based on estimates of the output covariances. A reference model represented by the left
null subspace of a Hankel matrix built from output covariances is confronted to the cor-
responding data Hankel matrix of the current system in a residual function. The residual
is asymptotically Gaussian distributed and due to orthogonality it has zero mean if the
tested system equals the reference system. If the system has changed the residual will be
different from zero. In a subsequent hypothesis test, it is evaluated if the change is statis-
tically significant, which is done by means of generalized likelihood ratio test considering
the asymptotic local approach [BBM87]. The test takes into account the uncertainties of
the tested data resulting from measurement and excitation noise by means of the resid-
ual covariance matrix. Thus, the statistical properties of the damage detection test can
be characterized precisely in the healthy and the damaged scenarios. The test statistic is
used as damage indicator, where a threshold is derived from empirical studies on reference
data.

Such methods based on direct model-data matching are particularly appealing for au-
tomated damage diagnosis, since the model estimation step is avoided in the testing state.
Successful applications in the field of vibration monitoring are reported e.g. in [JMM12;
DM13b; Döh+14; Ven+14]. Moreover, the approach provides solutions for damage local-
ization and quantification in connection with finite element model-based sensitivities in
the same framework [Bal+08b; DMZ16; All+19]. While having a strong theoretical back-
ground, some robustness under real application is still missing in this framework, which
is the subject of this work.

**Contributions to the research field**

The developments in this thesis aim at robust and automated methods for vibration-
based SHM of engineering structures. The output-only damage detection with SSDD
has been proven to be suitable for this task, since on one hand the dynamical model
building part is performed only in the reference state and avoided during the monitoring
phase, which enables fast performance. On the other hand the statistical properties of
the damage indicator can be described precisely for the reference and the testing system
state. However, some fundamental requirements for robust application are not met in the
existing method, namely regarding the treatment of uncertainties related to the reference
data, as well as the consideration of changes in the environmental variables that may
perturb damage detection. Appropriate developments are in the scope of this work. The
contribution to the research field can be summarized as follows:

**(1) Damage detection with uncertain reference** The stochastic character of the vibration data leads to natural variability in the residual and consequently in the resulting damage indicator. So far in the stochastic subspace-based residual the uncertainties related to the testing data are considered, and it is assumed that the reference is perfectly known without uncertainties. This is true in ideal cases but does not hold in real world applications, since the reference usually is estimated from data and thus can contain considerable uncertainties. When neglecting these uncertainties in the residual evaluation, the observed statistical properties do not match the theoretical expectations and the test behavior becomes unpredictable. Since a correct statistical characterization is missing in this case, thresholds can only be set empirically from studies on reference data.

To close this gap, the statistical properties of the residual function is evaluated for the more realistic case when the reference is estimated from data. A theory-based formulation of the statistical properties is derived leading to a residual covariance, which considers both the uncertainties related to the reference and the testing data. The associated hypothesis test shows very good accordance to the theoretical expected behavior. The results are predictable and independent from data length in the reference and the testing state.

The impact of these developments is not only of academic value. The crucial part of defining a threshold is emphasized repeatedly in the literature, e.g. [BB15; GN09]. With the new residual evaluation for SSDD a threshold now can be derived apriori based on theoretical considerations without extensive empirical studies on reference data. Moreover, neglecting the uncertainty of the reference null space in the previous procedure leads to an inherent bias in the test. This is eliminated with the new formulation leading to a significantly higher test sensitivity.

**(2) Damage detection under changing environmental parameters** A very important part of the monitoring procedure is the formulation of a reference model adequate for the testing situation, in particular when environmental conditions are changing. Many automated vibration-based damage detection methods fail due to such variability, since it affects the dynamical behavior of the structures in a similar way as damages [Soh07]. As a consequence either the damage detection becomes insensitive to damages due to masking effects or high false alarm rates occur.

The basic residual in the SSDD is sensitive to changes in the dynamical system and
so far in this context only few works [Bal+08a; FMG03; FKB13] have been made to consider if these changes are due to damages or result from environmental effects such as temperature. However, these methods are developed empirically. Some of them are not optimal when the temperature in the testing state deviates from the mean temperature of the reference data sets, or they require reference data sets at as many temperature scenarios as possible.

In this work a robust approach is derived based on system theoretic considerations that is inspired by a model interpolation approach for linear parameter varying (LPV) systems from [Zha18], assuming that the dynamical behavior of an engineering structure under changes in an external parameter (temperature) can be approximated as a linear time-invariant system at a fixed temperature point. The new approach provides the possibility to compute a parameter value adapted reference null space for any testing parameter value from data in the undamaged state at only few reference parameter values. Three strategies are derived, which are based on interpolation or regression analysis. The corresponding damage detection test is developed, considering the uncertainties related to the reference null space based on interpolation or regression. The resulting test is shown to be robust and still sensitive to damages.

For the experimental validation a multi-sensor equipped testing bridge is constructed within this thesis on an outside testing field of the Federal Institute for Material Research and Testing Germany (BAM).

## Outline of the work

The present thesis comprises three parts, I Preliminaries, II Contribution to theory, and III Applications.

Part I provides the background for the developments within this work. First, the state of the art is presented in chapter 1 proposing a classification concept for vibration-based damage detection procedures and introducing developments that were made (a) in the context of uncertainty considerations in the damage detection procedure and (b) for methods robust to changes in the environmental conditions. In chapter 2 the background theory is summarized comprising in particular the background of the SSDD and also theory on uncertainty computation and model interpolation.

Part II presents the contributions of this thesis to the methodology. In chapter 3 the damage detection test based on the residual for the SSDD with estimated reference is
derived. The concept is proven in a numerical study of a simple structure and the results
are discussed with respect to the strengths and limits. In chapter 4 the development of
the method robust to changes in environmental conditions is carried out considering three
different strategies for the reference computation. Again, the theoretical part is followed
by a proof of concept on the basis of numerical simulations. The considerations of the
reference setup, from which the reference null spaces can be derived, are discussed.

Finally, the approaches are validated in three application studies in part III. The
behavior of the new SSDD with estimated reference evaluated under stable environmental
conditions is discussed in chapter 5. The damage detection approach robust to temperature
changes is evaluated in chapters 6 and 7 on a concrete beam tested in a climate chamber
as well as on an outdoor testing bridge that was constructed related to this thesis on an
outside testing field and is exposed to realistic environmental conditions.
# PART I

## Preliminaries
In this part the reader is introduced to the subject of the thesis. An overview of the existing research in the field of vibration-based damage detection is given, in particular regarding the consideration of uncertainties and the effects of changing environmental conditions. In a second chapter, the theoretical background of the chosen methods is summarized.
STATE OF THE ART

**1.1 Introduction**

In the past decades many vibration-based structural health monitoring (SHM) tech-
niques have been proposed that aim at automated damage diagnosis procedures in en-
gineering structures. Classically, the procedure comprises the following steps [FDN01]:
the definition of the monitoring task and observation of the structural vibration data,
the extraction of damage sensitive data features, and the evaluation of the features to
classify the structure into undamaged and damaged. The latter is equivalent to the task
of finding feature deterioration and is often considered as a pattern recognition process
[SF01; CF04; FW12].

In the present work developments are made on the feature extraction and evaluation part within SHM. Based on the fact that damage affects the stiffness, mass or damping properties of the structure, and consequently its dynamic properties [FW07], a damage sensitive feature is extracted from the vibration data and evaluated subsequently in four stages with increasing complexity: (1) damage detection, (2) damage localization, (3) damage quantification, and (4) lifetime prediction [Ryt93]. The detection of damage is the fundamental task for SHM, before further levels of the damage diagnosis can be reached. What follows is a summary of the literature mainly related to the damage detection task on linear time-invariant systems that this thesis deals with, although some methods may provide the possibility to be modified for the localization, quantification, or prediction stage. Due to the usually unmeasured excitation of engineering structures, such as wind and traffic, output-only methods are of particular interest.

After a brief overview on vibration-based damage diagnosis concepts and the asso-
ciated damage sensitive features, evaluation techniques and their capability to consider
uncertainties is given. Subsequently methods are collected that are proposed to handle
changes in the environmental conditions.
## 1.2 Vibration-based damage diagnosis concepts

Many overviews on vibration-based damage diagnosis methods can be found that use different classification strategies for the methods, e.g. distinguishing between parametric and non-parametric methods [FK13] with or without machine learning basis [Avc+21], or categories such as basic modal properties, dynamical flexibility, and updating structural model parameters [DFP98; CF04]. Others discriminate between model-based and response-based techniques [Lim19] where the latter are classified in a comprehensive overview in [Lim+21] into modal-based methods, time series methods, and time-frequency methods. The authors in [FW12] propose to furthermore distinguish between feature extraction and feature evaluation step, pointing out the important property of the dimension of the extracted feature, where small dimensions simplify the subsequent evaluation step.

Following the Fundamental Axioms of SHM in [Wor+07] and the data-based decision making in [FKK13], the damage detection procedure can be described as illustrated in Fig. 1.1. First, the feature is extracted from data by some feature extraction operation **FE**. Then, usually an additional metric **ME** is applied with the purpose of quantifying the difference between features in the undamaged reference state and the current monitoring state. Subsequently the comparable quantity is evaluated, where either the quantity itself is considered to be the damage indicator or a statistical evaluation procedure **SE** leads to such an indicator. In the latter case it is possible to take into account uncertainties that are inherent in the damage feature due to the underlying noisy data. For decision making the distribution properties of the damage indicator under undamaged conditions are derived or computed empirically by means of a Monte Carlo study **MC** leading to a threshold to which the damage indicator in the testing state is compared.

This leads to the classification scheme in Fig. 1.2 for the following literature overview listing some prominent categories, but without claim on completeness. The extraction of the feature from data can be interpreted as a compression process of the vibration data with different complexity, similarly to the compression effect in system identification described in [CF04]. There the system model in the time or frequency domain is considered as a compression of the initial time history, and the modal representation is the compressed version of the system model. The feature extraction approaches are directly inspired by the available information at a particular system identification stage. Not all existing methods may fit into the scheme and since the boundaries are blurred some of the presented feature metrics may be meant to be part of the feature extraction or feature evaluation.
<figure><img src="image_8.png" /><figcaption>Figure 1.1 – Damage detection procedure.</figcaption></figure>

**Modal parameters and associated parameters** Initially, damage detection approaches were proposed using modal parameters such as frequency [Ada+78], mode shapes, and modal damping, since they are functions of the physical properties of the structure [DFP98] which are affected by damage. Those features are derived e.g. by non-parametric methods like peak picking in the simplest form or parametric methods like polyreference least-squares complex frequency-domain (LSCF) [Pee+04] or subspace identification methods [VD96; PD99; DM12]. Many approaches and application studies exist that are based on these fundamental dynamical characteristics even in recent publications. A shortcoming of modal parameter based damage detection is the requirement to identify the modal parameters in each analyzed data set, and to track and match them, which can be difficult in the presence of weakly excited or close modes.

An application to wind turbines can be found in [Oli+18] using the difference of the system frequencies and it is emphasized that characteristic situations of such structures can be taken into account. The authors of [RF10] provide a damage detection algorithm that uses frequencies and mode shapes, which are less sensitive to environmental effects [Der+08]. The mode shapes are estimated with an automated modal parameter identification procedure deploying the modal assurance criterion (MAC). The MAC is introduced in [AB82] aiming at a better comparability of mode shapes through a measure of correlation where a value close to 1 shows good correlation. It is pointed out that with the MAC and the modal scale factor the quality of the modal vector can be evaluated [All03]. The MAC is also used in [Ram+10] to discriminate the physical modal parameters of stationary structures. A finite element (FE) model updating method is applied for comparison
<table><thead><tr><th colspan="4">System representation</th></tr><tr><th>Re-<br>presentation</th><th>Signal</th><th>Model in time or frequency domain</th><th>Model in modal domain</th></tr></thead><tbody><tr><th rowspan="2">Available<br>information</th><th>Signal characteristics<br>(non-parametric models)</th><td colspan="2">Mathematical representations of the signal (parametric models)</td></tr><tr><td>Mean & Autocorrelation function,<br>Mean & Power Spectral density</td><td colspan="2">e.g. ARMA, State Space, Kalman filters</td></tr></tbody></table>

<table><caption>Figure 1.2 – Classification concept.</caption><thead><tr><th colspan="7">Damage detection</th></tr><tr><th></th><th>Signal and signal<br>functions</th><th>Characteristics in<br>frequency domain</th><th>Parameters of<br>time series models</th><th>Modal parameters</th><th>Modal parameter<br>based features</th><th></th></tr></thead><tbody><tr><th rowspan="5">Extraction<br>of<br>damage<br>sensitive<br>feature</th><th>Time history</th><td>Frequency Response<br>Function</td><td>Characteristics of<br>autoregressive term</td><td>Frequencies</td><td>Mode shape<br>curvature</td><td></td></tr><tr><th>Covariance function</th><td>Wavelet transforms</td><td>Kalman models /<br>filters</td><td>Mode shapes</td><td>Modal flexibility</td><td></td></tr><tr><th>Autocorrelation and<br>Cross-correlation</th><td>Power Spectral<br>Density</td><td>State-space matrices</td><td>Damping ratios</td><td>Flexibility matrix</td><td></td></tr><tr><th>Signal statistics</th><td>Cross-Spectral<br>Density</td><td></td><td>Ritz vectors</td><td>Modal strain energy</td><td></td></tr><tr><th></th><th colspan="2">Difference</th><th colspan="3">Feature functions / vectors</th><th rowspan="4">(FE) model updating<br>Machine learning concepts</th></tr><tr><th rowspan="2">Feature<br>metric</th><th>Delta</th><th>Kalman filter<br>innovations</th><th>MAC</th><th>modal force<br>vectors</th><th>subspace<br>residual</th></tr><tr><th>Ratio</th><th>Waveform<br>comparison</th><th>COMAC</th><th>statistical load<br>vector</th><th>GMM</th></tr><tr><th>Feature<br>evaluation</th><th colspan="2">Statistical process control (e.g.<br>control charts)</th><th colspan="3">Classical outlier analysis (e.g. MD,<br>MVEE, gap statistics)</th><th></th></tr></tbody></table>

of system states. In [PBS91] a damage detection method is derived that evaluates the changes in curvature modes shapes, since they are known to be more sensitive to damage than the MAC and the coordinate modal assurance criterion (COMAC). The procedure carries the potential to be expanded to damage localization and even quantification.

Damage detection based on modal curvatures is applied e.g. in [Ad99] to detect damage in a pre-stressed concrete bridge introducing the curvature damage factor, which summarizes the difference in curvature mode shape for all modes in one quantity at each sensor point.

Changes in the flexibility matrix estimated from modal parameters are used in [PB94]. Another approach that is based on those changes is derived in [Ber02] and enhanced to the Stochastic Dynamic Damage Locating Vector (SDDLV) for output-only data in [Ber10].
Damaged detection is done by considering the internal forces under a set of load vectors, which are zero at the damaged elements. In [Döh+15] the test performance could be improved by taking the uncertainties of the underlying identified modal parameters into account. The flexibility and mode shape based methods provide the possibility to be used for the damage localization step.

**Parameters of time series models** Features that are derived from the parameters of time series models methods are usually of higher dimension than the basic modal parameters and require some effort to provide comparability between the reference and a testing state. However, the modal identification step is avoided, making the methods more suitable for automated procedures.

One group of such features are based on autoregressive (AR) models such as the AR coefficients or autoregressive moving average (ARMA). In [SCF00] and [FSF01] a one dimensional feature projected from AR model coefficients using principal component analysis (PCA) is proposed, where the projection is mentioned to improve the subsequent evaluation procedure with control charts. In [CB08] time series responses are fitted with ARMA models. The resulting ARMA coefficients are fed to a classifier, which is capable of learning in an unsupervised manner and of forming new classes when the structural response exhibits change. The authors in [NK07] use the vector of AR coefficients and apply the method of Gaussian Mixture Models (GMM) evaluating the significance of a distances with gap statistic. The method provides also information on the damage extend from the Mahalanobis distance (MD) of GMMs and is found to be a useful framework for data fusion.

Features based on Kalman filter [Kál63] are capable to include statistical noise and other inaccuracies. The damage detection in [YDG04] is based on a Kalman filter using stochastic subspace identification (SSI). A prediction error is derived from the filter and testing data and its MD is evaluated. When a Kalman filter fits to the assumptions differences between measurements and filter prediction, denoted as innovations, are realization of white noise. An purely data-driven approach based on this property is proposed by [Ber13] evaluating the whiteness of the innovations in a hypothesis test. The authors of [Gre+21a] propose an approach that is capable to reject the influence of the periodic input by estimating the periodic information with a non-steady state Kalman filter.
**Characteristics in the frequency domain** When models in the frequency domain are used the associated model characteristics such as the Frequency Response Functions (FRF) may be applied for damage detection. The authors of [FC03] propose an approach where a single FRF measured at several frequencies is used together with an analytic model of the reference structure to detect damage. In [SMS99] a FRF curvature method is developed that is an extension of the method in [PBS91] to all frequencies in the measurement range. Similarly, the authors of [Lim10] propose the interpolated damage detection method (IDDM) that uses the deviation of the deformed shape from smooth behavior expressed by the variation of an interpolation error. This method is shown to be appropriate when changes in environmental conditions are supposed enjoying that curvature characteristics are less effected by environmental conditions [Der+08].

The Power Spectral Density (PSD) is used in [WMF00], where the PSD peaks are evaluated with outlier analysis of their MD. The authors of [VSF14] show in a comparative study together with other time series methods that features based on PSD are quite sensitive. In [TA18] an approach is introduced that relies on the vibration energy in acceleration records using the normalized cumulative PSD. The method is able to account for deviations from perfect stationary white noise excitation.

Some research is made on the application of wavelet transforms, which decomposes a function into a set of wavelets helping to extract local spectral and temporal information simultaneously. [Ulr+16] introduce an approach for continuous wavelet transform that captures local mode shape irregularities by significantly magnified transform coefficients leading to a damage indicator inspired by MAC that expresses the similarity wavelet transforms in the reference and the testing state.

**Signal and its characteristics and associated signal processing functions** The simplest form is to use the signal directly without any compression procedure which leads to high feature dimension and difficult comparability. A method based on Moving Average Filter (MAF) is used in [KXY18] smoothing the signal with MAF to make deviations measurable. Many approaches that rely directly on the signal are based on machine learning concepts, to which some remarks are made in the following section.

An easier comparability is reached when signal statistics are used. In [FW12] some of them are evaluated in a numerical study presenting damage induced changes in the signal peak amplitude, signal mean and variance, and other signal statistics. It is pointed out that their sensitivity towards outliers limits the direct application.
With the goal to reduce the dimension of large data sets for easier interpretation PCA is applied providing a procedure with minimized loss of information. The data is projected to a lower dimensional space creating new uncorrelated variables, the principal components [JC16]. Using PCA e.g. in [Tan05] a feature extraction/de-noising method is developed that is based on frequency response signals. The authors of [Kul03] used PCA for dimensionality reduction of the modal parameters for the subsequent evaluation in control charts.

The application of functions from signal processing allows a different perspective on the data and is widely used. The authors of [FJ97] prove that the cross-correlation functions of the responses excited by ambient noise have similar characteristics as the structural impulse response functions and can be adopted for identifying structural dynamic properties. An approach based on the MD of autocorrelation coefficients is used in [HO14].

Another group of damage sensitive features is based on the data covariance. In [BAB00] an approach is proposed that is inspired by the covariance-driven SSI [BF85], which identifies the state-space representation of the observed system from the covariance Hankel matrix. For damage detection, a Gaussian-distributed residual is computed by confronting a covariance Hankel matrix in the testing state to the left null subspace of the covariance Hankel matrix in the reference state. Together with the statistical evaluation in a hypothesis test the method is referred to as stochastic subspace-based damage detection (SSDD). Based on the concept of SSDD further developments have been made to meet specific requirements for the application to engineering structures: the authors of [Bal+08a] derive an approach that considers changes in the environmental conditions. In [DM13b; DMH14] a residual robust to changes in the excitation covariance is derived and the authors of [DMZ16] enhance the method for damage localization and quantification. Other residual formulations on the basis of the covariance-driven SSI are derived in [YG06] considering that the relevant structural responses are comprised in the active subspace of the covariance Hankel matrix. An approach that uses the MD of the output covariance Hankel matrices can be found in [Gre+21b], where it is shown to be robust to changes in excitation. Comparison of the SSDD to other methods is made in [KF16] and [ZR09]. In the first, the sensitivity of the classical SSDD is compared to a method that considers modal parameters identified with SSI applying an automated mode identification algorithm to find stable poles. The latter compares the SSDD method to damage detection using modal parameters.
## 1.3 Feature evaluation and related uncertainties

When the damage sensitive feature is extracted and a comparable quantity is derived, the value of this quantity is evaluated with respect to its nominal value in an undamaged reference state. Since the feature usually is computed from noisy data, it is subject to statistical variability and so is the considered damage indicator (cf. Fig. 1.1). In the evaluation the statistical properties must be taken into account in order to define a threshold, from which it is decided if a change in the dynamical system is significant or results from natural variations due to noise. The threshold is based on the distribution properties of the damage indicator and an allowed false positive error probability. When the distribution is not known, it is a common procedure to determine it empirically from several reference data sets in the healthy state in a Monte Carlo study.

However, recently the proper consideration and quantification of the feature uncertainties have become an important part of research on vibration-based damage detection. When the statistical properties of the underlying feature are known, the properties of the resulting damage indicator can often be derived from them, e.g. with a first-order perturbation approach using the Delta Method [CB02]. In the following uncertainty quantification approaches for modal parameters, system matrices and related damage indicators are recalled and some groups of evaluation techniques are presented.

**Uncertainty quantification** First, the statistical distribution properties of a feature need to be derived for its uncertainty quantification. When the distribution is (asymptotically) Gaussian uncertainty quantification boils down to computation of the respective covariance. For example, computation procedures for uncertainty quantification of modal parameters estimated with SSI are in the focus in [RPD08; PGS07; MDM16] enhanced by efficient algorithms in [DM13a]. For modal parameters computed with an Eigensystem-Realization-Algorithm (ERA) the corresponding uncertainties are derived in [LM11] and when they are computed on the basis of Bayesian spectral density an approach can be found in [YK19]. In [Bro+18] the prior computation of the expected modal parameter uncertainty is used to optimize the measurement setups and their necessary duration. It is pointed out in [BB15] that the proper uncertainty quantification becomes particularly difficult when only few or short data sets are available.

Approaches computing the uncertainty of system matrices estimated with output-only and input-output subspace identification procedures and the resulting transfer function are derived in [DM13a; GDM20]. Formulations of the dominant parts of the estimation
errors of system matrices in a shift invariance and a state approach are presented in
[Ike14].

Uncertainty bounds related to the MAC estimates are derived in [GDM21], discrim-
inating between the MAC computed from the same mode shape and the one computed
from two different mode shape. Distribution properties of the Modal Phase Collinear-
ity (MPC), which is used to discriminate between real and complex-valued mode shape
vectors, are computed in [Gre+21c].

Uncertainty in data-driven residuals is treated e.g. in [DVG12] in the context of null space-based additive fault detection. Statistical properties of a Hankel matrix-based Mahalanobis distance for damage detection are investigated from the authors of [Gre+21b]. An approach for uncertainty computation of updated model parameters in the context of a FE model updating procedure is proposed in [Gau+17]. Other authors focus on explicit uncertainty quantification models of damage sensitive features such as Frequency Response Functions and modal parameters [Mao12; CM11].

**Control charts** In order to monitor the damage indicator with respect to its expected nominal value and nominal distribution, control charts developed for statistical process control by Walter Shewhart in the 1920s are used as graphical tools. The quantity of interest is evaluated in real-time by illustrating its current value in relation to the expected statistics. Applications can be found e.g. in [SCF00; YDG04; FSF01], where the mean of one-dimensional damage indicators are monitored by means of the X-Bar chart.

When the feature dimension is greater than one, multi-variate charts can be used, e.g. Hotelling $T^2$ or Multivariate Exponentially Weighted Moving Average (MEWMA) charts, which combine information from several variables in one scalar. The authors of [Kul03] evaluate the sensitivity of univariate versus multivariate control charts. They emphasize that in cases where it is not known in advance which parts of the damage sensitive feature are relevant it can be advantageous to use high-dimensional feature vectors and reduce the dimension subsequently by PCA, providing higher sensitivity of the control charts.

**Classical outlier or novelty analysis** Feature evaluation based on classical outlier or novelty analysis investigates the statistical properties of the damage sensitive feature. A widely used approach for multivariate quantities is the Mahalanobis Distance (MD), which is the distance between a point and a distribution accounting for the respective variance in the direction of the evaluated point by considering the covariance matrix. The
authors of [WMF00] introduce a systematic study on the application of outlier analysis
by MD pointing out that this approach assumes at least near-Gaussian distribution of the
undamaged feature and refers the reader to Kernel Density Estimation approaches if the
distribution is multi-modal.

In [Der+15] several outlier analysis concepts are studied comparing the MD to the Minimum Covariance Determinant (MCD) and the Maximum Volume Enclosing Ellipsoid (MVEE), emphasizing that the MCD is more robust to outliers in the training data.

The goal of the work in [BB15] is to investigate the MD with different versions of the covariance matrix and the authors of [Gre+17] compare a damage evaluation approach using the MD with statistical hypothesis testing in the context of SSDD.

The application of gap statistics as in [NK07] may belong to this group of evaluation
techniques. Usually applied to measure the variation of data within a cluster, in this
approach the difference between GMM is evaluated to derive the number of existing
mixtures.

**Statistical hypothesis tests** Statistical hypothesis testing is another evaluation tech-
nique that is widely used to classify a monitored system. In fact, it is tested whether
a given sample of data is drawn from a given probability distribution. Therefore, two
hypotheses are defined assuming different underlying distributions. The null hypothesis
holds when the tested data belongs to the unchanged reference system, the alternate hy-
pothesis applies when the data is recorded from a changed system. Various test statistics
exist, which take the covariance of the tested feature into account, and the choice depends
on the feature properties.

A hypothesis test for change detection in is derived in [BAB00] using the asymptotic
local approach to change detection from [BBM87]. The test is designed for the subspace-
based residual introduced in the previous section. Since the residual is asymptotically
Gaussian, general statistical framework, namely Generalized Likelihood Ratio (GLR) test,
can be applied and the resulting test statistic follows a χ²-distribution. The test is able to
take into account the variability of the evaluated feature due to noisy test data by means
of the feature covariance matrix. Assuming that a dynamical system can be described
by a set of system parameters, the damage detection can be done with respect to these
system parameters, without determining the parameter anew in the testing state. This
parametric version has the benefit that the theoretical distribution of the test statistic
can be described precisely in the undamaged and in the damaged state.
Other hypothesis tests are proposed e.g. in [SF01] where the standard deviation ratio of residual errors from an ARMA model is evaluated statistically in a Box-Andersen test; in [Ber13] where a Lag Shifted Whiteness Test (LSWT) is designed that evaluates the whiteness of the Kalman filter innovations; in [SM13] where Bayesian tests are used deriving the uncertainties in the detection, localization, and quantification step; in [Döh+15] where a statistical $\chi^2$-test was used in a damage localization procedure; in [DVG12] where data-driven residuals are evaluated in a $\chi^2$-test providing an analytic design of the damage detection procedure considering measurement noise as well as identification uncertainty. Some applications to engineering structures can be found e.g. in [Döh+14; GK17; Ven+14].

For generalization of hypothesis testing for multiclass classification problems the authors of [FW12] provide suggestions.

**Model updating** When model updating procedures [FM95] are applied, the metric formulation and the evaluation step in Fig. 1.2 often are combined. Such approaches are inverse methods [Fri07] deriving updated FE models from data measured in the testing state. A function representing the distance between the measured data or data feature and the predicted data from the FE model is minimized by means of optimization methods. Many model updating methods are developed in the past decades. In [Wae+16] three techniques are compared that are based on the classical Tikhonov regularization, Constitutive Relation Error and a Bayesian approach using physical parameters such as the Young modulus as updating parameter. The Bayesian model updating is a suitable method handling uncertainties as emphasized by the authors of [SDL15]. They give an overview how to consider the uncertainties of the measurements and those related to the numerical model properly in the damage detection procedure on the basis of Bayesian and fuzzy approaches. Others provide approaches for uncertainty considerations, which use automated fault estimation based on interval arithmetic [KWA19] or which are based on the covariance estimation procedure for the updated model parameters [Gau+17]. More detailed overviews can be found in [Bro+01] or [FW12].

**Supervised learning** An increasing number of approaches come from the vibrant field of supervised learning methods, such as artificial neural networks, where the steps of the damage detection as displayed in Fig. 1.1 usually are hidden within the end-to-end algorithms. The stochastic character of the data is considered within the algorithms. These methods process the raw data of the time history itself or even plots of them
[Tan+19], as well as previously derived features requiring some classification information
of the already processed features [Mig+12]. A more detailed discussion of some concepts
for vibration-based damage detection can be found in [FW12] and [Avc+21].

1.4 Data-based damage detection under changing en-
vironmental conditions

The application of automated vibration-based damage diagnosis methods often fails due to changing environmental conditions. This is e.g. observed by the authors of [Far+94; PD00; Roh+00a] applying such SHM methods over long time periods to bridge structures. The dynamical behavior of the structures is significantly affected by environmental effects such as temperature, changes in the excitation statistics, humidity, mass loading, or soil conditions [Soh07; Ube+17; WSF02].

Since the impact of these effects to the dynamics is quite high compared to small
damages [PMD01; LD07; Far+00], automated methods are perturbed in two ways. Either
the data that is considered for the baseline contains several environmental nominal states.
In this case high variability of the damage indicator in the normal state is assumed, which
masks the effect of damages leading to poor sensitivity of the test. In the other case, when
the testing data is recorded at a normal state with fixed environmental conditions, changes
in the structural dynamics due to new conditions during the testing phase lead to high
false alarm rates.

Several approaches exist that aim at robust damage detection procedures by data
normalization in the sense of separating signal changes caused by environmental varia-
tions from structural changes [Soh07]. The normalization step can be performed at several
stages of the damage detection procedure, i.e. during the feature extraction by choosing
environmental insensitive features or prepare them adequately, during the computation
of the damage indicator, or during the evaluation step. In the following, two groups of
concepts are distinguished: The first group is based on the assumption that the environ-
mental conditions are measured and the other group works without further knowledge
about the detailed environmental conditions.
1.4.1 Methods with relation to measured environmental param-
eters

In this section methods are recalled, which assume that the environmental conditions
can be characterized uniquely by measurable external parameters and the related changes
of the feature or damage indicator distribution can be expressed as a function of this
parameter. Concepts that are based on regression or interpolation are widely used. Also
approaches exist that consider classification of the testing situation or parameter depen
dent sensitivities.

**Regression and Interpolation** In [PD00] a study on one-year data from the Z24 bridge in Switzerland is presented showing a strong dependency of the frequencies from temperature while other environmental conditions have no significant impact. The temperature dependency is not linear over time, but has a declination gradient change at a temperature of 0°C. To model the temperature dependency a black-box autoregressive model with exogenous terms (ARX) is fitted to the frequency-temperature relation obtained from data using temperatures as input and delivering the frequencies. The prediction error is used as damage indicator and compared to an expected distribution property. Similarly in [Soh+98] a linear adaptive model is used to reproduce the natural variability of the frequencies of a bridge leading to a baseline for discrimination of deviations in data sets at new testing temperatures. In [Hu+15] regression functions are used to remove the temperature effects on the frequencies and the residual is evaluated by outlier analysis. The authors of [Oli+18] propose a SHM procedure for application to wind turbines with removing environmental effects with regression from modal parameters. An approach considering a SSI for modal identification is presented in [KF16], deriving regression models for the modal parameters and their distribution properties. The authors of [MCC12] apply static and dynamic regression models complemented by a PCA. In [Ram+10] masonry structures are evaluated regarding the modal parameters and MAC using regression analysis to reject changes in environmental conditions.

Approaches in [WSF02] use AR coefficients as damage sensitive features, where dif-
ferent dimensions of the AR model are considered. A regression analysis is performed on
the damage indicator – the MD of the coefficients – on each component of its statistics. It
turned out that with this approach the eigenvalues of covariance matrix become negative
in some cases, hence an interpolation approach is also derived. In this interpolation the
statistics of the damage indicator is computed from the neighboring statistics applying
linear interpolation. The approach provides some capability for generalization for multi-parameter situations. The results are compared to a previous approach when large training data sets spanning all environmental conditions of interest were used, showing that the regression method is much more sensitive.

**Classification** Some approaches rely on the idea to classify the current testing feature to one particular reference feature of classification. In [FMG03] an approach is derived within the setup of a smart structure concept. It is based on the SSDD that uses a reference null space to set up a comparable residual. Several references are stored in a data base and for a particular testing temperature the closest reference is used. The method is developed further in [FKB13] computing a reference model based on weighted references.

An approach based on Bayesian system identification and model class selection is proposed in [Sim+20], where a set of competing parameterized probabilistic structural models are postulated, which explicitly account for the effect of varying ambient temperatures on the mechanical properties of the system. Information provided by the structural models and by recorded vibration data leads to the parameters of the temperature-dependent structural models. From this the plausible state of the system is obtained.

**Parameter dependent sensitivities** Statistical rejection approaches using parameter dependent sensitivities to compute residuals that are tested subsequently in hypothesis tests are derived e.g. in [Bal+09]. From data sets at some reference temperatures a combined system and temperature effect model is obtained and changes due to temperature are rejected as nuisance by using corresponding sensitivities. A temperature adjusted reference that is based on temperature related modal parameters computed from temperature models is derived in [Bas+10].

In [Bhu+19] an approach based on load vectors is proposed, where a sensitivity-based correction of the identified modal parameters in the damaged state is performed with respect to the temperature field and the temperature dependent parameters of the finite element model in the reference state.

### 1.4.2 Methods with unmeasured environmental conditions

In cases where the environmental conditions are unmeasured or cannot be described uniquely robust features that can be formulated insensitive to changes in the environmental conditions or large normal condition models are applied.
**Robust feature or damage indicator** It is found that e.g. mode shapes are less sensitive to changes in the environmental conditions than frequencies are [Der+08]. This motivates for approaches that are based on features or damage indicators which can be formulated insensitive to such changes.

In the first stage, features can be employed that are less affected by environmental conditions. This includes e.g. the IDDM in [Lim10] that is based on FRF and uses the deviation of the deformed shape from smooth behavior, and is shown to be adequate for damage detection under changes in the environmental conditions.

In other approaches the damage sensitive features are processed to be insensitive to changes in the environmental parameters. E.g. PCA is used to remove the temperature affected parts of a feature. In [Yan+05a] PCA is used for linear cases, where the environmental parameters are taken into account as embedded variables. The same authors [Yan+05b] propose a two-step procedure for non-linear cases – clustering of the data space into several regions and application of PCA to each local region – which allows performing a piecewise linearization. In another publication [KBD13] PCA is used to project the data into a subspace where the variance from the environment is small. Nonlinear PCA with autoassociative neural networks (AANN) is proposed in [Kra91] generalizing the mapping into feature space to allow arbitrary nonlinear functionalities for multivariate data analysis. The authors of [Der+08] propose a feature that is based on modal filtering tuned to a particular mode observed in frequency domain output. In the damage detection the fact is used that damage produces spurious peaks in the frequency domain output of the modal filter while environmental changes shift the peak in the frequency domain without changing the general shape. The authors of [Van+05] introduce damage detection based on robust singular value decomposition (SVD), which distinguishes between changes in the working conditions and damage by evaluating changes in the rank of a feature vector matrix. In [Cro+12] also a strategy is proposed that prepares the data for robust feature extraction by using cointegration, which is based on finding a stationary linear combination of non-stationary time series to purge the data set of its common trends leading to a robust residual.

In the context of robustness towards changes in the excitation the authors of [DM13b] and [Gre+21b] propose robust residuals which are based on output covariance Hankel matrices either using insensitive subspaces of the Hankel matrix or normalizing them before the damage indicator is computed.

The intention in [Man02] is to derive univariate damage indicators. The method is
based on MD, where the approach is applied to lamb-wave data. Similarly the authors of
[Der+15] introduce a robust outlier analysis that uses the MCD, which is able to remove
the masking effect and search inclusively for multiple outliers. An approach based on
GMM is explained in [Kul14]. Non-linear effects from environmental variations can be
removed by considering Minimum Mean Square Error (MMSE) and applying PCA.

**Large normal condition models** Another concept when no information is available on the current environmental state is based on large normal condition models. One vibrant field is the application of machine learning approaches to environmentally affected engineering structures, such as in [GGW17; Tan+19]. In [SW10] a negative selection algorithm is used for cases of damage detection where the data under normal conditions is significantly non-Gaussian or environmental conditions lead to variability. A combination of time series analysis, neural networks, and statistical inference approaches is developed in [SWF02]. The effect of damage on the extracted feature is separated from environmental effects by means of an AR model and AR-ARX model.

Other methods derive global models that represent the effect of the environmental pa-
rameters to the statistical properties of damage sensitive features. The authors of [RWD14]
propose an approach that identifies damage sensitive features in a reference state from
each part of a training data broken up into time sequences that are short compared to the
time of parameter variations. Using kernel PCA a data-driven global nonlinear system
model is derived from the identified features. The features in the testing state are then
compared to the feature predicted from the global model. In [SCS16] polynomial chaos
expansion (PCE) is used to find a model that represents the extracted features in relation
to the measured operational conditions.

Approaches that are based on average models are developed e.g. in [Bal+08a]. Here
a global reference model is built from a merged data set smoothing out the temperature
effect by the averaging operation. In fact, this can be understood as a kind of interpola-
tion with no respect to the testing temperature. Similarly, the authors of [HF09] provide
an approach that uses Constant Coefficient Pooled Vector Autogressive with Exogenous
Variables (CCP-VARX) representations, which provide averaged descriptions of the struc-
tural dynamics over temperature, while the functionally pooled FP-VARX counterparts
consider the explicit modeling of temperature dependence of the dynamics.

These approaches provide some robustness in practice, since they do not require ex-
plicit measurements of the environmental conditions. However, if the reference measure-
ments are not representative the approaches may fail, since no further relation between
the current measurement in the testing state and the large normal model is taken into
account.
<figure><img src="image_9.png" /></figure>
# BACKGROUND THEORY

## 2.1 Introduction

The developments in this thesis are based on the stochastic subspace-based damage detection (SSDD) and the associated residual vector. The residual uses subspace properties of an output covariance Hankel matrix. It is defined as the product of the left null space of the Hankel matrix in the reference state and the Hankel matrix estimated from current measurement data in the test state [BAB00; DMH14]. The statistical evaluation of the residual in a hypothesis test yields a damage detection test statistic that considers the uncertainties of the residual.

In the context of this thesis the theory on complete consideration of the residual uncertainties is developed, using perturbation theory, and a SSDD approach robust to changes in the environmental conditions is derived from local model interpolation. This chapter recalls the used theories and methods.

After introducing the models that describe the dynamics of a mechanical structure in section 2.2, section 2.3 gives an overview on how those models and its modal parameters can be extracted from vibration data in time domain by means of the stochastic system identification (SSI). Then, the current state of the SSDD method is presented in section 2.4. The mathematical perturbation theory including the delta method and the basics of a local model interpolation approach for linear parameter varying (LPV) systems are summarized in sections 2.5 and 2.6.
## 2.2 Structural dynamics of output-only systems

The dynamics of a linear time-invariant (LTI) mechanical structure can be described in continuous time $t$ by the following equation of motion

$$M\ddot{q}(t) + C\dot{q}(t) + Kq(t) = \nu(t). \quad (2.1)$$

The vector $q \in \mathbb{R}^m$ collects the displacements at the $m$ degrees of freedom (DOF) of the system and denote the derivatives with respect to time. $M$, $C$, and $K \in \mathbb{R}^{m \times m}$ are the mass, damping, and stiffness matrices of the system, respectively. $\nu$ is the excitation function that is represented by a white noise term in output-only systems with unknown excitation.

With the system state

$$x(t) = \begin{bmatrix} q(t) \\ \dot{q}(t) \end{bmatrix} \in \mathbb{R}^n, \quad (2.2)$$

where $n = 2m$ is the model order, a first-order system in state space can be derived from (2.1) [Jua94]. Sampled at rate $1/\tau$ such that $t = k\tau$, the stochastic state-space model of an output-only system at time step $k$ writes as

$$\begin{cases} x_{k+1} = Ax_k + w_k \\ y_k = Cx_k + v_k. \end{cases} \quad (2.3)$$

$y_k \in \mathbb{R}^r$ is the output of the system measured at $r$ sensor positions with measurement noise $v_k \in \mathbb{R}^r$. The state noise term $w_k \in \mathbb{R}^n$ is related to the unknown excitation. Both noise terms are unmeasured, and usually assumed to be stationary, centered, white, and having finite fourth moment. $A$ and $C$ are the state transition matrix and the observation matrix of the discrete time system. They are directly linked to the mechanical system matrices from (2.1) with

$$A = \exp \left( \begin{bmatrix} 0 & I \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix} \tau \right) \in \mathbb{R}^{n \times n}$$ $$C = [L_d - L_a \mathbf{M}^{-1}\mathbf{K} \quad L_v - L_a \mathbf{M}^{-1}\mathbf{C}] \in \mathbb{R}^{r \times n}. \quad (2.4)$$

$L_d$, $L_v$, and $L_a \in \{0, 1\}^{r \times m}$ are selection matrices denoting the sensor locations at the
corresponding DOF in $q$, where displacements $d$, velocities $v$, and accelerations $a$ are recorded.

## 2.3 System identification

Several approaches exist that are used to identify the systems in (2.1) and (2.3) from data. The determination of the particular system matrices in (2.1) from data only is impossible and is usually treated in the context of FE model updating [FM95]. Nevertheless, the dynamic properties of the system can be retrieved by the identification of its eigenstructure, namely the eigenfrequencies, damping ratios, and mode shapes. The state-space matrices in (2.3) can be identified directly from data, e.g. with subspace identification methods. From the identified matrices the eigenvalues, the eigenvectors, and the corresponding modal parameters can be derived in a subsequent step.

It has to be stated that the matrices and parameters, which are identified from data, are estimates and usually afflicted with uncertainties. This is due to the stochastic character of the measured data caused by the noise terms in (2.3). The estimates of the true variable are denoted by $\hat{\cdot}$ in the following.

In this thesis a SSI algorithm is used to identify modal parameters from data. The concept of the SSI and in particular the covariance-driven subspace identification [BF85] are recalled in the following and it is shown how to derive the eigenstructure of the dynamical systems. The theoretical background gives some context for the later introduced SSDD method.

### 2.3.1 Concept of stochastic subspace identification

The following concept of the SSI is inherent in many algorithms, as e.g. in [VD96; PD99; BM07]. It is based on a matrix $H_{p+1,q}$ computed from data projections and/or data covariances, whose column space is the observability matrix $O$ of the system arising in the context of control theory [Kál63]. The observability matrix writes as

$$O_{p+1} = \begin{bmatrix} C \\ CA \\ \vdots \\ CA^p \end{bmatrix} \in \mathbb{R}^{(p+1)r \times n}, \qquad (2.5)$$
where $p$ is a size parameter denoting the block rows, and $A$ and $C$ are the state-space matrices introduced in section 2.2. For system identification the two matrices have to be derived from the observability matrix. $C$ follows as the first block-row of (2.5). Due to the shift invariance property of $O_{p+1}$, which states that $A$ is the over the time steps consistent term in $O_{p+1}$, the state transition matrix can be approximated with a least square approach from

$$O_{p+1}^{\uparrow} A = O_{p+1}^{\downarrow}, \qquad (2.6)$$

where

$$O_{p+1}^{\uparrow} = \begin{bmatrix} C \\ CA \\ \vdots \\ CA^{p-1} \end{bmatrix}, \quad O_{p+1}^{\downarrow} = \begin{bmatrix} CA \\ CA^2 \\ \vdots \\ CA^p \end{bmatrix} \in \mathbb{R}^{pr \times n}. \qquad (2.7)$$

The observability matrix is derived in the SSI as follows. The matrix $\mathcal{H}_{p+1,q}$ holds the factorization property

$$\mathcal{H}_{p+1,q} = W O_{p+1} Z_q \in \mathbb{R}^{(p+1)r \times qr}. \qquad (2.8)$$

Herein $W$ is an invertible weighting function, often defined to be the identity matrix. The definition of $Z_q$ depends on the chosen SSI algorithm and the corresponding type of $\mathcal{H}_{p+1,q}$. $q$ denotes the block columns, with usually $q = p + 1$, and it is $\min\{pr, qr\} > n$. The observability matrix can be computed from the image subspace of the matrix $\mathcal{H}_{p+1,q}$ with

$$O_{p+1} = W^{-1} U_1 D_1^{1/2}, \qquad (2.9)$$

Image $U_1$ and the singular values in $D_1$ are computed by a singular value decomposition (SVD) of matrix $\mathcal{H}_{p+1,q}$, truncated at the assumed model order $n$, leading to

$$\mathcal{H}_{p+1,q} = U D V^T = \begin{bmatrix} U_1 & U_2 \end{bmatrix} \begin{bmatrix} D_1 & 0 \\ 0 & D_2 \end{bmatrix} \begin{bmatrix} V_1^T \\ V_2^T \end{bmatrix}. \qquad (2.10)$$

$D_1 \in \mathbb{R}^{n \times n}$ denotes the diagonal matrix of $n$ singular values and $D_2 \in \mathbb{R}^{(p+1)r-n \times qr-n}$ is assumed to be zero, or containing noise in the case of the matrix estimate. $U_1 \in \mathbb{R}^{(p+1)r \times n}$ and $U_2 \in \mathbb{R}^{(p+1)r \times (p+1)r-n}$ are the column space and the left null space, respectively. $V_1 \in \mathbb{R}^{qr \times n}$ spans the row space and $V_2 \in \mathbb{R}^{qr \times qr-n}$ denotes the null space of matrix $\mathcal{H}_{p+1,q}$.

A consistent estimate of the observability matrix is obtained, when $\mathcal{H}_{p+1,q}$ is replaced
by an estimate $\hat{\mathcal{H}}_{p+1,q}$.

## 2.3.2 Covariance-driven subspace identification

A well-known SSI algorithm is the covariance-driven subspace identification [BF85]. It uses the output covariances to built matrix $\mathcal{H}_{p+1,q}$ and is applied in this work.

Let $\mathcal{Y}_N = \{y_1, y_2, ..., y_N\}$ be a set of $N$ output data samples at $r$ sensor locations. The cross-covariances between the states and the output are $G = \mathbf{E}(x_{k+1}y_k^T)$ and the output covariances are computed to [VD96]

$$R_i = \mathbf{E}(y_{k+i} y_k^T) = CA^{i-1}G. \qquad (2.11)$$

Matrix $\mathcal{H}_{p+1,q}$ is the block Hankel matrix of output covariances with

$$\mathcal{H}_{p+1,q} \stackrel{\text{def}}{=} \begin{bmatrix} R_1 & R_2 & \dots & R_q \\ R_2 & R_3 & \dots & R_{q+1} \\ \vdots & \vdots & \ddots & \vdots \\ R_{p+1} & R_{p+2} & \dots & R_{p+q} \end{bmatrix} \stackrel{\text{def}}{=} \text{Hank}(R_i). \qquad (2.12)$$

The factorization property in (2.8) and $W = I$ yield

$$\mathcal{H}_{p+1,q} = \mathcal{O}_{p+1,q}C_q, \qquad (2.13)$$

since it becomes $Z_q = [G \ A G \ \dots \ A^{q-1}G]$, which is exactly the well-known controllability matrix $C_q$.

A consistent estimate of the block Hankel matrix $\hat{\mathcal{H}}_{p+1,q} \stackrel{\text{def}}{=} \text{Hank}(\hat{R}_i)$ can be obtained from data by using the empirical output covariances

$$\hat{R}_i = \frac{1}{N} \sum_{k=i+1}^{N} y_{k+i} y_k^T. \qquad (2.14)$$

## 2.3.3 Identification of the eigenstructure

With the previous SSI algorithm the state transition and the observation matrix can be derived on the basis of (2.5). From them the eigenstructure of system (2.3) follows with the eigenvalue problem of $A$

$$A\phi_l = \lambda_l\phi_l, \qquad (2.15)$$
where $\lambda_l$ and $\phi_l$ denote the l-th eigenvalue and eigenvector, respectively.

The l-th eigenfrequency $f_l$ and damping ratio $\xi_l$ of the system in continuous time in (2.1) are computed from $\lambda_l$ of the discrete time system as follows. The relation between $\lambda_l$ and the eigenvalue of the continuous-time system $\mu_l$ is defined by

$$ \lambda_l = e^{\mu_l \tau}. \qquad (2.16) $$

From this and regarding $\omega_l = 2\pi f_l$ and the damping ratio $\xi_l$ it holds

$$ \mu_l = -\omega_l \xi_l \pm \omega_l \sqrt{1 - \xi_l^2} i, \qquad (2.17) $$

leading to

$$ f_l = \frac{|\mu_l|}{2\pi}, \quad \xi_l = \frac{-\mathbf{Re}(\mu_l)}{|\mu_l|} \qquad (2.18) $$

The observable eigenvectors or mode shapes $\varphi_l$ can be derived with the observation matrix

$$ \Psi = C\Phi, \qquad (2.19) $$

with $\Psi = [\varphi_1 \varphi_2 \dots \varphi_n]$ and $\Phi = [\phi_1 \phi_2 \dots \phi_n]$ being matrices that contain the mode shapes and the eigenvectors in its columns, respectively.

Consistent estimates of the modal parameters can be derived by using estimates of $\hat{A}$ and $\hat{C}$, which result from an observability matrix that is computed on a Hankel matrix estimate in (2.10).

## 2.4 Stochastic subspace-based damage detection

The stochastic subspace-based damage detection (SSDD) method as presented in [BAB00] makes use of the previous described covariance-driven SSI. With a data-driven, asymptotically Gaussian residual changes in a parameter vector are uncovered. The parameter vector is a collection of system parameters, which describe the dynamical behavior of a system completely. Changes in the parameters may be caused by damages [FW07] and can be detected by evaluating the residual statistically. The evaluation is done by means of the asymptotic local approach to change detection. In the following, the theoretical background is recalled, introducing the residual function and its statistical properties first, followed by a summary of the statistical damage testing.
2.4.1 Subspace-based residual function

It is assumed that the dynamical properties of the system in (2.1) or (2.3) can be completely described by some system parameters, collected in a parameter vector $\theta$. Common parametrizations of mechanical structures are physical parameters such as mass, damping, and stiffness parameters, or the eigenstructure. The nominal value of the parameter vector in the undamaged state is denoted by $\theta^*$ and detecting changes in $\theta$ becomes the monitoring task. In the following notation it is assumed that theta contains the modal parameters and is defined as

$$ \theta = \begin{bmatrix} \Lambda \\ \mathrm{vec}(\Phi) \end{bmatrix}, \qquad (2.20) $$

where $\Lambda = [\lambda_1, \lambda_2, ..., \lambda_n]^T$ is a vector containing all eigenvalues.

In accordance to the covariance-driven SSI in section 2.3.2 a block Hankel matrix $\mathcal{H}_{p+1,q}^*$ (2.12) is considered that is built from the output covariances $R_i$ in an undamaged reference state. $\mathcal{H}_{p+1,q}^*$ can be factorized into the observability and the controllability matrix (2.13), which contain the state-space matrices $A$ and $C$ (2.5) and implicitly parameter vector $\theta$. Consequently, a residual function on the basis of $\mathcal{O}_{p+1}$ allows for monitoring the parameter vector $\theta$. The observability matrix in the modal basis can be built from the system parameters in $\theta^*$ with

$$ \mathcal{O}_{p+1}(\theta^*) = \begin{bmatrix} \Phi \\ \Phi\Delta \\ \vdots \\ \Phi\Delta^{p-1} \end{bmatrix} \in \mathbb{C}^{(p+1)r \times n}, \qquad (2.21) $$

with the diagonal matrix $\Delta \stackrel{\text{def}}{=} \mathrm{diag}(\Lambda)$.

The residual function is based on the left null space $S(\theta^*)$ of $\mathcal{O}_{p+1}(\theta^*)$ that is obtained from a SVD analogous to (2.10) with $S(\theta^*) = U_2$. It holds

$$ S(\theta^*)^T \mathcal{O}_{p+1}(\theta^*) = 0. \qquad (2.22) $$

Due to factorization property (2.13), $\mathcal{H}_{p+1,q}^*$ and $\mathcal{O}_{p+1}(\theta^*)$ share the same left null space. Thus, with a Hankel matrix that is computed under $\theta = \theta^*$ it also holds

$$ S(\theta^*)^T \mathcal{H}_{p+1,q}^* = 0. \qquad (2.23) $$
In the residual function a covariance Hankel matrix estimated from a data set $\mathcal{Y}_N$ of length $N$ of the current monitoring state is confronted to the reference matrix $S(\theta^*)$. The residual function is defined by

$$ \zeta(\theta^*, \mathcal{Y}_N) \stackrel{\text{def}}{=} \sqrt{N} \operatorname{vec}(S(\theta^*)^T \hat{\mathcal{H}}_{p+1,q}), \qquad (2.24) $$

with

$$ \mathbf{E}_{\theta}(\zeta(\theta^*, \mathcal{Y}_N)) = 0 \quad \text{iff} \quad \theta = \theta^*. \qquad (2.25) $$

$\mathbf{E}_{\theta}$ denotes the expectation when $\mathcal{Y}_N$ is recorded from a system with system parameter $\theta$. Hence, the residual mean is zero when the system is in the reference state, since the testing Hankel matrix has the same left null space that was determined in the reference state. The mean value deviates from zero when the system is damaged.

Although the idea for the residual relies on (2.22), usually the reference matrix $S(\theta^*)$ is computed directly as the left null space of the covariance Hankel matrix in the reference state. This avoids the system identification step to derive $\theta^*$. In this case, condition (2.25) still applies, since $\mathcal{H}_{p+1,q}$ and $\mathcal{O}_{p+1}$ share the same left null space, as remarked in the context of (2.23).

## 2.4.2 Statistical damage testing

The statistical evaluation of the residual in (2.24) requires its distribution properties, which are unknown in general. To solve this problem, the asymptotic local approach to change detection [BBM87] is used. This approach is based on the concept of local asymptotic normality [Le 56] and compares the statistical distributions under close hypotheses.

The close hypotheses are formulated as

$$ \begin{aligned} \mathbf{H}_0: & \theta = \theta^* && (\text{reference system}) \\ \mathbf{H}_1: & \theta = \theta^* + \delta/\sqrt{N} && (\text{damaged system}), \end{aligned} \qquad (2.26) $$

and discriminate between the parameter vector in the undamaged and the damaged state. The parameter change $\delta$ is unknown but fixed and consequently very small damages can be detected by increasing the number of data samples $N$.

With the close hypotheses, the central limit theorem applies for the residual function, i.e. the residual is asymptotically Gaussian distributed as shown in [BAB00]. For $N \to \infty$
it holds

$$ \zeta(\theta^*, \mathcal{Y}_N) \xrightarrow{d} \begin{cases} \mathcal{N}(0, \Sigma_\zeta(\theta^*)) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\mathcal{J}_\zeta(\theta^*) \delta, \Sigma_\zeta(\theta^*)) & \text{under } \mathbf{H}_1, \end{cases} \qquad (2.27) $$

where “d” denotes convergence in distribution. $\mathcal{J}_\zeta(\theta^*)$ is a Jacobian matrix and denotes the asymptotic sensitivity of the residual with respect to the system parameters in the reference state. For a system that is parametrized with respect to its modal parameters, the asymptotic residual sensitivity was derived in [BAB00; BMG04]. It is written as

$$ \mathcal{J}_\zeta(\theta^*) = (\mathcal{O}_{p+1}(\theta^*)^\dagger \mathcal{H}_{p+1,q} \otimes S(\theta^*))^T \mathcal{O}'_{p+1}(\theta^*). \qquad (2.28) $$

The computation of the derivative $\mathcal{O}'_{p+1}$ is recalled in section 2.5. $\Sigma_\zeta(\theta^*)$ is the residual covariance. When it can be assumed that the noise properties in the dynamical system do not change, there is no need to recompute the covariance in the testing state [BAB00]. A detailed computation procedure of the residual covariance $\Sigma_\zeta(\theta^*)$ is presented in the subsequent section.

With the residual being asymptotically Gaussian distributed (2.27) the damage detection task is to decide if the mean of a Gaussian variable is significantly different from zero. Assuming $\mathcal{J}_\zeta(\theta^*)$ to be of full rank, the residual can be tested in a generalized likelihood ratio test [BAB00]. The test writes as

$$ t = \zeta(\theta^*, \mathcal{Y}_N)^T \Sigma_\zeta(\theta^*)^{-1} \mathcal{J}_\zeta(\theta^*) \left( \mathcal{J}_\zeta(\theta^*)^T \Sigma_\zeta(\theta^*)^{-1} \mathcal{J}_\zeta(\theta^*) \right)^{-1} \mathcal{J}_\zeta(\theta^*)^T \Sigma_\zeta(\theta^*)^{-1} \zeta(\theta^*, \mathcal{Y}_N). \qquad (2.29) $$

The Gaussian distribution properties of the residual yield an asymptotically $\chi^2$-distributed test statistic $t$. The Gaussian residual with zero mean in the reference state leads to a central $\chi^2$-distribution. From the mean value $\mathcal{J}_\zeta\delta$ in the damaged case a non-central $\chi^2$-distribution follows and it holds

$$ t \xrightarrow{d} \begin{cases} \chi^2(\nu, 0) & \text{under } \mathbf{H}_0 \\ \chi^2(\nu, \lambda) & \text{under } \mathbf{H}_1. \end{cases} \qquad (2.30) $$

For a full rank $\mathcal{J}_\zeta(\theta^*)$, the degree of freedom $\nu$ and the non-centrality parameter $\lambda$ in the damaged case are computed with

$$ \nu = \operatorname{rank}(\mathcal{J}_\zeta(\theta^*)) , \quad \lambda = \delta^T (\mathcal{J}_\zeta^T(\theta^*) \Sigma_\zeta^{-1}(\theta^*) \mathcal{J}_\zeta(\theta^*)) \delta. \qquad (2.31) $$

For cases, where the sensitivity matrix is not of full rank, the reader is referred to [Men20].
A special case are the non-parametric tests [Bal+08b], which are applied without focusing on a particular system parametrization $\theta$. This corresponds to a direct detection of changes in the residual and $J = I$ [DMH14; MDV21]. They are attractive due to their simplicity, since they avoid a detailed system identification even in the reference state. The central limit theorem from (2.27) becomes

$$ \zeta(\theta^*, \mathcal{Y}_N) \rightarrow \begin{cases} \mathcal{N}(0, \Sigma_\zeta(\theta^*)) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\gamma, \Sigma_\zeta(\theta^*)) & \text{under } \mathbf{H}_1, \end{cases} \quad (2.32) $$

where $\gamma$ is a change in the residual mean due to damage. The test writes as

$$ t = \zeta(\theta^*, \mathcal{Y}_N)^T \Sigma_\zeta(\theta^*)^{-1} \zeta(\theta^*, \mathcal{Y}_N), \quad (2.33) $$

and $\nu$ and $\lambda$ become

$$ \nu = \mathrm{rank}(\Sigma_\zeta(\theta^*)) , \quad \lambda = \gamma^T \Sigma_\zeta(\theta^*)^{-1} \gamma. \quad (2.34) $$

### 2.4.3 Asymptotic residual covariance

For the statistical evaluation of the residual in (2.29) its covariance is required. This covariance cannot be calculated directly, but has to be derived from the covariance of the underlying Hankel matrix $\Sigma_H$. The computation algorithm for $\Sigma_H$ depends on the computation method for $H$. For the covariance-driven SSI an efficient computation algorithm of the Hankel covariance is given in [DM13a], taking problems due to small number of data sets into account. This algorithm is recalled in the following.

Let $\tilde{\mathcal{H}}_{p+1,q}$ be estimated on a data set of length $N$. The asymptotic covariance of the vectorized Hankel matrix $\Sigma_H$ is defined as

$$ \Sigma_H \stackrel{\mathrm{def}}{=} \lim_{N \to \infty} \mathrm{cov}(\sqrt{N} \mathrm{vec} \tilde{\mathcal{H}}_{p+1,q}). \quad (2.35) $$

To estimate the covariance from one reference data set, the data is cut into $n_b$ blocks, each of length $N_b$ such that $N = N_b \cdot n_b$. From each data block of length $N_b$ a Hankel matrix $\tilde{\mathcal{H}}_{p+1,q}^{(s)}$, with $s = 1, \dots, n_b$ is computed. With $\tilde{\mathcal{H}}_{p+1,q} = \frac{1}{n_b} \sum_{s=1}^{n_b} \tilde{\mathcal{H}}_{p+1,q}^{(s)}$ the well-known computation procedure of the sample covariance leads to

$$ \hat{\Sigma}_H = \frac{N_b}{n_b - 1} \sum_{s=1}^{n_b} \mathrm{vec}(\tilde{\mathcal{H}}_{p+1,q}^{(s)} - \tilde{\mathcal{H}}_{p+1,q}) \mathrm{vec}(\tilde{\mathcal{H}}_{p+1,q}^{(s)} - \tilde{\mathcal{H}}_{p+1,q})^T. \quad (2.36) $$
The asymptotic residual covariance is derived from $\hat{\Sigma}_H$ using the statistical delta method to propagate the uncertainties of the underlying Hankel matrix to the residual. The residual covariance estimate writes as [DMH14]

$$ \hat{\Sigma}_{\zeta}(\theta^*) = (I \otimes S^T(\theta^*)) \hat{\Sigma}_H (I \otimes S(\theta^*)), \quad (2.37) $$

In case of the non-parametric test (2.33) the parameter $\theta^*$ is omitted.

### 2.4.4 Decision making from test results

To draw conclusions from the statistical tests in (2.29) or (2.33) the distribution of $t$ in the undamaged state is considered in order to define a threshold for decision making.

In theory $t$ follows a $\chi^2$-distribution defined by two parameters, $\nu$ denoting the DOF and the non-centrality parameter $\lambda$ (which is 0 in the undamaged state). The number of DOF of the distribution is affected by the system parametrization considered in parameter vector $\theta$. From (2.31) it follows that for a system parametrized with respect to its modal parameters, it holds $\nu = 2mr$. When e.g. a number of $\bar{m}$ structural parameters are used for a full parametrization of the system, it becomes $\nu = \bar{m}$. The non-centrality parameter $\lambda$ in (2.31) depends on the damage (changed parameter and extent of change), as it is affected by $\delta$ only [MDV21]. In the undamaged state $\delta$ is zero and the distribution of the test is the central $\chi^2$-distribution. However, a theoretical value of $\lambda$ cannot be derived, since the residual covariance cannot be obtained from theory but is estimated from data.

In applications the statistical properties usually are evaluated empirically on the basis of a Monte Carlo study. From the empirical distribution of the test values a threshold is derived such that false alarms occur with a desired probability, e.g. with less than 1% as illustrated in Fig.2.1 (left). A test value in the testing phase above the threshold will classify the structure to be damaged.

Fig. 2.1 (right) shows the probability density function (pdf) for a $\chi^2$-distributed random variable for two different DOF in the central case and for $\lambda = 50$. For a small number of DOF the distribution is narrow. The spread increases with larger DOF resulting in a more significant overlap of the central and the non-central distribution. In the context of damage detection based on the test statistic $t$ this means that the smaller the dimension of the parametrization, the smaller damages can be detected [MDV21].

Very common characteristics for the evaluation of the damage detection test perfor-
<figure><img src="image_10.png" /><figcaption>Figure 2.1 – χ²-distribution for different parameters (left) and definition of threshold (right).</figcaption></figure>

mance are the probability of detection (POD) as a measure of the test sensitivity, and the probability of false alarms (PFA), which usually is part of the conceptualization on the damage detection task and leads to the choice of the damage detection threshold. They are defined as

$$POD = \frac{\text{true positives detected}}{\text{total true positive}}, \quad PFA = \frac{\text{false positive detected}}{\text{total true negative}}. \qquad (2.38)$$

The PFA depends on the spread of the test values in the reference state and the chosen threshold, but not on the damage size. The POD is affected by the spread of the test values, the threshold, and the damage size, where big damages lead to high POD. Consequently, a high threshold avoids false alarms, but depending on the damage size many damages will be undetected. A low threshold leads to a higher sensitivity, but also results in a higher PFA.

To evaluate the test performance independently from the choice of threshold, the relation between POD and PFA can be considered in a receiver operating characteristic (ROC) curve, e.g. described in [SCN09] in the context of noisy non-destructive testing devices. For this curve, couples of POD and PFA are computed for different thresholds at a fixed damage size. A typical ROC curve is demonstrated in Fig. 2.2. Hence, the best performance of a test is when no false alarms occur and the POD is 100%. The more the curve approaches this point, the better it performs. When the curve approaches 100% POD for any PFA, this can be considered as an indication of decoupled POD and PFA. The line of no performance denotes the cases, where POD equals the PFA.
<figure><img src="image_11.png" /><figcaption>Figure 2.2 – Typical receiver operating characteristic curve.</figcaption></figure>

## 2.5 Considerations for uncertainty computation

In the previously described damage detection methodology the decision over the structure's state to be damaged or undamaged relies on a statistical test (2.29). In the test the uncertainties related to the residual are considered in terms of the residual covariance $\Sigma_{\zeta}(\theta^*)$. The covariance is of particular importance for the test behavior. In general, the covariance can be estimated by the sample covariance. However, the direct computation of such a sample covariance is often inconvenient, and instead a sensitivity-based propagation of the (sample) covariance of an underlying random vector is performed. In this thesis such sensitivities are derived by means of perturbation theory. In this section the concept of perturbation propagation is introduced referring to the delta method. Some relevant uncertainty computations from literature are recalled subsequently.

### 2.5.1 Mathematical remarks

In the following, some mathematical considerations are summarized that are used in the subsequent developments.

Let $A$, $B$, $C$, and $D$ be matrices, such that the matrix multiplication in the following equations applies. For the Kronecker product "$\otimes$" it holds

$$AC \otimes BD = (A \otimes B)(C \otimes D). \qquad (2.39)$$

With "vec" denoting the column stacking operator it is

$$\mathrm{vec}(ABC) = (C^T \otimes A) \mathrm{vec}(B). \qquad (2.40)$$
**Permutation matrix** $\mathcal{P}$ A vectorized matrix $X \in \mathbb{R}^{a \times b}$ can be linked to the vectorization of its transposed by a permutation matrix $\mathcal{P}_{a,b}$ such that

$$ \mathrm{vec}(X^T) = \mathcal{P}_{a,b} \mathrm{vec}(X). \qquad (2.41) $$

The permutation matrix $\mathcal{P}_{a,b}$ is defined as

$$ \mathcal{P}_{a,b} \stackrel{\mathrm{def}}{=} \sum_{k=1}^{a} \sum_{l=1}^{b} E_{k,l}^{a,b} \otimes E_{l,k}^{b,a} = [I_a \otimes e_1 \ I_a \otimes e_2 \ \dots \ I_a \otimes e_b], \qquad (2.42) $$

with properties

$$ \mathcal{P}_{a,b}^{T} \mathcal{P}_{a,b} = \mathcal{P}_{a,b} \mathcal{P}_{a,b}^{T} = I_{ab} $$

$$ \mathcal{P}_{a,b}^{T} = \mathcal{P}_{b,a}. $$

$E_{k,l}^{a,b} \in \mathbb{R}^{a \times b}$ are matrices which are equal to 1 at entry $(k, l)$ and zero elsewhere, $I_a$ is the identity matrix of size $a \times a$, and $e_i \in \mathbb{R}^b$ are unit vectors with the entry 1 at position $i$ and zero elsewhere. With matrix $H \in \mathbb{R}^{c \times d}$ it holds

$$ (X \otimes H) \mathcal{P}_{b,d} = \mathcal{P}_{a,c}(H \otimes X). \qquad (2.43) $$

## 2.5.2 Concept of uncertainty propagation

The perturbation of a function of a random vector is computed from the perturbation of its underlying random vector by means of the delta method. Let $\hat{X}$ be a random vector whose true value is $X$ and with covariance $\Sigma(\hat{X})$. The first-order Taylor approximation of vector valued functions $g(\hat{X})$ yields

$$ g(\hat{X}) \approx g(X) + \Delta g(X) = g(X) + \mathcal{J}_{g,X}(\hat{X} - X), \qquad (2.44) $$

where the derivative or sensitivity $\mathcal{J}_{g,X}$ is the Jacobian matrix in the multivariate case. The Jacobian writes as

$$ \mathcal{J}_{g,X} = \frac{\partial g}{\partial X} = \begin{bmatrix} \frac{\partial g_1}{\partial x_1} & \cdots & \frac{\partial g_1}{\partial x_k} \\ \vdots & \ddots & \vdots \\ \frac{\partial g_j}{\partial x_1} & \cdots & \frac{\partial g_j}{\partial x_k} \end{bmatrix} \qquad (2.45) $$
for function $g : \mathbb{R}^k \to \mathbb{R}^j$, and is evaluated at $X$. $\Delta g(X)$ is the perturbation of the function and thus the covariance of the vector valued function follows from $\mathrm{vec}\Delta g(\mathrm{vec}\Delta g)^T$. It can be approximated with

$$ \Sigma(g(\hat{X})) \approx \mathcal{J}_{g,X} \Sigma(\hat{X}) \mathcal{J}_{g,X}^T. \qquad (2.46) $$

With the delta method [CB02] it is shown that the approximated covariance of a function of a random variable is correct at the limit, i.e. when the number of samples $N$ goes to infinity. For a sequence of random vectors $\hat{X}_N$ converging to a true value $X$ for $N \to \infty$ and satisfying the central limit theorem

$$ \sqrt{N}(\hat{X}_N - X) \xrightarrow{d} \mathcal{N}(0, \Sigma) \qquad (2.47) $$

it holds

$$ \sqrt{N}(g(\hat{X}_N) - g(X)) \xrightarrow{d} (0, \mathcal{J}_{g,X} \Sigma \mathcal{J}_{g,X}^T). \qquad (2.48) $$

In most cases $g(X)$ is unknown and the Jacobian (2.45) has to be derived by means of mathematical perturbation theory using first-order perturbations denoted by $\Delta$. This is convenient for asymptotic Gaussian variables, e.g. [MDM16]. The basic idea is that $\mathcal{J}_{g,X}$ can be obtained from the relation between $\Delta g(X)$ and $\Delta X$, due to $\Delta g(X) \approx \mathcal{J}_{g,X} \Delta X$ (2.44), where $\Delta X = \hat{X} - X$ for $\hat{X}$ close to $X$. In the application case, a perturbation $\Delta$ is introduced to the vector valued function. The $\Delta$-terms are isolated, leading to the needed sensitivity by the remaining part of the perturbed function.

### 2.5.3 Relevant uncertainty computations from literature

This section recalls some uncertainty computation procedures from literature, which are used later in the contribution part of the thesis.

**Residual sensitivity $\mathcal{J}_{\zeta,H}$ related to the testing Hankel matrix** The concept of uncertainty propagation is firstly applied to a simple example deriving the sensitivity of the residual as it is used in (2.37). The distribution properties of the conventional residual are derived from the distribution of the underlying Hankel estimate $\hat{\mathcal{H}}_{p+1,q}$. The index denoting the block row and block column sizes is dropped in the following for simplicity. The distribution of the Hankel matrix is not exactly known but is shown to be
asymptotically Gaussian [Han70] and writes as

$$ \sqrt{N}\mathrm{vec}(\hat{\mathcal{H}} - \mathcal{H}) \xrightarrow{d} \mathcal{N}(0, \Sigma_{\mathcal{H}}). \qquad (2.49) $$

With the delta method for a function $f(\hat{\mathcal{H}}) = \mathrm{vec}(S^T \hat{\mathcal{H}})$ it follows

$$ \sqrt{N}\mathrm{vec}(f(\hat{\mathcal{H}}) - f(\mathcal{H})) \xrightarrow{d} \mathcal{N}(0, \mathcal{J}_f \Sigma_{\mathcal{H}} \mathcal{J}_f^T). \qquad (2.50) $$

Due to the definition of $S$ as the left null space of $\mathcal{H}$ and assuming it is deterministic, the function $f$ is zero at the true value of $\mathcal{H}$. (2.50) becomes

$$ \sqrt{N}\mathrm{vec}(S^T \hat{\mathcal{H}}) \xrightarrow{d} \mathcal{N}(0, \mathcal{J}_f \Sigma_{\mathcal{H}} \mathcal{J}_f^T), \qquad (2.51) $$

which is exactly the residual (2.24) and its distribution in the undamaged state. Derivative $\mathcal{J}_f$ is derived by introducing a perturbation into function $f$ and with (2.40) it yields

$$ \Delta f(\mathcal{H}) = \mathrm{vec}(S^T \Delta \mathcal{H}) = (I \otimes S^T)\mathrm{vec}(\Delta \mathcal{H}) = \mathcal{J}_f \mathrm{vec}(\Delta \mathcal{H}). \qquad (2.52) $$

### Perturbation propagation to modal parameters from SSI
An algorithm for perturbation propagation from the estimated Hankel matrix to the modal parameters with the SSI procedure is given in [PGS07; RPD08]. A numerical efficient computation is derived in [DM13a]. The computation of the sensitivities is recalled in the following.

The propagation follows the computation procedure of the modal parameters from the Hankel matrix that was described in section 2.3. With a SVD of the covariance Hankel matrix $\mathcal{H}$ the observability matrix $\mathcal{O}$ is derived (2.9) leading to the state-space matrices $A$ and $C$ (2.5). Subsequently, the modal parameters $\lambda$ and $\varphi$ in $\theta^*$ are computed from $A$ and $C$ (section 2.3.3). The same procedure is considered for the propagation of the perturbations from the Hankel matrix to the identified modal parameters, leading to the associated sensitivity $\mathcal{J}_{\theta^*, \mathcal{H}} = \mathcal{J}_{\theta^*, A|C} \mathcal{J}_{A|C, \mathcal{O}} \mathcal{J}_{\mathcal{O}, \mathcal{H}}$.

Let $I_a$ denote a $a \times a$ identity matrix and $0_{a,b}$ be the $a \times b$ zero matrix. $\mathcal{P}$ is the permutation matrix in (2.42). The sensitivity $\mathcal{J}_{\mathcal{O}, \mathcal{H}} \in \mathbb{R}^{(p+1)^{rn \times (p+1)rqr}}$ of the observability matrix writes as

$$ \mathcal{J}_{\mathcal{O}, \mathcal{H}} \stackrel{\mathrm{def}}{=} \frac{1}{2} (I_n \otimes U_1 \Sigma_1^{-1/2}) F_1 \begin{bmatrix} (v_1 \otimes u_1)^T \\ \vdots \\ (v_n \otimes u_n)^T \end{bmatrix} + (\Sigma_1^{-1/2} \otimes [I_{(p+1)r} \quad 0_{(p+1)r,qr}]) \begin{bmatrix} B_1^\dagger C_1 \\ \vdots \\ B_n^\dagger C_n \end{bmatrix}, \qquad (2.53) $$
where

$$B_i \stackrel{\text{def}}{=} \begin{bmatrix} I_{(p+1)r} & -\frac{1}{\sigma_i} \mathcal{H}_{p+1,q} \\ -\frac{1}{\sigma_i} \mathcal{H}^{0\ T} & I_{qr} \end{bmatrix}, \quad C_i \stackrel{\text{def}}{=} \frac{1}{\sigma_i} \begin{bmatrix} v_i^T \otimes (I_{(p+1)r} - u_i u_i^T) \\ (u_i^T \otimes (I_{qr} - v_i v_i^T)) \mathcal{P}_{(p+1)r,qr} \end{bmatrix}. \quad (2.54)$$

$\sigma_i$, $u_i$, and $v_i$ are the i-th singular values, left and right singular vectors of the Hankel matrix. The selection matrix $S_1$ is

$$F_1 \stackrel{\text{def}}{=} \sum_{k=1}^{n} E_{(k-1)n+k,k}^{n^2,n} \quad (2.55)$$

The sensitivities of system matrices A and C with respect to the observability matrix $\mathcal{J}_{A,O} \in \mathbb{R}^{n^2 \times (p+1)rn}$ and $\mathcal{J}_{C,O} \in \mathbb{R}^{rn \times (p+1)rn}$ are derived to

$$\mathcal{J}_{A,O} \stackrel{\text{def}}{=} (I_n \otimes \mathcal{O}^{\uparrow\dagger} F_2) - (A^T \otimes \mathcal{O}^{\uparrow\dagger} F_3) + ((\mathcal{O}^{\downarrow\uparrow} F_2 - A^T \mathcal{O}^{\uparrow\dagger} F_2) \otimes (\mathcal{O}^{\uparrow\dagger} \mathcal{O}^{\uparrow})^{-1}) \mathcal{P}_{(p+1)r,n}, \quad (2.56)$$

with $\mathcal{O} = \mathcal{O}_{p+1}$, the selection matrices

$$F_2 \stackrel{\text{def}}{=} [0_{pr,r} \quad I_{pr}], \quad F_3 \stackrel{\text{def}}{=} [I_{pr} \quad 0_{pr,r}], \quad (2.57)$$

and

$$\mathcal{J}_{C,O} \stackrel{\text{def}}{=} I_n \otimes [I_r \quad 0_{r,pr}]. \quad (2.58)$$

For the sensitivities of the eigenvalues and mode shapes $\mathcal{J}_{\lambda_i,A} \in \mathbb{C}^{1 \times n^2}$ and $\mathcal{J}_{\varphi_i,A/C} \in \mathbb{C}^{r \times (n^2+rn)}$ we have

$$\mathcal{J}_{\lambda_i,A} \stackrel{\text{def}}{=} \frac{1}{\chi_i^* \phi_i} (\phi_i^T \otimes \chi_i^*) \quad (2.59)$$

$$\mathcal{J}_{\varphi_i,A/C} \stackrel{\text{def}}{=} \frac{1}{(C\phi_i)_k} (I_r - [0_{r,k-1} \quad \phi_i \quad 0_{r,r-k}]) [C \mathcal{J}_{\phi_i,A} \quad \phi_i^T \otimes I_r], \quad (2.60)$$

where $\lambda_i$, $\phi_i$, $\chi_i$ and $\varphi_i$ are the i-th eigenvalue, right eigenvector, left eigenvector of A and the i-th mode shape respectively. * denotes the complex conjugate transpose, and $\mathcal{J}_{\phi_i,A}$ is defined as

$$\mathcal{J}_{\phi_i,A} \stackrel{\text{def}}{=} (\lambda_i I_n - A)^{\dagger} \left( \phi_i^T \otimes (I_n - \frac{\phi_i \chi_i^*}{\chi_i^* \phi_i}) \right). \quad (2.61)$$
**Perturbation propagation to parametric observability matrix** In [DMH14] the sensitivity of the parametric observability matrix in the modal basis $O_{p+1}(\theta^*)$ with respect to the modal parameters in $\theta^*$ is derived. It is computed as the derivative of the real-valued parametric observability matrix $O'_{p+1}(\theta^*)$.

The system parameter $\theta^*$ is separated into its conjugate complex pairs $\theta_c^*$ and $\theta_c^*$. Corresponding to this $\phi_c$ and $\Delta_c$ are computed such that $\phi = [\phi_c \ \phi_c]$ and $\Lambda = [\Lambda_c^T \ \Lambda_c^T]$ and $\Delta_c = \text{diag}(\Lambda_c)$.

The real-valued parametric observability matrix writes as

$$O_{p+1}(\theta^*) \stackrel{\text{def}}{=} \begin{bmatrix} \Re(\phi_c) & \Im(\phi_c) \\ \Re(\phi_c \Delta_c) & \Im(\phi_c \Delta_c) \\ \vdots & \vdots \\ \Re(\phi_c \Delta_c^p) & \Im(\phi_c \Delta_c^p) \end{bmatrix}, \qquad (2.62)$$

leading to the derivative

$$\mathcal{J}_{\mathcal{O}, \theta^*} \stackrel{\text{def}}{=} \begin{bmatrix} \Re(\mathcal{O}'_{p+1}(\theta_c)) & -\Im(\mathcal{O}'_{p+1}(\theta_c)) \\ \Im(\mathcal{O}'_{p+1}(\theta_c)) & \Re(\mathcal{O}'_{p+1}(\theta_c)) \end{bmatrix}. \qquad (2.63)$$

$\mathcal{O}'_{p+1}(\theta_c)$ is the complex-valued derivative of the vectorized observability matrix

$$\mathcal{O}'_{p+1}(\theta_c) \stackrel{\text{def}}{=} \begin{bmatrix} \Lambda_1^{(p)} \otimes \varphi_1 & 0 & \Lambda_1^{(p)} \otimes I_r & 0 \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \Lambda_m^{(p)} \otimes \varphi_m & 0 & \Lambda_m^{(p)} \otimes I_r \end{bmatrix}, \qquad (2.64)$$

where

$$\Lambda_i^{(p)} \stackrel{\text{def}}{=} [1 \ \lambda_i \ \lambda_i^2 \ \dots \ \lambda_i^p]^T, \quad \Lambda_i^{(p)} \stackrel{\text{def}}{=} [0 \ 1 \ 2\lambda_i \ \dots \ p\lambda_i^{p-1}]^T \qquad (2.65)$$

for $1 \le i \le m$.

## 2.6 Model interpolation of linear parameter varying systems

In this thesis a damage detection method is developed that is robust to environmental changes. The method relies on the stochastic subspace-based damage detection and uses
a direct model interpolation method from [Zha18] for linear parameter varying (LPV) systems that is recalled in this section.

## 2.6.1 Concept

LPV systems are systems whose dynamical behavior depends on an over the time varying scheduling parameter $P$. They can be approximated by LTI systems at fixed parameter working points. Assume a LPV system from which output data is available at several working points $P_j$ with $j = 1, 2, ..., u$. A model of the system at a particular parameter point $P$, where $P$ is in general different from the previous $P_j$, can be obtained by the direct model interpolation from [Zha18]. In this approach local state-space models are interpolated to a large global model, which is subsequently reduced and from which data $\mathcal{Y}(P)$ at parameter $P$ can be characterized. This data is shown to be equivalent to interpolated data from coherent local models.

## 2.6.2 Model interpolation

The LTI system at the $j$-th working point is represented by a local version of (2.3) with

$$ egin{cases} x_{j,k+1} = A_j x_{j,k} + w_{j,k} \\ y_{j,k} = C_j x_{j,k} + v_{j,k}. \end{cases} \tag{2.66} $$

Assume data sets in the undamaged state are available at parameter working points $P_j$ with $j = 1, ..., u$. For any parameter $P$ a global model can be derived as

$$ egin{aligned} \underline{x}_{k+1} &= \underline{A} \underline{x}_k + \underline{w}_k \\ y_k &= \underline{C}(P) \underline{x}_k + \underline{v}_k(P), \end{aligned} \tag{2.67} $$

at time step $k$ and with

$$ \underline{A} = \begin{bmatrix} A_1 & & & \\ & \ddots & & \\ & & \ddots & \\ & & & A_u \end{bmatrix}, \quad \underline{x}_k = \begin{bmatrix} x_{1,k} \\ \vdots \\ x_{u,k} \end{bmatrix}, \quad \underline{w}_k = \begin{bmatrix} w_{1,k} \\ \vdots \\ w_{u,k} \end{bmatrix}, $$

$$ \underline{C}(P) = [\rho_1(P)C_1 \cdots \rho_l(P)C_u], $$

$$ \underline{v}_k(P) = \sum_{j=1}^{u} \rho_j(P) v_{j,k}(t). $$
The state transition matrix and the state vector are not interpolated but are concatenated from all working points $j$.

Weighting functions $\rho(P)$ are chosen such that

$$ \sum_{j=1}^{u} \rho_j(P) = 1 \quad \text{for all } P \in \mathbb{P} \quad \text{and} \quad \rho_j(P) : \mathbb{P} \Rightarrow [0, 1] \qquad (2.68) $$

centered at $P$.

The global model (2.67) is of order $mn$ and subsequently reduced by means of the well-known balanced reduction method. In the linearly transformed state vector $\underline{x}_k$ the smallest Hankel singular values are removed. This leads to a state space model of order $n$, same as for the local models.

## 2.7 Summary

In this chapter the theories and methods that are used in the context of this thesis were recalled. In particular, the background of the stochastic subspace-based damage detection was described in detail, since this method forms the basis of the presented developments.
# PART II

## Contribution to theory
In this part the methodological contribution of the thesis to the research field of automated vibration-based damage detection for civil structures is presented. The developments are based on the covariance-driven stochastic subspace based damage detection (SSDD). The first chapter is related to distribution properties of the SSDD residual when the reference model is estimated from data and thus afflicted with uncertainties. The second chapter provides a SSDD method that is robust do changes in the environment. Each chapter contains a theory section together with a proof of concept, as well as a discussion section.
# STOCHASTIC SUBSPACE-BASED DAMAGE DETECTION WITH UNCERTAIN REFERENCE NULL SPACE

## 3.1 Introduction

The stochastic subspace-based damage detection (SSDD) method (section 2.4) first presented in [BAB00] works on a residual function, in which a reference null space is confronted to a data covariance Hankel estimate. The residual is evaluated statistically by means of hypothesis testing in a generalized likelihood ratio test, where the residual's covariance is an important part of the test. The covariance is derived from the uncertainty of the Hankel matrix, while the reference null space is assumed to be perfectly known without any impact on the residual uncertainty.

The assumption of a deterministic reference null space is true in ideal cases. However, in real world applications, the reference $S(\theta^*)$ usually is estimated from data. Hence, the residual is a function of two random variables, namely the Hankel matrix estimate in the testing state and the estimate of the reference null space. This has a significant impact on the residual covariance for the damage detection test (2.29), (2.33). Neglecting the uncertainties of the reference null space herein leads to test statistics that deviate from the theoretical values of the expected $\chi^2$-distribution.

This is not only a theoretical problem, but also affects the application in practice when it comes to the determination of a threshold for decision making. Often empirical thresholds are computed from data in the reference state by Monte Carlo studies and the discrepancy between the computed and the theoretical test values is discarded. It is assumed that enough data in the reference state is available. However, especially in cases when monitoring systems are applied retroactively due to indication of initial damage usually only few data of the reference system is available. Thresholds that are derived
from theory can confirm results from empirical studies with few realizations and thus improve the reliability of the test.

In this chapter the gap of the SSDD method for real-world applications with estimated reference is closed. A test that considers the uncertainties related to the reference matrix and the testing data is derived. The residual is formulated as a function of the two random variables $\hat{S}(\theta^*)$ and $\hat{H}$ to which the central limit theorem applies. The statistical delta method is used to derive the statistical properties of the residual and its asymptotic covariance. The associated hypothesis test is presented allowing for theoretical based threshold and predictable test behavior for damage detection.

In section 3.2 the theory for the two random variable dependent residual is derived, and the corresponding hypothesis test is presented. Subsequently, in section 3.3, the new test is evaluated on numerical data. The statistical behavior of the new test is verified towards theoretical expectations, and the test performance is discussed in detail. The impact of considering the uncertainty related to the reference null space is emphasized, and effects of data lengths, system parametrization, or reference null space computation procedure are analyzed.

## 3.2 Damage detection with uncertain reference

In this section the theory for a residual with an uncertain reference is derived. The developments are based on the SSDD explained in section 2.4. After putting the problem into a formalized context for SSDD, the distribution properties of a residual with uncertain reference are derived. An algorithm for computing the residual covariance is presented using perturbation theory introduced in section 2.5.

### 3.2.1 Formalized problem statement

The subspace-based residual of the conventional SSDD in section 2.4 is defined as a function of one random variable (cf. (2.24)) with

$$ \zeta(\theta^*, \mathcal{Y}_N) = \sqrt{N} \mathrm{vec}(S(\theta^*)^T \hat{H}_{p+1,q}). \quad (3.1) $$

$S(\theta^*)$ denotes the reference matrix of a system with system parameter $\theta^*$ and is assumed to be perfectly known. $\hat{H}_{p+1,q}$ is the covariance Hankel matrix estimated from a data set of length $N$ in the monitoring state. $p$ and $q$ are parameters for the block row and block
column size of the Hankel matrix and will be dropped in the following for simplicity. The statistical properties of the residual are described in (2.27) with zero mean in the reference state due to property $S(\theta^*)^T\mathcal{H} = 0$. The covariance of the residual is denoted by $\Sigma_\zeta(\theta^*)$ (2.37) and is derived from the uncertainties of the Hankel matrix estimate in the monitoring state.

Now assume that in real-world applications the reference matrix is computed in the reference state from a Hankel matrix $\hat{\mathcal{H}}^*$ that is estimated from a reference data set $\mathcal{Y}_M$ of length $M$. Then the reference matrix itself is a random variable afflicted with uncertainties. The fully data-driven residual must be defined as a two random variable function with

$$ \tilde{\zeta}(\theta^*, \mathcal{Y}_N) = \sqrt{N} \mathrm{vec}(\hat{S}(\theta^*)^T \hat{\mathcal{H}}), \quad (3.2) $$

and it follows that $\hat{S}(\theta^*)^T\mathcal{H}$ only approximates zero. Consequently, the distribution properties in (2.27) do not hold for a residual with an estimated reference null space.

The distribution properties of the conventional residual are derived from the distribution of the underlying Hankel estimate $\hat{\mathcal{H}}$ that is not exactly known. It is shown to be asymptotically Gaussian [Han70] and under the close hypotheses (2.26) it writes as

$$ \sqrt{N}\mathrm{vec}(\hat{\mathcal{H}} - \mathcal{H}^*) \xrightarrow{d} \begin{cases} \mathcal{N}(0, \Sigma_{\mathcal{H}}(\theta^*)) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\mathcal{J}_{\mathcal{H}}(\theta^*) \delta, \Sigma_{\mathcal{H}}(\theta^*)) & \text{under } \mathbf{H}_1 \end{cases} \quad (3.3) $$

$\delta$ is the change of the system parameters, $\mathcal{J}_{\mathcal{H}}(\theta^*)$ is the sensitivity of $\mathrm{vec}(\mathcal{H})$ with respect to $\theta$ evaluated at $\theta^*$, and $\Sigma_{\mathcal{H}}(\theta^*)$ is the Hankel matrix covariance that is the same as in the reference state under the assumption of unchanged noise properties [BBM87; DMH14] and follows from (2.35).

Now remark that the central limit theorem recalled from (2.27) with

$$ \zeta(\theta^*, \mathcal{Y}_N) \xrightarrow{d} \begin{cases} \mathcal{N}(0, \Sigma_\zeta(\theta^*)) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\mathcal{J}_\zeta(\theta^*) \delta, \Sigma_\zeta(\theta^*)) & \text{under } \mathbf{H}_1, \end{cases} \quad (3.4) $$

is a consequence of (3.3), since multiplication of (3.3) with deterministic matrix $S(\theta^*)^T$ yields $\sqrt{N}\mathrm{vec}(S(\theta^*)^T(\hat{\mathcal{H}}-\mathcal{H}^*))$ on the left. This is equal to residual $\zeta$ due to $S(\theta^*)^T\mathcal{H}^* = 0$, which still is an asymptotically Gaussian vector (with different mean and covariance). This allows to state the residual distribution (2.27) without any knowledge on the theoretical limit value of $\mathcal{H}^*$ that is unknown in practice.

However, when replacing $S(\theta^*)$ by its estimate $\hat{S}(\theta^*)$, this does not hold anymore
since $\hat{S}(\theta^*)^T \mathcal{H}^* \neq 0$, and the central limit theorem in (2.27) would become incorrect. Since $\hat{S}(\theta^*)$ is a consistent estimate of $S(\theta^*)$, the asymptotic mean of the residual is not affected, but the fact that $\hat{S}(\theta^*)^T \mathcal{H}^*$ is only approximately but not exactly zero leads to a modification of the residual covariance. This computational inaccuracy has often been discarded in previous works [Döh+14], and empirical thresholds had to be used for the hypothesis test. For a correct computation of the corresponding damage detection test the distribution properties of $\tilde{\zeta}$ need to be derived including the exact expression for its covariance.

### 3.2.2 Asymptotic distribution of the residual

The asymptotic distribution of the residual is derived from the joint distribution of the testing Hankel matrix $\mathcal{H}$ and the reference null space Hankel matrix $\mathcal{H}^*$, considering the effect of the data lengths $M$ and $N$. $N$ is the length of the test data set in the monitoring state to which the residual is normalized. The distribution properties of the conventional residual in (2.27) are derived with respect to this length from the Hankel matrix distribution (3.3). However, the reference null space estimate $\hat{S}(\theta^*)$ is computed on a data set of length $M$. Consequently, the error that is induced in the residual by the uncertainty afflicted reference null space depends on this $M$. The distribution of the Hankel matrix in the reference state $\hat{\mathcal{H}}^*$ is asymptotically Gaussian with

$$ \sqrt{M} \mathrm{vec}(\hat{\mathcal{H}}^* - \mathcal{H}^*) \xrightarrow{d} \mathcal{N}(0, \Sigma_{\mathcal{H}}(\theta^*)) \quad (3.5) $$

for $M \to \infty$, similar as in (3.3). $\Sigma_{\mathcal{H}}(\theta^*)$ is the asymptotic covariance of the data Hankel matrix (2.35) and is the same in the reference and the damaged state under the assumption of unchanged noise properties [BBM87; DMH14]. When the data-driven residual is analyzed as a two random variable function depending on $\hat{S}(\theta^*)$ and $\hat{\mathcal{H}}$ both random variables contribute to the distribution properties of the residual and the joint asymptotic distribution of them is needed.

Denote the two random variable function with

$$ g(\hat{h}) = \mathrm{vec}(\hat{S}(\theta^*)^T \hat{\mathcal{H}}) \quad (3.6) $$

where

$$ \hat{h} = \begin{bmatrix} \mathrm{vec}(\hat{\mathcal{H}}^*) \\ \mathrm{vec}(\hat{\mathcal{H}}) \end{bmatrix}, \quad \text{and} \quad h = \begin{bmatrix} \mathrm{vec}(\mathcal{H}^*) \\ \mathrm{vec}(\mathcal{H}^*) \end{bmatrix}, $$
and $g(h) = 0$. With the statistical delta method [CB02] from section 2.5 the asymptotically Gaussian distribution properties of $\hat{h}$ are propagated to $g(\hat{h})$, i.e. to the residual written as $\tilde{\zeta} = \sqrt{N}(g(\hat{h}) - g(h))$. Since $\hat{\mathcal{H}}^*$ is computed on a data set of length $M$ but the residual is normalized with data length $N$ of the test state, the ratio $c = \frac{N}{M}$ is introduced. Multiplication of $\sqrt{c}$ to the asymptotic distribution properties of $\hat{\mathcal{H}}^*$ in (3.3) leads to

$$ \sqrt{M} \sqrt{\frac{N}{M}} \mathrm{vec}(\hat{\mathcal{H}}^* - \mathcal{H}^*) \xrightarrow{d} \mathcal{N}(0, c\Sigma_{\mathcal{H}}(\theta^*)). \quad (3.7) $$

It can be stated that $\hat{\mathcal{H}}$ and $\hat{\mathcal{H}}^*$ are independent, as they are computed from different data sets. Consequently, the cross-covariances between the two Hankel matrices will be zero. With the distribution properties of the Hankel matrix under the close hypotheses (3.3) and from (3.7) the joint distribution yields for $\sqrt{N}(\hat{h} - h)$ as

$$ \sqrt{N} \left( \begin{bmatrix} \mathrm{vec}(\hat{\mathcal{H}}^*) \\ \mathrm{vec}(\hat{\mathcal{H}}) \end{bmatrix} - \begin{bmatrix} \mathrm{vec}(\mathcal{H}^*) \\ \mathrm{vec}(\mathcal{H}^*) \end{bmatrix} \right) \xrightarrow{d} \begin{cases} \mathcal{N}\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} c\Sigma_{\mathcal{H}}(\theta^*) & 0 \\ 0 & \Sigma_{\mathcal{H}}(\theta^*) \end{bmatrix}\right) & \text{under } \mathbf{H}_0 \\ \mathcal{N}\left(\begin{bmatrix} 0 \\ \mathcal{J}_{\mathcal{H}}(\theta^*)\delta \end{bmatrix}, \begin{bmatrix} c\Sigma_{\mathcal{H}}(\theta^*) & 0 \\ 0 & \Sigma_{\mathcal{H}}(\theta^*) \end{bmatrix}\right) & \text{under } \mathbf{H}_1 \end{cases} \quad (3.8) $$

From the joint distribution and the first order Taylor expansion of $g(\hat{h})$

$$ g(\hat{h}) = g(h) + \mathcal{J}_g(\hat{h} - h), \quad \mathcal{J}_g = \frac{\partial g(h)}{\partial h}, \quad (3.9) $$

it can be concluded according to (2.48) based on the delta method that $\tilde{\zeta} = \sqrt{N}(g(\hat{h}) - g(h))$ is asymptotically Gaussian with

$$ \tilde{\zeta}(\theta^*, \mathcal{Y}_N) \xrightarrow{d} \begin{cases} \mathcal{N}(0, \tilde{\Sigma}_\zeta(\theta^*)) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\tilde{\mathcal{J}}_\zeta(\theta^*)\delta, \tilde{\Sigma}_\zeta(\theta^*)) & \text{under } \mathbf{H}_1. \end{cases} \quad (3.10) $$

$\delta$ and $\tilde{\mathcal{J}}_\zeta(\theta^*) = \mathcal{J}_\zeta(\theta^*)$ have not changed with respect to (2.27). The close hypotheses definition (2.26) yields that the asymptotic residual covariance is the same under both hypotheses writing as

$$ \tilde{\Sigma}_\zeta(\theta^*) = \mathcal{J}_g \begin{bmatrix} c \Sigma_{\mathcal{H}}(\theta^*) & 0 \\ 0 & \Sigma_{\mathcal{H}}(\theta^*) \end{bmatrix} \mathcal{J}_g^T. \quad (3.11) $$

The sensitivity $\mathcal{J}_g$ is obviously of form $\mathcal{J}_g = [\mathcal{J}_1 \ \mathcal{J}_2]$, where $\mathcal{J}_1$ is related to the uncer-
tainties of the reference null space and $J_2$ considers the perturbation propagated from the Hankel matrix in the testing state. The sensitivity is derived in the following section, leading to the asymptotic residual covariance.

In the non-parametric case the residual distribution is analogous to (2.32)

$$ \tilde{\zeta}(\theta^*, \mathcal{Y}_N) \xrightarrow{d} \begin{cases} \mathcal{N}(0, \tilde{\Sigma}_\zeta(\theta^*)) & \text{under } \mathbf{H}_0 \\ \mathcal{N}(\tilde{\gamma}, \tilde{\Sigma}_\zeta(\theta^*)) & \text{under } \mathbf{H}_1, \end{cases} \quad (3.12) $$

where $\tilde{\gamma}$ is a change in the residual mean due to damage and has not changed with respect to (2.32).

### 3.2.3 Asymptotic covariance of the residual

The asymptotic residual covariance $\tilde{\Sigma}_\zeta(\theta^*)$ follows from the Hankel matrix covariance (2.36) with sensitivity $J_g$. Perturbation theory (section 2.5) is used to propagate the uncertainties from the Hankel matrices $H^*$ and $H$ to the residual.

Considering the definition of $g(\hat{h})$ (3.6), the sensitivity $J_g$ can be computed as the derivative of $\text{vec}(S(\theta^*)^T H)$ with respect to $\text{vec}(H^*)$, since $S(\theta^*)$ depends on $H^*$, and to $\text{vec}(H)$. It is obtained through first-order perturbation of the asymptotically Gaussian variables. From the first-order Taylor approximation (3.9) the perturbation in $g$ follows as $\Delta g = J_g \Delta h$. With (3.6) this writes as

$$ \text{vec}(\Delta(S(\theta^*)^T H)) = J_g \begin{bmatrix} \text{vec}(\Delta H^*) \\ \text{vec}(\Delta H) \end{bmatrix}. \quad (3.13) $$

In (3.13) the perturbation of the residual (left side) is related to the reference matrix $S(\theta^*)$ and the tested Hankel matrix $H$. Being a sum of both perturbations it writes as

$$ \text{vec}(\Delta(S(\theta^*)^T H)) = \text{vec}(\Delta S(\theta^*)^T H) + \text{vec}(S(\theta^*)^T \Delta H). \quad (3.14) $$

With matrix operations in (2.40) and (2.41) and perturbation matrix $P$ from (2.42) it can be transformed to

$$ \text{vec}(\Delta(S(\theta^*)^T H)) = (H^T \otimes I_s)P_{t,s} \text{vec}(\Delta S(\theta^*)) + I_{qr} \otimes S(\theta^*)^T \text{vec}(\Delta H), \quad (3.15) $$

which yields the partial sensitivities of the residual with respect to the reference null space
and with respect to the tested Hankel matrix

$$
\begin{align}
\mathcal{J}_{\tilde{\zeta},s} &= (\mathcal{H}^T \otimes I_s)\mathcal{P}_{t,s} \\
\mathcal{J}_{\tilde{\zeta},\mathcal{H}} &= I_{qr} \otimes S(\theta^*)^T. \tag{3.16}
\end{align}
$$

t and s are the number of rows and columns of S(θ*), and J̃* is the sensitivity that was
used for the covariance of the conventional residual in (2.37).

The sensitivity $\mathcal{J}_g$ in (3.13) is related to the Hankel matrices $\mathcal{H}^*$ and $\mathcal{H}$. Thus, with (3.16) the link between the reference matrix $S(\theta^*)$ and the Hankel matrix in the reference state $\mathcal{H}^*$ is still missing. It is denoted by the sensitivity $\mathcal{J}_{S,\mathcal{H}^*}$ and $\mathcal{J}_g$ becomes

$$
\mathcal{J}_g = [\mathcal{J}_{\tilde{\zeta},S} \ \mathcal{J}_{S,\mathcal{H}^*} \quad \mathcal{J}_{\tilde{\zeta},\mathcal{H}}]. \qquad (3.17)
$$

$\mathcal{J}_{S,\mathcal{H}^*}$ depends on the computation procedure of the reference matrix. The reference null space is computed with a SVD either from the Hankel matrix or from the observability matrix in the modal basis, where the identification of the modal parameters is required in the reference state (see section 2.4.1). The respective derivatives are derived in the following two paragraphs, and finally the computation of the residual covariance estimate is presented.

**Sensitivity of a reference null space computed from Hankel matrix** When the reference matrix is computed from the Hankel matrix directly, the relation between them is defined through the SVD of $\mathcal{H}^*$ in accordance to (2.10). The reference matrix is the left null space $U_2$ of the Hankel matrix. This relation is used to derive the corresponding sensitivity $\mathcal{J}_{S,\mathcal{H}^*}$.

In a first step, the perturbation $\Delta \mathcal{H}^*$ is propagated to its column space $U_1$ in the SVD.
This was published in [LLM08] and writes as

$$
\Delta U_1 = U_1 R + U_2 U_2^T \Delta \mathcal{H}^* V_1 D_1^{-1}. \quad (3.18)
$$

Matrices $U_1$, $U_2$, and $V_1$ are the column space, the left null space, and the row space, respectively. $D_1$ is the matrix with singular values, which expected values are nonzero. Matrix $R$ will be eliminated later. Applying matrix operation (2.40) it follows a vectorized
form with

$$ \mathrm{vec}(\Delta U_1) = (I_n \otimes U_1) \mathrm{vec}(R) + (D_1^{-1} V_1^T \otimes U_2 U_2^T) \mathrm{vec}(\Delta \mathcal{H}^*). \quad (3.19) $$

This vectorized perturbation of the column space is propagated to the left null space $U_2$ as follows. First it is considered that $U_1^T U_2 = 0$ holds. Introduce a perturbation $\Delta(\cdot)$ to this yields $\Delta U_1^T U_2 + U_1^T \Delta U_2 = 0$. With identity matrices $I$ as defined in section 2.5.1, permutation matrix $\mathcal{P}$ (2.42), and (2.40) to (2.43) this can be vectorized and converted to

$$ (I_s \otimes U_1^T) \mathrm{vec}(\Delta U_2) = -(U_2^T \otimes I_n) \mathcal{P}_{t,n} \mathrm{vec}(\Delta U_1). \quad (3.20) $$

Then, with $U_2^T U_2 = I$ and thus $\Delta(U_2^T U_2) = 0$, it holds $\Delta U_2^T U_2 + U_2^T \Delta U_2 = 0$ leading with the permutation matrix $\mathcal{P}$ (2.42) and (2.40) to (2.43) to

$$ \mathcal{P}_{s,s}(I_s \otimes U_2^T) \mathrm{vec}(\Delta U_2) + (I_s \otimes U_2^T) \mathrm{vec}(\Delta U_2) = 0. \quad (3.21) $$

A particular solution for $\mathrm{vec}(\Delta U_2)$ follows from (3.20) and (3.21) to

$$ \begin{aligned} \mathrm{vec}(\Delta U_2) &= -(I_s \otimes U_1)(U_2^T \otimes I_n) \mathcal{P}_{t,n} \mathrm{vec}(\Delta U_1) \\ &= -\mathcal{P}_{s,t}(U_1 \otimes U_2^T) \mathrm{vec}(\Delta U_1), \end{aligned} \quad (3.22) $$

using properties of the permutation matrix $\mathcal{P}$ in (2.42).

The sensitivity of the left null space to its source matrix is derived by substituting (3.19) into (3.22). This leads to

$$ \mathrm{vec}(\Delta U_2) = -\mathcal{P}_{s,t}(U_1 D_1^{-1} V_1^T \otimes I_s) (I_{qr} \otimes U_2^T) \mathrm{vec}(\Delta \mathcal{H}^*). \quad (3.23) $$

And with (2.43) it is

$$ \mathcal{J}_{U_2, \mathcal{H}^*} = -\mathcal{P}_{s,t}(U_1 D_1^{-1} V_1^T \otimes U_1^T). \quad (3.24) $$

Now coming back to the reference null space as left null space of the Hankel matrix $U_2 = S(\theta^*)$. Since the asymptotic covariance can be derived at $\mathcal{H} = \mathcal{H}^*$ in the reference
state, with (3.24) and (3.16) it is

$$
\begin{align}
\mathcal{J}_{\tilde{\zeta}, \mathcal{H}^*} &= \mathcal{J}_{\tilde{\zeta}, S} \mathcal{J}_{S, \mathcal{H}^*} \nonumber \\
&= (\mathcal{H}^T \otimes I_s) \mathcal{P}_{t,s} (-\mathcal{P}_{s,t}) (U_1 D_1^{-1} V_1^T \otimes S(\theta^*)^T) \nonumber \\
&= -\mathcal{H}^T U_1 D_1^{-1} V_1^T \otimes S(\theta^*)^T \nonumber \\
&= -V_1 V_1^T \otimes S(\theta^*)^T. \tag{3.25}
\end{align}
$$

When the reference null space is computed from the Hankel matrix, sensitivity $\mathcal{J}_g$ in
(3.17) writes as

$$
\mathcal{J}_g = [\mathcal{J}_{\tilde{\zeta}, \mathcal{H}^*} \quad \mathcal{J}_{\tilde{\zeta}, \mathcal{H}}] = [-V_1 V_1^T \otimes S(\theta^*)^T \quad I_{qr} \otimes S(\theta^*)^T]. \quad (3.26)
$$

**Sensitivity of a reference null space computed from observability matrix** Deriving the reference null space from the observability matrix in modal basis includes several steps: First the modal parameters are identified from the Hankel matrix $\mathcal{H}^*$. Let them be collected in parameter vector $\theta^*$. From this the parametric observability matrix $\mathcal{O}(\theta^*)$ is built. In a third step, the left null space of $\mathcal{O}(\theta^*)$ is computed with SVD. The sensitivity of the reference matrix to the Hankel matrix $\mathcal{J}_{S,\mathcal{H}^*}$ is derived by propagate the perturbation of the Hankel matrix, respectively. It is defined as $\mathcal{J}_{S,\mathcal{H}^*} = \mathcal{J}_{S,O} \mathcal{J}_{O,\theta^*} \mathcal{J}_{\theta^*,\mathcal{H}^*}$, where all derivatives are already known from the sections before.

$\mathcal{J}_{S,O}$ is the sensitivity of a left null space $U_2 = S(\theta^*)$ with respect to its source matrix, as in (3.24). Here the source matrix is the parametric observability matrix, and $U_1, D_1$, and $V_1$ are the respective matrices from the SVD of $\mathcal{O}(\theta^*)$. The derivative $\mathcal{J}_{O,\theta^*}$, denoting the sensitivity of the parametric observability matrix with respect to these modal parameters, is defined in (2.63). The uncertainties in the modal parameters that are induced by the SSI are propagated by $\mathcal{J}_{\theta^*,\mathcal{H}^*}$, composed from (2.53), (2.56), and (2.58) to (2.60).

When the reference null space is computed from the parametric observability matrix,
sensitivity $\mathcal{J}_g$ in (3.17) writes with $\mathcal{J}_{\tilde{\zeta},S}$ and $\mathcal{J}_{\tilde{\zeta},H}$ from (3.16) as

$$
\mathcal{J}_g = [\mathcal{J}_{\tilde{\zeta}, \mathcal{H}^*} \quad \mathcal{J}_{\tilde{\zeta}, \mathcal{H}}] = [\mathcal{J}_{\tilde{\zeta}, S} \mathcal{J}_{S, O} \mathcal{J}_{O, \theta^*} \mathcal{J}_{\theta^*, \mathcal{H}^*} \quad \mathcal{J}_{\tilde{\zeta}, \mathcal{H}}]. \quad (3.27)
$$

**Residual covariance estimate** From the formulations of the derivative $\mathcal{J}_g$ in (3.26) or (3.27) the residual covariance $\tilde{\Sigma}_{\zeta}$ can be computed as

$$
\tilde{\Sigma}_{\zeta} = c\Sigma_{1} + \Sigma_{2}, \qquad (3.28)
$$
with the two partial covariances

$$
\begin{align}
\Sigma_1 &= \mathcal{J}_{\zeta, \mathcal{H}^*} \Sigma_{\mathcal{H}} \mathcal{J}_{\zeta, \mathcal{H}^*}^T \nonumber \\
\Sigma_2 &= \mathcal{J}_{\zeta, \mathcal{H}} \Sigma_{\mathcal{H}} \mathcal{J}_{\zeta, \mathcal{H}}^T. \tag{3.29}
\end{align}
$$

$\Sigma_1$ is related to the uncertainty in $\hat{S}(\theta^*)$, and $\Sigma_2$ to the uncertainty of $\hat{\mathcal{H}}$ in the test state.
$\Sigma_{\mathcal{H}}$ is the Hankel matrix covariance, computed in the reference state.

An estimate of the residual covariance can be obtained, when instead of the theoretical Hankel matrix covariance $\Sigma_{\mathcal{H}}$ a consistent estimate $\hat{\Sigma}_{\mathcal{H}}$ (2.36) is used for the partial covariances in (3.29). With data sets of length $M$ for the reference setup and length $N$ in the testing state, the estimated residual covariance is

$$
\hat{\Sigma}_{\zeta} = \frac{N}{M} \hat{\Sigma}_{1} + \hat{\Sigma}_{2}. \tag{3.30}
$$

It follows that for long data sets in the reference state in relation to the data length *N* in the testing state, the effect of the uncertainty of the testing Hankel matrix in $\hat{\Sigma}_2$ to the complete residual covariance is much stronger, than the contribution of the uncertainties related to the reference null space in $\hat{\Sigma}_1$. This is reflected in a small ratio *c* = $N/M$. In this case the testing Hankel matrix is the more uncertain part of the residual. On the other hand, when for the reference setup very few data is used compared to the amount of data in the testing state, *c* becomes large and the uncertain reference has a high impact to the residual covariance.

3.2.4 Test statistic

Analogous to (2.29) the damage detection test for the new residual writes in the parametric case as

$$
\tilde{t} = \tilde{\zeta}(\theta^*, \mathcal{Y}_N)^T \tilde{\Sigma}_\zeta (\theta^*)^{-1} \mathcal{J}_\zeta(\theta^*) \left( \mathcal{J}_\zeta(\theta^*)^T \tilde{\Sigma}_\zeta (\theta^*)^{-1} \mathcal{J}_\zeta(\theta^*) \right)^{-1} \mathcal{J}_\zeta(\theta^*)^T \tilde{\Sigma}_\zeta (\theta^*)^{-1} \tilde{\zeta}(\theta^*, \mathcal{Y}_N), \quad (3.31)
$$

with non-centrality parameter $\tilde{\lambda}$ in the damaged case defined as

$$
\tilde{\lambda} = \delta^T (\mathcal{J}_{\zeta}(\theta^*)^T \tilde{\Sigma}_{\zeta} (\theta^*)^{-1} \mathcal{J}_{\zeta}(\theta^*)) \delta. \quad (3.32)
$$
The non-parametric version of the test becomes equivalent to (2.33)

$$ \tilde{t} = \tilde{\zeta}(\theta^*, \mathcal{Y}_N)^T \tilde{\Sigma}_{\zeta}(\theta^*)^{-1} \tilde{\zeta}(\theta^*, \mathcal{Y}_N), \quad (3.33) $$

and non-centrality parameter $\tilde{\lambda}$ is

$$ \tilde{\lambda} = \gamma^T \tilde{\Sigma}_{\zeta}(\theta^*)^{-1} \gamma. \quad (3.34) $$

The test statistic is consistently computed when $\mathcal{J}_{\zeta}(\theta^*)$ and $\tilde{\Sigma}_{\zeta}(\theta^*)$ are replaced by their respective estimates.

## 3.3 Proof of concept and discussion

In this section the previous developments on theory are validated on numerical simulation data and the results are discussed. The statistical distribution properties of the new residual are evaluated by drawing conclusions from the distribution properties of $\tilde{t}$ (3.31). Since $\tilde{t}$ is a scalar quantity, it can be evaluated easily, e.g. with histograms from a Monte Carlo study. The performance of the new test is represented by the probability of detection (POD) and the receiver operating characteristic (ROC) curves that deliver information on the ratio between POD and PFA for thresholds computed for different PFA limits.

With the main goal to emphasize the impact of the uncertainty in the reference null space $\hat{S}(\theta^*)$, results from the new test (3.31) are compared to those from the conventional test (2.29), which neglects these uncertainties, i.e. $\Sigma_1 = 0$. In this context, the effect of different data lengths is analyzed, since they have an impact on the accuracy of the estimated matrices and consequently the residual covariance matrix (3.30). However, the estimation of the asymptotic Hankel matrix covariance $\hat{\Sigma}_H(\theta^*)$ is not in the focus of this thesis and is assumed to be consistently estimated from (2.36). Other effects are examined, such as the choice of system parametrization, which is interesting when the system is parametrized with respect to its structural parameters for damage localization with SSDD or when a fully automated damage detection procedure uses the simpler unparametrized test version (3.33).

In the following, the setup of the numerical Monte Carlo study is introduced first. Then, the statistical properties and test performance is evaluated. This is followed by four subsections on the effects of data length $M$, data length $N$, system parametrization,
and null space computation procedure.

### 3.3.1 Study setup

Previously, the theory of the stochastic subspace-based damage detection was derived when the reference null space is estimated from data and thus afflicted with uncertainties. The findings are applied in the following on numerical simulation data. A simple system is chosen, which allows for a clear definition of the correct theoretical values to which the results from the numerical study are compared.

Simulation data of the mass-spring-damper system in Fig. 3.1 is used. The system is defined by eight masses of $m_1 = m_3 = m_5 = m_7 = 1$ and $m_2 = m_4 = m_6 = m_8 = 2$. The stiffnesses are set to $k_1 = k_3 = k_5 = k_7 = 200$ and $k_2 = k_4 = k_6 = k_8 = 100$. In the damaged case the stiffness at element 3 is decreased. Classical damping with a damping ratio of 2% is considered for all modes. The system is excited with random white noise at elements 1, 3, 5, and 7. At the same positions the velocities are recorded sampled at 20 Hz. White measurement noise having 5% standard deviation of the signal is added.

<figure><img src="image_12.png" /><figcaption>Figure 3.1 – Mass-spring chain with four sensors.</figcaption></figure>

For the evaluation of the residual behavior a Monte Carlo study with 1,000 repetitions is performed. The test data and also the data for the reference null space is generated anew in each repetition. The reference null space is computed directly from the Hankel matrix estimate in the reference state.

Let the system be parametrized with respect to its modal parameters. Then the parameter vector is $\theta = [\Lambda^T \text{vec}(\Phi)^T]^T$. In this case, the asymptotic residual sensitivity $J_{\zeta}(\theta^*)$ is computed in accordance to (2.28). Assuming the sensitivity is of full rank, it holds $\text{rank}(J_{\zeta}(\theta^*)) = 2mr = 64$, with $m = 8$ degrees of freedom of the mechanical system and $r = 4$ sensors. Thus, theoretically the test statistics follows a $\chi^2$-distribution with 64 DOF (2.31) and the expected test value mean in the reference state is 64. It should be noted that the test values in the damaged state depend on the non-centrality parameter (3.32), that contains the asymptotic covariance matrix estimate. Consequently, no theoretically expected value can be derived for the test statistics of the damaged system.
3.3.2 Statistical properties and test performance

The test values in the reference state and the corresponding theoretical χ²-distribution
are presented in Fig. 3.2. The null space is estimated from data with length *M* = 50,000
and the data set for the test Hankel matrix is of length *N* = 100,000. Obviously, consid-
erable uncertainties related to the reference null space are neglected in the conventional
test (2.29) with Σ₁ = 0. The conventional test statistics differs clearly from the theoret-
ical value. The central limit theorem in (2.27) does not apply, since it is formulated for
a residual with perfectly known reference null space. A bias is introduced in the residual
leading to a non-zero mean value *b* ≠ 0 in the reference state and consequently to a non-
centrality of the test distribution already in the reference state. The resulting test values
are significantly higher than the theoretical value. This implies also a larger spread due to
the χ²-distribution properties (Fig. 2.1, right) reflecting the uncertainty in the reference
null space. When the uncertainty in the reference null space is considered correctly with
the new test (3.31) the values are well centered around the expected value and the spread
follows the theoretical χ²-distribution with 64 DOF.

<figure><img src="image_13.png" /><figcaption>Figure 3.2 – Theoretical χ²-distribution and histograms of the new test statistics and the conventional test with Σ₁ = 0, M = 50,000, N = 100,000.</figcaption></figure>

In the damaged state, the distributions in both tests are shifted from the reference val-
ues as demonstrated in Fig. 3.3. The damage is simulated by decreasing the third element
stiffness by 2%. The shift reflects the non-centrality parameter $\lambda$ in the conventional test
or $\tilde{\lambda}$ (3.32) in the new test. Notice, that due to neglecting the uncertainty of the reference
null space the non-centrality parameter $\lambda$ of the conventional test is affected by the for-
mer described bias, too. It cannot be predicted by equation (2.31) in presence of this bias.
With the conventional test (left) the values have wide spread in both states resulting in
a significant overlap. In contrary, the new test (right) leads to narrow distributions of the
test statistics and a clear separation of both states.
<figure><img src="image_14.png" /><figcaption>Figure 3.3 – Histogram in the reference and the damaged state of the conventional (left) and the new test (right), $M = 50,000$, $N = 100,000$, damage = 2%.</figcaption></figure>

Two thresholds for damage detection are shown in Fig. 3.3, derived as explained in Fig. 2.1 (left). An empirical threshold that is based on the distribution of the test values in the undamaged state received by Monte Carlo simulations. The other threshold is derived from the theoretical $\chi^2$-distribution in the reference state. The thresholds are chosen, allowing for a maximum probability of false alarms (PFA) of 0.1%. In cases where the test value exceeds this threshold, damage is assumed. For the conventional test in the left plot, the empirical threshold is much larger than the theoretical one, similar to the reference test values. Using the theoretical threshold will lead to 100% PFA, since all reference values are above the threshold. When the uncertainties of the reference null space are considered correctly in the new test (right), the empirical threshold coincides with the theoretical one. This allows for threshold computation a priori that is purely based on theory.

From the distributions in the reference and the damaged state, the probability of detection (POD) and the receiver operating characteristic (ROC) curves are derived in Fig. 3.4 demonstrating the impact of the new test computation on the test performance. The left plot shows the POD of both tests computed with the empirical thresholds and data length $M = 50,000$. The conventional test is not very sensitive towards small damages. This is due to the overlap between the two states that is apparent in Fig. 3.3 (left). Bigger damages, i.e. big change vectors $\delta$ and thus big non-centrality parameters, lead to a better separation of the states and consequently to a better POD. The new test performs well even with very small changes in the stiffness, due to the well comprised uncertainties of the reference null space and thus good separation of the system states (Fig. 3.3, right).

Regarding the relation between PFA and POD, the ROC curves in Fig. 3.4 (right)
<figure><img src="image_15.png" /><figcaption>Figure 3.4 – POD (left) and ROC diagram with damage of 2% (right) of the conventional and the new test, $M = 50,000$, $N = 100,000$.</figcaption></figure>

for the conventional and the new test can be derived. To computed them, the maximum allowed number of false alarms is varied leading to a corresponding POD and PFA. The ROC curve of the conventional test is much closer to the line of no performance, while the new test approaches the best performance point. In the latter the POD is less depending on the allowed PFA.

It has to be emphasized that the conventional test may show unpredictable test performance due to the former described bias in the residual. In the reference state a non-centrality exists due to the residual mean $b$, and in the damaged state the non-centrality parameter is a combination of terms with $b$ and change vector $\delta$, including mixed terms. Since this mixed terms are not necessarily positive, this may lead to a bigger or a smaller non-centrality parameter in one or the other case.

In the above evaluation it could be shown that the new test with correct consideration of the uncertainties related to $\hat{S}(\theta^*)$ is in accordance to the theoretically expected behavior. For this test, a threshold for damage detection can be computed on a theoretical basis. The test is sensitive towards small damages and show a very good POD-PFA relation. In the following, several effects are analyzed regarding their impact on the test behavior.

### 3.3.3 Data length in the reference state

It can be assumed that the bias and the uncertainties of the reference null space decrease when better estimates of $\hat{S}(\theta^*)$ are achieved from longer data sets approaching the theoretical value $\hat{S}(\theta^*) \to S(\theta^*)$. This is illustrated by the empirical mean values of both tests, plotted with respect to the reference data length $M$ in Fig. 3.5. The conventional test that neglects the uncertainty related to the reference null space approaches the values
<figure><img src="image_16.png" /><figcaption>Figure 3.5 – Empirical mean of conventional and new tests for different data lengths M, $N = 100,000$.</figcaption></figure>

of the new test and thus the theoretical value for longer data sets, i.e. when the reference null space is well estimated. The test statistics of the new test is mostly independent from the data length and is in accordance to the theoretical value, since the uncertainties are considered correctly in any case.

The distributions of both tests in the reference state for three different data lengths $M$ in Fig. 3.6 correspond to these findings. The reference test values of the conventional test (left) approach the theoretical value with less spread for larger $M$, since the bias in the residual mean decreases for longer data sets due to less uncertainty in the reference null space. When the correct distribution properties are considered in the new test (right), the results in the reference state are nearly independent from the reference data length and coincide with the theoretical values for any $M$. This is in particular interesting for the threshold computation: While for the conventional test the threshold depends on the data length $M$ in the reference state, the new test has one fix (theoretically based) threshold.

However, in the damaged state the test values depend on $M$ for both tests. In the conventional test (left) the damaged test values decrease for bigger $M$ and have less spread. This is due to less uncertainty in $\hat{S}(\theta^*)$ and consequently smaller bias, similar as for the test values in the reference state. The shift between the mean of both states is similar for any $M$. As discussed before, this shift does not correspond to the non-centrality parameter in (2.31) for small $M$ due to the bias in the residual mean, and is unpredictable. The test values of the new test in the damaged state increase with larger $M$. Less uncertainty on the reference null space is reflected in the decreasing covariance contribution $\frac{N}{M}\Sigma_1$ for larger $M$. This leads qualitatively to a smaller residual covariance $\tilde{\Sigma}_{\zeta}$ in (3.30), a larger inverse $\tilde{\Sigma}_{\zeta}^{-1}$, and thus a larger non-centrality parameter $\tilde{\lambda}$ in (3.32).
<figure><img src="image_17.png" /><figcaption>Figure 3.6 – Histogram in the reference and the damaged state of the conventional (left) and the new test (right) for different *M*, *N* = 100,000, damage = 2%.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 3.7 – POD of conventional and new test for different *M*, *N* = 100,000, damage = 1.4%.</figcaption></figure>

The resulting test values in the damaged state become larger.

In both tests, the effect of length *M* leads to a better separation of the damaged from the reference state. Hence, the POD of both tests improves for larger reference data length *M* as shown in Fig. 3.7 for a damage of 1.4%.

When the reference null space is well estimated from a long data set with *M* = 500,000, the conventional and the new test perform more similar, due to small uncertainty related to the null space. This is demonstrated in Fig. 3.8.

It can be concluded that for a well estimated reference null space, the conventional test leads to results similar to the new test. Both tests have a better test performance for long reference data sets. The new test, however, does not rely on *M* in the reference
<figure><img src="image_19.png" /><figcaption>Figure 3.8 – POD (left) and ROC diagram with damage of 2% (right) of the conventional and the new test, $M = 500,000$, $N = 100,000$.</figcaption></figure>

state allowing for a theory-based threshold for damage detection and a predictable test performance.

### 3.3.4 Data length of the testing data set

So far, the effect of the data length $M$ in the reference state was analyzed, affecting the accuracy of $\hat{S}(\theta^*)$ directly. In the following the impact of the testing data length $N$ will be evaluated.

Fig. 3.9 presents the histograms of both tests for a fix damage of 2% and data length $M = 100,000$ for three different testing data lengths $N$. The reference values of the conventional test (left) increase with larger $N$. As in the context of the investigations to data length $M$, this is due to the bias that is introduced to the conventional residual by neglecting the uncertainties in $\hat{S}(\theta^*)$ and the resulting mean $b$ of the residual in the reference state. Since the residual is normalized to testing data length $N$, $b$ depends on $N$ and so the resulting non-centrality parameter in the reference state does. This leads to larger reference test values for longer testing data sets. A damage detection threshold for the conventional test must be computed on empirical distribution properties. Moreover, the testing data length $N$ cannot be changed at any time, since the reference test values depend on it. The new test (right) yields reference test values that are in good accordance to the theoretical values and independent from $N$ allowing for the definition of a fixed threshold based on theory.

In the damaged case, both tests are affected by the testing data length $N$. For longer test data sets but under the same physical parameter change $(\theta - \theta^*)$, the change vector $\delta$ increases as defined in the hypotheses in (2.26). This leads to larger non-centrality pa-
<figure><img src="image_20.png" /><figcaption>Figure 3.9 – Histogram in the reference and the damaged state of the conventional (left) and the new test (right) for different $N$, $M = 100,000$, damage = 2%.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 3.10 – POD of conventional and new test for different $N$, $M = 100,000$, damage = 2%.</figcaption></figure>

parameter $\lambda$ or $\tilde{\lambda}$ (3.32). As a consequence, smaller damages can be detected with increasing data length $N$, which is a well known property of the SSDD. This is reflected in Fig. 3.10, where the POD for 2% damage is plotted for different data length $N$.

These results show that the testing data length $N$ affects the values of the conventional test in the reference and the damaged state, and changes in data length $N$ perturb the setup of even an empirical threshold. In the new test the effect of data length $N$ is in accordance to theory, where longer testing data sets increase the non-centrality parameter and thus improve the POD, while $N$ has no impact on the reference values and thus on the computation of a damage detection threshold.
### 3.3.5 System parametrization

It has been shown that for a system parametrized with respect to its modal parame-
ters the new test statistics follows the theoretical expected distributions and that the test
performance is reliable. However, different system parametrizations may be used. This is
the case, when e.g. the SSDD is applied to damage localization or quantification [DMZ16]
and a system parametrized with respect to its structural parameters is required. Or the
damage detection tests are applied without focusing on a particular system parametriza-
tion $\theta$, which corresponds to a direct detection of changes in the residual and $\mathcal{J} = I$
[DMH14; MDV21]. These non-parametric tests are attractive due to their simplicity and
may be preferred e.g. when the identification of the modal parameters in the reference
state of the system is inconvenient. In the following the consistency of the new residual
for different system parametrizations is examined.

Assume the system in Fig. 3.1 is parametrized with respect to its structural parameters, namely to its eight stiffnesses such that $\theta = [\mathbf{k}_1 \mathbf{k}_2 \dots \mathbf{k}_8]^T$. Than the theoretical test value is 8. In Fig. 3.11 the corresponding test values from the conventional test with $\Sigma_1 = 0$ and the new test are presented for different data lengths $M$. The reference test values of the new test are nearly independent from $M$, while the conventional test approaches the theoretical value only for long reference data sets yielding low uncertainty of the reference null space. This is comparable to the modal parametrized system in Fig. 3.5.

The POD of both tests is presented in Fig. 3.12 computed with an empirical threshold,
data lengths $N = 100,000$, $M = 50,000$ (left), and $M = 500,000$ (right). The POD is
affected by data length $M$ similarly to the modal parametrized system: The performance
of the new test is better than of the conventional test for small $M$ (left) and the POD
of both tests is the same when $\hat{S}(\theta^*)$ is well estimated from long data sets (right). Both
tests are more sensitive when less uncertainties are related to the reference null space in
case of large $M$ due to the effect explained in the context of Fig. 3.6.

In the non-parametric case, the evaluation of the mean value towards a theoretical
value is omitted in the context of the numerical study, since the theoretical value cannot
be determined. It is the theoretical rank of the residual, which cannot be derived clearly
from an estimate of the covariance. Only an upper bound for the number of degrees of
freedom of the $\chi^2$-distribution can be found from the covariance matrix size. Similarly
to the parametric cases, the POD in Fig. 3.13 computed with the same data settings as
before is better for the new test compared to the conventional test. Both tests show a
better and closer performance for large $M$.
<figure><img src="image_22.png" /><figcaption>Figure 3.11 – Empirical mean of conventional and new tests for different data lengths M for a system parametrized with respect to structural parameters, *N* = 100,000.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 3.12 – POD of conventional and new test for a system parametrized with respect to structural parameters for *M* = 50,000 (left) and *M* = 500,000 (right), *N* = 100,000.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 3.13 – POD of conventional and new test with non parametric test for $M = 50,000$ (left) and $M = 500,000$ (right), $N = 100,000$.</figcaption></figure>

When comparing the test performance of the tests with different parametization, the POD in Figs. 3.4, 3.8, 3.12, and 3.13 decreases with higher DOF of the theoretical $\chi^2$-distribution. This is coherent with the findings in [MDV21] and is due to the wider spread
of the corresponding distribution. It is not an impact of the new residual.

Independently from the parametrization, the behavior of the new test is in accordance
to theory. For any parametrization, the new test yields a better POD than the conventional
test in presence of considerable uncertainty in the reference null space.

3.3.6 Reference null space computed from observability matrix

The previous evaluations are based on test values that are computed with a reference
null space derived from the covariance Hankel matrix in the reference state. It is possible
to compute the reference matrix with a SVD from the parametric observablity matrix
O(θ*), analogous to (2.10). O(θ*) is built from the m conjugate complex pairs of the
system modes (2.21). Thus, the model order n = 2m is known, leading to a left null space
that is correctly separated from the column space.

The test value means of the conventional and the new test with a reference null space computed from the parametric observability matrix are presented in Fig. 3.14. The test values show good accordance to Fig. 3.5, where the new test matches the theoretical value even for short data and is independent from the data length. The conventional test depends on the data length and approaches the theoretical value only for longer data sets when less uncertainty is related to the reference null space. However, the test values approach the theoretical values faster than in the former cases. It is concluded that a reference null space, which is computed from the parametric observability matrix, is less afflicted with uncertainties, leading to a smaller bias in the conventional residual.

The test performance of the conventional and the new test with a reference null space computed from $O(\theta^*)$ is shown in Fig. 3.15 for $M = 50,000$ (left) and $M = 500,000$

<figure><img src="image_25.png" /><figcaption>Figure 3.14 – Empirical mean of conventional and new tests for different data lengths M when the reference null space is computed from observability matrix, N = 100,000.</figcaption></figure>
(right), and $N = 100,000$. While the null space computation procedure has no impact on the POD of the new test (cf. Figs. 3.4 and 3.8), the conventional test performs better than before when the reference null space was computed from the covariance Hankel matrix directly. In fact, the POD of the conventional test coincides with the results of the new test for any data length. This corresponds to the conclusion from above, that a null space computed from the observability matrix leads to smaller bias, which is here reflected in a similar POD. And it confirms the remarks that are made in the context of Fig. 3.4 regarding the unpredictable non-centrality and POD of the conventional test.

<figure><img src="image_26.png" /><figcaption>Figure 3.15 – POD of conventional and new test when the reference null space is computed from observability matrix for $M = 50,000$ (left) and $M = 500,000$ (right), $N = 100,000$.</figcaption></figure>

It could be shown that the results of the new test are in good accordance to the theoretical expectations and that the new test performs reliably, for any null space computation procedure. The conventional test, however, performs unpredictable for different null space computations.

## 3.4 Summary

In this section the theory for a covariance-driven subspace-based residual with an uncertain reference has been derived. Previously, only statistical uncertainties related to test data have been considered in the formulation of the residual, and the reference matrix was assumed to be perfectly known. However, this is not true in real world applications when the reference is estimated from data. In this case the estimation uncertainties stemming from reference data must be taken into account for the residual and the associated statistical test. Now, the residual has been designed as a function of the uncertain reference and the testing data. The central limit theorem has been applied to the two random
variable residual and the statistical distribution properties have been derived. With the
statistical delta method the perturbations from both variables were propagated to the
residual, leading to its asymptotic residual covariance. The corresponding statistical test
has been formulated.

Correct and predictable statistical behavior of the residual and the new test could
be shown by means of application to numerical data. It has been demonstrated that the
estimation uncertainties related to the reference data are non-negligible especially when
the reference data is of short length.

Neglecting these uncertainties can lead to a considerable deviation from the theoreti-
cally expected test behavior. As a consequence, the conventional damage detection relies
on empirical evaluated distribution properties, and the test behavior is not predictable
as long as considerable uncertainty is related to the reference null space. It depends e.g.
on the choice of null space computation. Stable results can only be achieved from suffi-
cient long reference data sets leading to well estimated null spaces. Parameters such as
the testing length *N* cannot be changed during the monitoring phase, since it affects the
reference values. All this can lead to false alarms and makes it impossible to set thresholds
between healthy and damaged states based on the expected statistical distribution of the
test.

In contrary, the new test considering the uncertainties of the reference null space
correctly makes the damaged detection procedure predictable and shows good accordance
with the theoretically expected statistical behavior. The results are not affected by the
choice of system parametrization nor by the null space computation procedure. Thresholds
for damage detection can be derived from theoretical distribution properties a priori in
the reference state and are not depending on any of the investigated parameters.

Results from an application of these developments to experimental data can be found
in chapter 5 and 7.

**3.5 Dissemination**

Parts of this chapter have been published in

[Vie+17]: E. Viefhues, M. Döhler, F. Hille, and L. Mevel. « Stochastic subspace-based damage detection with uncertainty in the reference null space ». In: *IWSHM - 11th International Workshop on Structural Health Monitoring*. Stanford, United States, 2017
[Vie+18]: E. Viefhues, M. Döhler, F. Hille, and L. Mevel. « Asymptotic analysis of subspace-based data-driven residual for fault detection with uncertain reference ». In: *IFAC-PapersOnLine* 51.24 (2018). 10th IFAC Symposium on Fault Detection, Supervision and Safety for Technical Processes SAFEPROCESS 2018, pp. 414–419

[Vie+22]: E. Viefhues, M. Döhler, F. Hille, and L. Mevel. « Statistical subspace-based damage detection with estimated reference ». In: *Mechanical Systems and Signal Processing* 164 (2022)
<figure><img src="image_27.png" /></figure>
# STOCHASTIC SUBSPACE-BASED
DAMAGED DETECTION UNDER CHANGING
ENVIRONMENTAL CONDITIONS

## 4.1 Introduction

The application of automated vibration-based damage diagnosis to structural health monitoring (SHM) of civil engineering structures often fails due to changing environmental conditions, since the dynamical behavior of the structures is affected [Soh07] leading to high false alarm rates. The subspace-based damage detection (SSDD) method (section 2.4 and chapter 3) is based on a residual that is sensitive to changes in the dynamical system. So far it is difficult to classify if these changes are due to damages or result from environmental effects. In this chapter a SSDD approach robust to such changes in the environmental conditions and suitable for automated methods is derived.

The developments are motivated by the model interpolation approach for linear pa-
rameter varying (LPV) systems from [Zha18] introduced in section 2.6. A mechanical
system that is affected by e.g. temperature, can be interpreted as a LPV system depend-
ing on an external, measurable parameter. To derive an appropriate reference null space
for the residual in (3.2) that is valid at a certain (testing) parameter, a global state-space
model is computed. It is based on local reference models, which are obtained from data at
some parameter points in the reference state. The global model describes the system with
its dynamical properties at the parameter point of interest, i.e. the testing temperature.

The new damage detection approach is appropriate for automated vibration based
SHM with high robustness to the effects of a varying measurable environmental parameter.
The method reduces the probability of false alarms and is sensitive to small damages. The
notation in this chapter refers to the residual from chapter 3 that takes the uncertainties
of an interpolated reference null space into account, but the developments hold for the
conventional residual from section 2.4 as well.

In the following, the theory for SSDD with a reference null space that is adapted to the current testing temperature is derived pursuing three different strategies in section 4.2 and the associated testing procedure is presented including the covariance computation for a residual with interpolated reference null space. In a numerical study in section 4.3 the concept is evaluated regarding its performance. The results are compared to the performance of an existing approach from [Bal+08a] that defines a global reference model built from a merged data set smoothing out the temperature effect by the averaging operation.

## 4.2 Damage detection robust to temperature effects

In this section the theory for the subspace-based damage detection method robust to changes in the environmental conditions is derived. First, the formalized problem statement is given. The theory for a parameter adapted reference null space for a covariance driven subspace-based residual follows considering three different computation strategies: computation of a null space from an interpolated output covariance Hankel matrix, computation from an observability matrix built from an interpolated eigenstructure, and a null space computed from an observability matrix that is based on modal parameters derived by means of regression analysis. Afterwards, the residual function with the interpolated reference null space and the associated statistical test procedure are presented and the appropriate residual covariance is derived. The section closes with comparing the principles of the new approach with previous temperature rejection approaches for SSDD.

### 4.2.1 Formalized problem statement

When a mechanical system, such as a civil engineering structure, is affected by an external scalar parameter $P$, such as the temperature, it is represented for any parameter working point $P = P_j$ by a local version of the equation of motion from 2.1 writing as

$$ M_j \ddot{q}_j(t) + C_j \dot{q}_j(t) + K_j q_j(t) = v_j(t). \quad (4.1) $$

Matrices $M_j$, $C_j$, and $K_j \in \mathbb{R}^{m \times m}$ contain the masses, dampings, and stiffnesses at the $m$ degrees of freedom arising at the working point $P = P_j$. In accordance to (2.1) $t$ denotes the continuous time, vector $q_j \in \mathbb{R}^m$ contains the displacements at the degrees of
freedom of the structure, and the excitation force $v_j \in \mathbb{R}^m$ is considered to be unknown. Consequently, the state-space representation at time instance $k$ in (2.3) derived from the mechanical system matrices is similarly affected by the external parameter. This results in an individual state-space model for any parameter working point $P = P_j$ with state transition matrix $A_j$ and observation matrix $C_j$ writing as

$$ \left\{ \begin{aligned} x_{j,k+1} &= A_j x_{j,k} + w_{j,k} \\ y_{j,k} &= C_j x_{j,k} + v_{j,k}. \end{aligned} \right. \qquad (4.2) $$

As defined for (2.3) $x_{j,k} \in \mathbb{R}^n$ is the state vector of the system at time instance $k$, $n = 2m$ is the model order, and $y_{j,k} \in \mathbb{R}^r$ contains the measured outputs. The noise terms $w_{j,k}$ and $v_{j,k}$ are related to the unknown excitation and the measurement noise, and are considered to be stationary and white. This is what is called the *local model* in the model interpolation approach introduced in section 2.6.

To each local model a parameter vector $\theta_j$ with a nominal value $\theta_j^*$ in the reference state is associated. The considered damage detection method described in section 2.4 and further developed in chapter 3 is based on a subspace-based feature vector that evaluates changes in the monitored system with respect to an undamaged reference state statistically on the basis of the parameter vector $\theta_j$ [BAB00]. A local version of the covariance-driven subspace-based damage detection method can be derived for a fixed working point as follows.

The block Hankel matrix in (2.12) is built from theoretical output covariances at $P = P_j$, which are defined as

$$ R_{j,i} = \mathbf{E}(y_{j,k+i} y_{j,k}^T) = C_j A_j^{i-1} G_j, \qquad (4.3) $$

with $G_j = \mathbf{E}(x_{j,k+1} y_{j,k}^T)$. The local Hankel matrix $\mathcal{H}_j \in \mathbb{R}^{(p+1)r \times qr}$ writes as

$$ \mathcal{H}_j \stackrel{\text{def}}{=} \begin{bmatrix} R_{j,1} & R_{j,2} & \dots & R_{j,q} \\ R_{j,2} & R_{j,3} & \dots & R_{j,q+1} \\ \vdots & \vdots & \ddots & \vdots \\ R_{j,p+1} & R_{j,p+2} & \dots & R_{j,p+q} \end{bmatrix}. \qquad (4.4) $$

With the factorization property of the Hankel matrix into the observability matrix $\mathcal{O}_j$
and the controllability matrix $\mathcal{C}_j$ (cf. (2.13))

$$ \mathcal{H}_j = \mathcal{O}_j \mathcal{C}_j, \quad \text{where} \quad \mathcal{O}_j = \begin{bmatrix} C_j \\ C_j A_j \\ \vdots \\ C_j A_j^p \end{bmatrix}, \quad \mathcal{C}_j = [G_j \quad A_j G_j \quad \dots \quad A_j^{q-1} G_j], \quad (4.5) $$

the underlying property in the reference state for the residual formulation in (2.23) and (2.22) becomes

$$ S_j^T \mathcal{H}_j = 0 \quad \text{and} \quad S_j^T \mathcal{O}_j = 0. \quad (4.6) $$

$S_j$ denotes the reference null space for parameter point $P_j$ computed either from the Hankel matrix or from the observability matrix $\mathcal{O}_j$, e.g. by means of the singular value decomposition.

With a data set $\mathcal{Y}_{j,N} = \{y_{j,1}, \dots, y_{j,N}\}$ of length $N$ in the test state, for $P = P_j$ the residual in (3.2) is defined as

$$ \tilde{\zeta}(\theta_j^*, \mathcal{Y}_{j,N}) \stackrel{\text{def}}{=} \sqrt{N} \operatorname{vec}(\hat{S}_j^T \hat{\mathcal{H}}_j). \quad (4.7) $$

The close hypotheses write as

$$ \begin{aligned} \mathbf{H}_0^j &: \theta = \theta_j^* \quad (\text{reference state}), \\ \mathbf{H}_1^j &: \theta \neq \theta_j^* \quad (\text{damaged state}), \end{aligned} \quad (4.8) $$

and the distribution properties of the residual in the non-parametric version correspond to (3.12) with

$$ \tilde{\zeta}(\theta_j^*, \mathcal{Y}_{j,N}) \rightarrow \begin{cases} \mathcal{N}(0, \tilde{\Sigma}_{\zeta,j}) & \text{under } \mathbf{H}_0^j \\ \mathcal{N}(\tilde{\gamma}, \tilde{\Sigma}_{\zeta,j}) & \text{under } \mathbf{H}_1^j, \end{cases} \quad (4.9) $$

where $\tilde{\gamma}$ is a change in the residual mean due to damage, and $\tilde{\Sigma}_{\zeta,j}$ is the asymptotic residual covariance that can be easily estimated from data in the reference state with (3.28) at parameter point $P = P_j$. The associated statistical test is (cf. (3.33))

$$ \tilde{t} = \tilde{\zeta}(\theta_j^*, \mathcal{Y}_{j,N})^T \tilde{\Sigma}_{\zeta,j}^{-1} \tilde{\zeta}(\theta_j^*, \mathcal{Y}_{j,N}). \quad (4.10) $$

The described method works under the premise that the testing data $\mathcal{Y}_{j,N}$ is obtained under the same external parameter as the left reference null space $S_j$ in the definition of
the residual in (4.7). If the Hankel matrix is computed from testing data under a different external parameter than the left reference null space, condition (4.6) of the reference state may not hold anymore, since the external parameter may change the dynamic properties of the system without the presence of damage. Thus, the residual mean is not zero anymore in the reference state under a change in the external parameter leading to false alarms of the damage detection test (4.10).

Suppose testing data is obtained under a known external parameter $P$ that is possibly different from a set of reference parameter points {$P_1, P_2, \dots, P_u$} for which measurements in the reference state of the structure are available. To avoid false alarms, a left reference null space is required that is appropriate for the environmental conditions of the testing data, i.e. for parameter $P$ such that condition (4.6) of the system in the reference state is actually satisfied. It is the goal to develop from the available reference measurements an appropriate reference null space $\underline{S}(P)$ and the associated damage detection test for data under parameter $P$.

### 4.2.2 Parameter adapted reference null space

The idea of computing a reference null space $\underline{S}(P)$ based on weighted references was already used in [FKB13]. In the following, an approach is developed based on system theoretic considerations, where $\underline{S}(P)$ is obtained by means of the local model interpolation approach from [Zha18] as recalled in section 2.6. With the goal to characterize data of a dynamical system at parameter $P$, in this approach a global reference model that is valid at external parameter $P$ is interpolated from local models (4.2) at parameters {$P_1, P_2, \dots, P_u$}. Then, the outputs of the system at parameter $P$ can be considered as an interpolated version of the outputs at parameters {$P_1, P_2, \dots, P_u$} [Zha18]. However, no data is generated explicitly. Instead, the matrices from which the reference null space is derived (2.23) or the underlying parameters for (2.22) are interpolated yielding a temperature adapted reference null space for the subspace-based residual at parameter $P$.

The global model recalled from (2.67) writes as

$$ 
\begin{aligned}
\underline{x}_{k+1} &= \underline{A} \underline{x}_k + \underline{w}_k \\
y_k &= \underset{\underline{C}(P)}{\overset{\text{def}}{\equiv}}} \underline{x}_k + \underline{v}_k,
\end{aligned} 
\quad (4.11) 
$$
with

$$
\begin{align*}
\underline{A} &\stackrel{\text{def}}{=} \begin{bmatrix} A_1 & \cdots & A_u \end{bmatrix}, \quad \underline{x}_k &\stackrel{\text{def}}{=} \begin{bmatrix} x_{1,k} \\ \vdots \\ x_{u,k} \end{bmatrix}, \quad \underline{w}_k &\stackrel{\text{def}}{=} \begin{bmatrix} w_{1,k} \\ \vdots \\ w_{u,k} \end{bmatrix}, \\
\underline{C}(P) &\stackrel{\text{def}}{=} [\rho_1(P)C_1 \ \cdots \ \rho_u(P)C_u], \quad \underline{v}_k &\stackrel{\text{def}}{=} \sum_{j=1}^u \rho_j(P) v_{j,k}(t).
\end{align*}
$$

The interpolation weights $\rho_j(P)$ for each of the $u$ reference models are determined for the current parameter $P$ in dependence of the reference parameters $P_1, \dots, P_u$. The weighting function can be chosen depending on prior knowledge of the qualitative behavior of the mechanical structure under changes in the external paramete, and yields

$$
\sum_{j=1}^{u} \rho_j(P) = 1 \qquad (4.12)
$$

for any $P$ with $0 \le \rho_j(P) \le 1$. For example, a bell-shaped (Gaussian) function centered at $P$ can be chosen, and the weights $\rho_j(P)$, $j = 1, \dots, u$, are the values of this function at $P_j$ after normalization by (4.12). Fig. 4.1 shows how the weights for the local models are obtained from a Gaussian weighting function for any $P$. Wile the local models near to the centered point have a high impact to the global model, the far away located models have a small weighting factor.

In the following three null space computation strategies are presented that are based
on or motivated by the local model interpolation approach from [Zha18]. It is derived how
the interpolated Hankel matrix at parameter $P$ can be obtained by means of the model

<figure><img src="image_28.png" /><figcaption>Figure 4.1 – Determination of weights at parameter working points for any *P* from a Gaussian function.</figcaption></figure>
interpolation approach but without actually generating reference output. With this idea in mind, interpolation of the reference modal parameters is used in the second strategy to obtain the left reference null space from the respective observability matrix, and the third approach pursues this strategy with regression of the reference modal parameters.

Note that the following notation is used:

$\mathcal{H}_j$ is the Hankel matrix in the reference state under parameter $P_j$

$\tilde{\mathcal{H}}_j$ is the consistent estimate of $\mathcal{H}_j$

$S_j$ is the left null space of $\mathcal{H}_j$

$\tilde{\mathcal{H}}(P)$ denotes the Hankel matrix estimate under parameter $P$ in the testing state

**Strategy 1: Hankel matrix interpolation** The outputs of the interpolated state-space model in (4.11) can be used to compute the respective reference Hankel matrix under parameter $P$, denoted by $\mathcal{H}(P)$, in order to obtain its left reference null space. Based on the properties of model (4.11), it is now shown that it is not necessary to compute the interpolated outputs explicitly, but that $\mathcal{H}(P)$ follows directly as an interpolated version of the reference Hankel matrices $\mathcal{H}_j, j = 1, \dots, u$, that can be estimated from measurements in the reference state at the respective reference parameters $P_j$.

From the output equation of model (4.11) it follows that $y_k = \sum_{j=1}^u \rho_j(P)y_{j,k}$. Hence, with (4.3) $\mathbf{E}(x_{j,k+1}y_k^T) = \rho_j(P)\mathbf{E}(x_{j,k+1}y_{j,k}^T) = \rho_j(P)G_j$ since the cross terms from measurements under different parameters and thus at different times are uncorrelated, and the cross-correlation between states and outputs of the interpolated model follows as $\underline{G}(P) = \mathbf{E}(x_{k+1}y_k^T) = [\rho_1(P)G_1^T \dots \rho_u(P)G_u^T]^T$. The output covariances of the interpolated model follow from (4.3) and (4.11) as

$$ \underline{R}_i(P) = \mathbf{E}(y_{k+i}y_k^T) = \underline{C}(P)\underline{A}^{i-1}\underline{G}(P) \\ = \sum_{j=1}^u \rho_j^2(P)C_jA_j^{i-1}G_j = \sum_{j=1}^u \rho_j^2(P)R_{j,i}, $$

The respective “interpolated” Hankel matrix is filled similarly as in (4.4) with these output covariances to

$$ \mathcal{H}(P) = \sum_{j=1}^{u} \rho_{j}^{2}(P) \mathcal{H}_{j}. \qquad (4.13) $$

The desired left reference null space $\underline{S}(P)$ can finally be computed as the left null space of the interpolated Hankel matrix $\underline{\mathcal{H}}(P)$.
**Strategy 2: Observability matrix interpolation** Suppose the reference data is recorded over a longer time period to obtain systems at different reference temperatures, changes in the excitation covariance must be suspected at the different reference points. In the context of civil structures such changes may result from different wind conditions or changes in the traffic loads affecting the output covariances from which the reference null space is derived. Due to general unknown excitation within the output-only framework those unmeasured changes are not covered in the system matrices of the state-space model and thus perturb the interpolation from strategy 1.

This motivates for computing the reference null space from the observability matrix in (2.62) exploiting the property $S_j^T O_j = 0$ (cf. (4.6)). A parametric version of the observability matrix uses the system matrices in (4.5) in the modal basis [DMH14] and can be obtained from the modal parameters of the system, which are not affected by the excitation covariance. A set of modal parameters under $P$ can be obtained by interpolating the eigenstructure under the different $P_j, j = 1, \dots, u$. Subsequently, the observability matrix $\mathcal{O}(P)$ and its left reference null space $\mathcal{S}(P)$ can be computed.

The local modal parameters are identified from the reference measurements at each of the parameters $P_j$, using e.g. subspace-based system identification as recalled in section 2.3. The eigenvalues $\lambda_{j,l}$, eigenvectors $\phi_{j,l}$ and mode shapes $\varphi_{j,l}$ of (4.2) yield in accordance to (2.15)

$$A_j \phi_{j,l} = \lambda_{j,l} \phi_{j,l}, \quad \varphi_{j,l} = C_j \phi_{j,l}, \quad l = 1, 2, \dots, n. \qquad (4.14)$$

The parameter vector $\theta_j^*$ can be defined as the collection of the modal parameters, or equivalently the eigenvalues and mode shapes of the system

$$\theta_j^* \stackrel{\text{def}}{=} \begin{bmatrix} \Lambda_j \\ \text{vec}(\Phi_j) \end{bmatrix} \qquad (4.15)$$

with $\Lambda_j = [\lambda_{j,1} \ \lambda_{j,2} \ \dots \ \lambda_{j,n}]^T$, $\Phi_j = [\varphi_{j,1} \ \varphi_{j,2} \ \dots \ \varphi_{j,n}]$. Notice that the mode shapes are only obtained up to a constant factor. Their comparability is required for the different reference parameters $P_j$ for interpolation. It is therefore supposed that the mode shapes are normalized such that one component is always one and that the chosen component is the same for the same mode under the different $j$.
The modal parameters at parameter P can then be interpolated as

$$ \underline{\theta}^*(P) = \sum_{j=1}^{u} \rho_j(P) \theta_j^* \qquad (4.16) $$

and used for the computation of the observability matrix $\underline{\mathcal{O}}(P)$. The desired left reference null space $\underline{\mathcal{S}}(P)$ can finally be computed as the left null space of $\underline{\mathcal{O}}(P)$.

**Strategy 3: Observability matrix with regression** Instead of interpolating the whole modal parameter vector (4.16) to obtain the modal parameters at P, they can also be analyzed individually in a regression, which allows a functional description of each modal parameter (eigenvalues and all mode shape components) in dependence of the external parameter. Let $\mathcal{G}$ denote the collection of functions generated by the regression analysis. The modal parameters at parameter P can then be obtained as

$$ \underline{\theta}^*(P) = \mathcal{G}(P) \qquad (4.17) $$

and used for the computation of the observability matrix $\underline{\mathcal{O}}(P)$, whose left reference null space is computed as $\underline{\mathcal{S}}(P)$.

### 4.2.3 Residual function and statistical evaluation

In the testing state at parameter P the Hankel matrix estimate $\hat{\mathcal{H}}(P)$ is computed from a data set $\mathcal{Y}_N(P)$ of length N. The Hankel matrix is confronted to the reference null space $\underline{\mathcal{S}}(P)$ derived from either of the three presented strategies. Replacing $\underline{\mathcal{S}}(P)$ by its consistent estimate $\hat{\mathcal{S}}(P)$, the subspace-based residual at parameter P in the testing state writes in accordance to (4.7) as

$$ \tilde{\zeta}(\underline{\theta}^*(P), \mathcal{Y}_N(P)) = \sqrt{N} \mathrm{vec}(\hat{\mathcal{S}}(P)^T \hat{\mathcal{H}}(P)). \qquad (4.18) $$

The test statistic analogous to (4.10) becomes

$$ \tilde{t} = \tilde{\zeta}(\underline{\theta}^*(P), \mathcal{Y}_N(P))^T \tilde{\Sigma}_{\zeta}(P)^{-1} \tilde{\zeta}(\underline{\theta}^*(P), \mathcal{Y}_N(P)), \qquad (4.19) $$

where $\tilde{\Sigma}_{\zeta}(P)$ is the covariance of the residual, which is detailed in section 4.2.4.