# Generative Surrogates for Fast Simulation: TPC Case

Fedor Ratnikov<EMAIL_ADDRESS>Artem Maevskiy
<EMAIL_ADDRESS>Alexander Zinchenko<EMAIL_ADDRESS>Victor
Riabov<EMAIL_ADDRESS>Alexey Sukhorosov<EMAIL_ADDRESS>Dmitrii
Evdokimov<EMAIL_ADDRESS>HSE University, 20 Myasnitskaya Ulitsa, Moscow,
Russia  Yandex School of Data Analysis, 11-2 Timura Frunze Street, Moscow,
Russia  Joint Institute for Nuclear Research, 6 Joliot-Curie St, Dubna, Moscow
Oblast, Russia  Petersburg Nuclear Physics Institute, 1, mkr. Orlova roshcha,
Gatchina, Leningradskaya Oblast, Russia

###### Abstract

Simulation of High Energy Physics experiments is widely used, necessary for
both detector and physics studies. Detailed Monte-Carlo simulation algorithms
are often limited due to the computational complexity of such methods, and
therefore faster approaches are desired. Generative Adversarial Networks
(GANs) are well suited for aggregating a number of detailed simulation steps
into a surrogate probability density estimator readily available for fast
sampling. In this work, we demonstrate the power of the GAN-based fast
simulation model on the use case of simulating the response for the Time
Projection Chamber (TPC) in the MPD experiment at the NICA accelerator
complex. We show that our model can generate high-fidelity TPC responses,
while accelerating the TPC simulation by at least an order of magnitude. We
describe alternative representation approaches for this problem and also
outline the roadmap for the deployment of our method into the software stack
of the experiment.

###### keywords:

fast simulation, time projection chamber, generative adversarial network

††journal: Nuclear Instruments and Methods in Physics Research A

##  1 Introduction

Efficient detector simulation algorithms are crucial for successful planning,
optimization and operation of High Energy Physics experiments [1]. The number
of samples that can be produced with the existing detailed Monte-Carlo
simulation approaches is often limited by the computational complexity of such
methods. Faster simulation can be achieved by aggregating a number of detailed
simulation steps into a single surrogate generative model that learns the
underlying probability density and allows for prompt sampling. Generative
Adversarial Networks (GANs) [2] is a deep learning method that is well-suited
for this goal [3]. In our previous work [4], we introduced a GAN-based fast
simulation model for the Time Projection Chamber (TPC) in the MPD experiment
at the NICA accelerator complex. Here we briefly review our fast simulation
approach along with its validation and then describe an alternative
representation approach for this problem. We also outline the deployment of
our method into the software stack of the experiment.

##  2 Fast simulation approach

Our fast simulation approach for the TPC detector consists of building a GAN
that generates the raw detector responses conditioned on the ground truth
event content [4]. For each heavy- ion collision, MPD TPC response contains
approximately $3\cdot 10^{7}$ numbers. To simplify the problem of generating
such a large amount of numbers, we factorize the response into an additive
decomposition of contributions from small track segments. Within our approach,
each track segment contributes to a geometrically localized window, making the
target response just an 8-by-16 matrix.

The built model is validated both at the low level with the raw detector
responses and at the high level using reconstructed track characteristics. The
low-level validation is done by calculating physically motivated projections
from the 8-by-16 response matrix: two barycenters (which are the proxies for
the reconstructed cluster coordinates), and the two widths and corresponding
covariance (which describe the shape of the signal and therefore define the
two-track resolution of the detector). The profiles for the distributions of
these quantities as a function of the conditional variables are then built and
compared between the model prediction and the detailed simulation, as shown in
Fig. 1 (left). According to these metrics, our model reproduces the detailed
simulation with high accuracy.

Pad barycenterPad coordinate fraction

Figure 1: (left) Low-level validation profiles. Thick lines correspond to (top
to bottom): $\mu+\sigma$, $\mu$, $\mu-\sigma$. Line thickness denotes the
statistical uncertainty. (right) Relative momentum resolution as a function of
the full momentum.

The high-level validation is done by integrating our model into the MPD
codebase and running the full reconstruction on the GAN-generated responses.
The reconstructed track characteristics are then compared with the detailed
simulation predictions. An example of such characteristic is shown in Fig. 1
(right). Most of the characteristics are reproduced with high fidelity. Some
minor discrepancies may be seen, which are well understood and expected to be
corrected for in the next iterations of the model development (see [4] for
more details). The TPC simulation with our model runs $\sim 12$ times faster
compared to the detailed simulation on the target hardware.

##  3 Alternative low-dimensional representation

Since the most important signal properties are the above-mentioned barycenters
and widths, one may train a GAN to generate these low-level characteristics
directly and then analytically and deterministically construct the desired pad
response matrices from them. In this way, we drastically reduce the
dimensionality of the problem and therefore may expect more stable and
efficient training, with quality at least as good as in our main approach.

$\widetilde{\chi^{2}}$ (lower better) Figure 2: Quality comparison between our
main model (the baseline) and the alternative model working in the low-
dimensional representation.

We have examined such low-dimensional approach with a choice of the
discretized 2D Gaussian function for constructing the response matrices.
Unfortunately, the discretization distorts the low-level characteristics, so,
in order to account for these distortions, we train the model in two stages.
At the first stage, which we call _pretraining_ , the model is trained to
reproduce the low-level characteristics directly. Then, at the second stage,
which we call _training_ , we introduce the distorting function in between the
generator and discriminator such that the generator learns to account for the
distortions automatically. The quality of our best such model, compared with
our main approach (the _baseline_), is shown in Fig. 2. Contrary to our
expectations, we find that this alternative model is less stable in training
and requires finer hyper-parameter tuning to achieve a reasonable quality
level.

##  4 Model deployment pipeline

The simulation configurations are expected to vary as MPD development
progresses and running conditions are adapted for. Each new configuration
requires our fast simulation model to be re-trained. It is therefore desired
to automate the training process as much as possible, including the following
steps: training data generation, model training, model evaluation and
selection. Since a multitude of models arise for different configurations, a
model library and database for storage and prompt model retrieval are
required. A preliminary scheme for such a pipeline is shown in Fig. 3.

Figure 3: The model deployment pipeline.

A prototype of this scheme has been implemented with the following framework
choices. We store the models in the ONNX111https://onnxruntime.ai/ format and
use ONNXruntime for inference. Airflow222https://airflow.apache.org/ is used
to manage the workflow with the calculation pipelines described as Directed
Acyclic Graphs (DAGs). We use the “Model Registry” component from
MLflow333https://mlflow.org/ as the model library, which allows us to use REST
API to download a model for inference directly in the MPD code.

##  5 Summary

In this work, we have demonstrated our GAN-based fast simulation model for the
Time Projection Chamber (TPC) in the MPD experiment at the NICA accelerator
complex. Our model can generate high-fidelity TPC responses with an order of
magnitude acceleration. We have also described an alternative low-dimensional
representation approach for this problem which, contrarily to our
expectations, requires more thorough hyper-parameter tuning to achieve a
reasonable quality level. Finally, we have shown the strategy for the
deployment of our method into the software stack of the experiment.

## References

  * [1] J. Albrecht, et al., A Roadmap for HEP Software and Computing R&D for the 2020s, Comput. Softw. Big Sci. 3 (1) (2019) 7.  arXiv:1712.06982, doi:10.1007/s41781-018-0018-8.
  * [2] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, Generative Adversarial Networks, Advances in Neural Information Processing Systems 27.  arXiv:1406.2661.
URL
https://proceedings.neurips.cc/paper/2014/hash/5ca3e9b122f61f8f06494c97b1afccf3-Abstract.html

  * [3] L. de Oliveira, M. Paganini, B. Nachman, Learning Particle Physics by Example: Location-Aware Generative Adversarial Networks for Physics Synthesis, Comput. Softw. Big Sci. 1 (1) (2017) 4.  arXiv:1701.05927, doi:10.1007/s41781-017-0004-6.
  * [4] A. Maevskiy, F. Ratnikov, A. Zinchenko, V. Riabov, Simulating the time projection chamber responses at the MPD detector using generative adversarial networks, Eur. Phys. J. C 81 (7) (2021) 599.  arXiv:2012.04595, doi:10.1140/epjc/s10052-021-09366-4.


# Modelling variability in vibration-based PBSHM via a generalised population
form

T.A. Dardeno<EMAIL_ADDRESS>L.A. Bull  R.S. Mills  N. Dervilis
K. Worden  Dynamics Research Group, Department of Mechanical Engineering,
University of Sheffield, Sheffield S1 3JD, UK  The Alan Turing Institute, The
British Library, London NW1 2DB, UK

###### Abstract

Structural health monitoring (SHM) has been an active research area for the
last three decades, and has accumulated a number of critical advances over
that period, as can be seen in the literature. However, SHM is still facing
challenges because of the paucity of damage-state data, operational and
environmental fluctuations, repeatability issues, and changes in boundary
conditions. These issues present as inconsistencies in the captured features
and can have a huge impact on the practical implementation, but more
critically, on the generalisation of the technology. Population-based SHM has
been designed to address some of these concerns by modelling and transferring
missing information using data collected from groups of similar structures.

In this work, vibration data were collected from four healthy, nominally-
identical, full-scale composite helicopter blades. Manufacturing differences
(e.g., slight differences in geometry and/or material properties), among the
blades presented as variability in their structural dynamics, which can be
very problematic for SHM based on machine learning from vibration data. This
work aims to address this variability by defining a general model for the
frequency response functions of the blades, called a _form_ , using mixtures
of Gaussian processes.

###### keywords:

repeatability , uncertainty , Gaussian processes , population-based SHM

††journal: Journal of Sound and Vibration

##  1 Introduction

Engineering structures are inherently uncertain, even when structures can be
considered _nominally-identical_ , because of manufacturing differences
causing variation in geometry and material properties. These differences, as
well as variations caused by ageing parts and changes in testing conditions
(e.g., operational and environmental fluctuations or changes in boundary
stiffness), make generalisation among structures difficult, particularly with
respect to their dynamic properties [1, 2]. The current work is focussed on
applying structural health monitoring (SHM), to a set of nominally-identical
(i.e., _homogeneous_ [3, 4]) structures for the purpose of damage detection.
(As opposed to a _heterogeneous_ population [5, 6], where the members are more
disparate, such as different designs of a suspension bridge). However, these
fluctuations unrelated to damage are challenging for SHM, as damage-sensitive
features may also be sensitive to normal variations, which may cause issues
with discerning between damaged and healthy states [7, 8, 9, 10]. For example,
it is well-known that structural damage is often associated with a reduction
in stiffness. Depending on the location and severity of the damage, and the
complexity of the examined mode, this damage may present as a decrease in
natural frequency. Increased ambient temperature, loosening of bolts at a
support, or similar situations may also reduce stiffness and can mimic damage.
Likewise, conditions that increase stiffness (e.g., decreased ambient
temperature) can increase natural frequency and mask damage.

This sensitivity of natural frequency to environmental/operational changes in
addition to damage has been demonstrated experimentally [8, 9, 10]. Alampalli
[8], performed impact modal testing on a small 6.76 x 5.26 metre bridge,
before and after damaging the bridge with a saw cut made across the bottom
flanges of both the main girders. They found that for temperatures above
freezing, the stiffness reduction caused by damage was reflected as a downward
shift in natural frequencies (-2% to -11%, depending on the mode), as
expected. However, they also found that for temperatures below freezing, the
measured frequencies increased significantly, primarily from freezing of
accumulated moisture within the supports (+25% to +58%, depending on the
mode), which masked the reduction caused by damage. In addition, uncertainty
among nominally-identical structures (e.g., slight variation in internal
structure and material properties), like those caused by environmental and
other fluctuations, may be erroneously flagged as damage or mask subtle damage
in a structure. Cawley [9], initiated a crack at the fixed end of a cantilever
beam and evaluated the natural frequencies as the beam length varied. He found
that the change in natural frequency resulting from a cut through 2% of the
beam depth was 40 times smaller than that caused by a 2% increase in the beam
length [9, 10]. Given that these small differences among structures and
testing conditions can have such a large impact on dynamics and damage
detection, one would ideally aim to collect new data for each structure under
each set of testing conditions. However, collecting comprehensive data
covering a range of operational and damaged conditions from each individual
structure is typically not feasible (a wind farm, for example, may contain
hundreds or even thousands of wind turbines). In such cases, population-based
SHM (PBSHM) can be used to make inferences between different members of a
_population_ , i.e., groups of similar systems [3, 4, 5, 6, 11].

PBSHM seeks to transfer valuable information across similar structures. In
certain cases, the behaviour of the group can be represented using a general
model. Bull _et al._ used Gaussian process (GP) regression of real [3, 4], and
imaginary [4], frequency response functions (FRFs) to develop a generic
representation, called a _form_ , of a population of nominally-identical, but
slightly different, eight degree-of-freedom (DOF) systems. The population was
comprised of multiple simulated healthy and damaged systems [3, 4], and one
experimental rig, which was used to test the form in [3]. Normal-condition
data from several members of the population were used to train the form, and
the posterior predictive distributions of the GP were used to aid novelty
detection across the population, using the Mahalanobis squared-distance
novelty index [12]. The form was found to predict accurately the healthy and
damaged states across the population [3, 4]. Lázaro-Gredilla _et al._ [13],
introduced the overlapping mixture of Gaussian processes (OMGP), which uses a
variational Bayesian approach to learn hyperparameters and labels by
clustering data according to various trajectories across the input space. This
technique was applied to the population form in [3, 11], using an OMGP to
infer multivalued wind-turbine power-curve data, with categorisation of the
data performed in an unsupervised manner. Novelty detection was performed
using the negative log-marginal-likelihood of the OMGP, and correctly
identified _good_ and _bad_ power curves for the majority of data tested [3,
11].

In this work, data were collected from four healthy, nominally-identical
helicopter blades. FRFs were computed from the measured time-domain data, and
these data were used to develop a form for the blades, first using a
supervised mixture of GPs, and second using an (unsupervised) OMGP. The
posterior predictive distributions from the OMGP were later applied to a
novelty detection approach to evaluate new (simulated) FRF data against the
form. As such, this work builds upon that performed in [3, 4], by extending
the concept of a population form to a set of nominally-identical helicopter
blades with a range of undamaged-condition FRF data, using mixtures of
probabilistic regression models [3, 11, 13], to account for the large
variability in the FRFs of the blades. The significant discrepancies in the
experimental data set are representative of the differences in the blades in
practice, with variability arising from small variations in the internal
structure and material properties of the helicopter blades, as well as slight
changes in boundary conditions. The small population of helicopter blades
provided a useful data set to develop and test the techniques presented
herein, such that these methods can be later applied to other homogeneous
populations such as a wind farm.

Although the specific application of this work relates to population-based
damage detection, this work is also applicable to the more general problem of
dynamic analysis of a population. Attempts to quantify uncertainty in
vibration data have included parametric [14], and nonparametric [2, 15],
approaches. Parametric models allow incorporation of engineering knowledge,
but can suffer from over-training (if highly complex), or over-simplification.
Nonparametric models, such as the approach used in [2], and the polynominal-
chaos expansion used in [15], in general do not incorporate prior engineering
knowledge, and this can lead to problematic solutions because of the high
model flexibility. Gaussian processes provide a good middle ground, as they
are nonparametric, but allow for the use of a parametric mean function where
prior knowledge about the data can be incorporated into the model. The (OM)GP
form provides a probabilistic model of the FRF that can capture uncertainty
caused by the variations described above as well as measurement noise.

This work formalises the use of a modal analysis approach to represent a group
of similar _dynamic_ systems, using a data-driven model (or representation)
referred to as the _population form_ [3, 4]. Unlike previous work [4], the
form here uses a mixture of regressions to represent _population
dynamics/FRFs_ \- rather than a single, averaged model (complete pooling). In
turn, the modelling procedure can approximate population data with more
complex variations and different modes in the posterior distribution; for
example, relating to distinct operating conditions, characteristic subgroups
in the population, or individual members. While mixture models have been used
to represent the _form_ in previous work, relating to wind farm power
prediction [3, 11], these models were not physics informed [16]; here, in
contrast, the proposed (semi-parametric) form is constructed around the
established theory of linear modal analysis – allowing for the natural
inclusion of domain (engineering) expertise via interpretable parameters
(e.g., natural frequencies, damping ratios, mode shapes).

The layout of this paper is as follows. Section 2 provides an overview of the
theory used in this work, including modal analysis, conventional GPs, and
OMGP. Section 3 describes the experimental campaign performed on the
population of nominally-identical helicopter blades and introduces the data
set used in the analyses. Section 4 discusses how a generalised normal
condition was determined for the helicopter blades, first using a supervised
mixture of GPs, and second using an (unsupervised) OMGP. Finally, Section 5
demonstrates how the OMGP can be used to inform damage detection, by computing
the marginal likelihood for a series of (simulated) damaged FRFs.

##  2 Background theory

In this work, a generalised normal condition was determined for a small
population of nominally-identical helicopter blades. A form was developed for
the population, first using a supervised mixture of GPs, and second using an
(unsupervised) OMGP. A modal-based FRF estimation technique was parameterised
and used to compute a non-zero prior mean function for the GPs. The marginal
likelihood (evidence) was then used along with the posterior predictive
distributions of the OMGP to evaluate (simulated) new data for novelty. A
theoretical background for these methodologies is provided below.

###  2.1 Modal analysis

The equation of motion for a multiple DOF system can be written as,

$\mathbf{M}\ddot{\mathbf{u}}(t)+\mathbf{C}\dot{\mathbf{u}}(t)+\mathbf{K}{\mathbf{u}}(t)=\mathbf{z}(t)$
(1)

where $\ddot{\mathbf{u}}(t)$, $\dot{\mathbf{u}}(t)$, ${\mathbf{u}}(t)$, and
$\mathbf{z}(t)$ are acceleration, velocity, displacement, and force,
respectively. In most cases, the mass $\mathbf{M}$, damping $\mathbf{C}$, and
stiffness $\mathbf{K}$ matrices are coupled. For linear systems, and in the
absence of viscous damping, the equation of motion can be decoupled, such that
the system is represented by multiple single degree-of-freedom (SDOF)
oscillators. This decoupling is performed via the eigenvalue expression,

$[\mathbf{K}-\omega_{n}^{2}\mathbf{M}]\boldsymbol{\Psi}=0$ (2)

and yields the natural frequencies in radians, $\omega_{n}$, and mode shapes,
$\boldsymbol{\Psi}$, of the system.

The physical equation of motion can then be cast in modal space to give the
uncoupled modal equations, with modal coordinates ${\mathbf{p}}(t)$, written
as,

$\boldsymbol{\Psi}^{\text{T}}\mathbf{M}\boldsymbol{\Psi}\ddot{\mathbf{p}}(t)+\boldsymbol{\Psi}^{\text{T}}\mathbf{C}\boldsymbol{\Psi}\dot{\mathbf{p}}(t)+\boldsymbol{\Psi}^{\text{T}}\mathbf{K}\boldsymbol{\Psi}{\mathbf{p}}(t)=\boldsymbol{\Psi}^{\text{T}}\mathbf{z}(t)$
(3)

where,

$\mathbf{u}(t)=\boldsymbol{\Psi}{\mathbf{p}}(t)$ (4)

Frequency response functions (FRFs) can be used to visualise the natural
frequency components of a system, and are computed by normalising the response
signal at a given location to the excitation force. This work used the
$\text{H}_{1}$ estimator to compute FRFs. For a response at location $i$
resulting from excitation at $j$, the $\text{H}_{1}$ estimator is computed as,

$\mathbf{H}_{ij}(\omega)=\frac{\mathbf{G}_{zu}(\omega)}{\mathbf{G}_{zz}(\omega)}$
(5)

where,

$\begin{split}&\mathbf{G}_{zu}(\omega)\overset{\Delta}{=}{\it\mathbb{E}}[\mathbf{S}_{u_{i}}(\omega)\mathbf{S}_{z_{j}}^{*}(\omega)]\\\
&\mathbf{G}_{zz}(\omega)\overset{\Delta}{=}{\it\mathbb{E}}[\mathbf{S}_{z_{j}}(\omega)\mathbf{S}_{z_{j}}^{*}(\omega)]\end{split}$
(6)

and,

$\begin{split}&\mathbf{S}_{z_{j}}(\omega)\overset{\Delta}{=}\mathcal{F}[{\mathbf{z}_{j}}(t)]\\\
&\mathbf{S}_{u_{i}}(\omega)\overset{\Delta}{=}\mathcal{F}[{\mathbf{u}_{i}}(t)]\end{split}$
(7)

The asterisk $*$ denotes complex conjugation, $\omega$ is frequency in
radians, and $\mathcal{F}$ is a discrete Fourier transform (this work used a
fast Fourier transform or FFT). The input force in the time domain at location
$j$ is ${\mathbf{z}_{j}}(t)$, and the output response (i.e., acceleration,
velocity, or displacement) in the time domain at $i$ is ${\mathbf{u}_{i}}(t)$.

With assumed linear behaviour and proportional damping, the accelerance FRF
(i.e., given acceleration response data) can also be estimated using modal
parameters,

$\mathbf{H}_{ij}(\omega)=-\omega^{2}\sum_{k=1}^{n}\frac{\text{A}_{ij}^{(k)}}{\omega_{nk}^{2}-\omega^{2}+2i\zeta_{k}\omega\omega_{nk}}$
(8)

where $\text{A}_{ij}^{(k)}$ is the residue for mode $k$, defined as the
product of the mass-normalised mode shapes at locations $i$ and $j$
($\text{A}_{ij}^{(k)}=\psi_{ik}\psi_{jk}$) [17]. The natural frequency
associated with mode $k$ is $\omega_{nk}$, and the modal damping associated
with mode $k$ is $\zeta_{k}$ [17]. In the examples below, Eq. (8) was used to
compute a non-zero prior mean function for the (OM)GP form.

###  2.2 Conventional Gaussian process regression

Gaussian processes provide a Bayesian nonparametric approach to
regression/classification and are a suitable basis for a population form
because of their robust uncertainty quantification. For a set of training data
$\mathcal{D}=\\{x_{i},y_{i}\\}_{i=1}^{N}=\\{{\mathbf{x},\mathbf{y}}\\}$, where
$x_{i}$ and $y_{i}$ are sets of $N$ inputs and outputs, respectively, the GP
regression attempts to obtain the predictive distribution for observation
$y_{*}$ given a new input $x_{*}$ and training data $\mathcal{D}$ [13]. The GP
regression model can be expressed as a noiseless latent function $f(x_{i})$
plus an independent noise term $\epsilon_{i}$ [13], written as,

$y_{i}=f(x_{i})+\epsilon_{i}\ $ (9)

A GP prior is placed over the latent function $f(x_{i})$ and a Gaussian prior
is placed over the noise term $\epsilon_{i}$,

$f(x_{i})\sim\mathcal{GP}(m(x_{i}),k(x_{i},x_{j})),\;\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2})\
$ (10)

where $m(x_{i})$ is a mean function, $k(x_{i},x_{j})$ is a covariance
function, and $\sigma^{2}$ is a noise variance term and the first
hyperparameter. Over a finite, abitrary set of inputs
$\mathbf{x}=\\{x_{1},...,x_{N}\\}$, the GP is a (joint) multivariate Gaussian
[11, 18],

$p\left(\mathbf{f}\;|\;\mathbf{x}\right)=\mathcal{N}\left(\mathbf{m},\mathbf{K_{xx}}\right)$
(11)

where $\mathbf{m}=\\{m(x_{i}),...,m(x_{N})\\}$,
$\mathbf{f}=\\{f(x_{i}),...,f(x_{N})\\}$, and $\mathbf{K_{xx}}$ is the
covariance matrix, such that $\mathbf{K_{xx}}[i,j]=k(x_{i},x_{j})\;\forall
i,j\in\\{1,...,N\\}$ [11]. Note that in certain instances within this paper,
square brackets are used to index matrices and vectors for readability.

The mean function $m(x_{i})$ is often assumed to be zero, given that the GP is
flexible enough to model arbitrary trends [18, 19, 11]. However, the GP fit
may be improved, by incorporating _a priori_ knowledge from engineering
judgement or an understanding of the system physics, into the GP prior by
choosing suitable mean and covariance functions [3, 11]. Coupling behaviour
between $y_{i}$ and $y_{j}$, and therefore important properties such as
process variance and smoothness, is specified by the covariance function [13].
Appropriate covariance functions are selected according to the shape and
smoothness of the data. This work applied the commonly-used squared
exponential,

$k(x_{i},x_{j})=\sigma_{f}^{2}\exp\left(-\frac{1}{2l^{2}}(x_{i}-x_{j})^{2}\right)\
$ (12)

where $\sigma_{f}^{2}$ and $l$ are second and third hyperparameters,
representing the process variance and length scale, respectively. The length
scale $l$ influences how fast the correlation between outputs decays according
to the input separation, and the process variance $\sigma_{f}^{2}$ determines
the variance of the modelled signal [13].

The remaining hyperparameters $\alpha$ are those required by the mean
function, if non-zero, giving the full model hyperparameters as
$\boldsymbol{\theta}=\\{l,\sigma_{f},\sigma,\boldsymbol{\alpha}\\}$. In this
work, the modal damping-based FRF estimation from Eq. (8) provided the mean
function for the GP regression. Modal damping ($\zeta_{k}$), residue
($\text{A}_{ij}^{(k)}$), and natural frequency ($\omega_{nk}$), provided the
mean-function hyperparameters. Keeping $\boldsymbol{\theta}$ fixed, the joint
distribution of the training data
$\mathcal{D}=\\{{\mathbf{x},\mathbf{y}}\\}=\\{x_{i},y_{i}\\}_{i=1}^{N}$ and
unknown function values at test locations under the prior
$\\{{\mathbf{x_{*}},\mathbf{y_{*}}}\\}=\\{\mathbf{x_{*}}[i],\mathbf{y_{*}}[i]\\}_{i=1}^{M}$
is a multivariate Gaussian, written as,

$\displaystyle\begin{bmatrix}[l]\mathbf{y}\\\
\mathbf{y_{*}}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}[l]\mathbf{m}\\\
\mathbf{m_{*}}\end{bmatrix},\begin{bmatrix}[c]\mathbf{K_{xx}}+\mathbf{R}&\mathbf{K_{xx_{*}}}\\\
\mathbf{K_{x_{*}x}}&\mathbf{K_{x_{*}x_{*}}}+\mathbf{R_{*}}\end{bmatrix}\right)$
(13)

$\begin{split}&\mathbf{R}\overset{\Delta}{=}\sigma^{2}\mathbf{I}_{N}\\\
&\mathbf{R_{*}}\overset{\Delta}{=}\sigma^{2}\mathbf{I}_{M}\end{split}$ (14)

where $\mathbf{I}$ is an identity matrix of dimension $N$ or $M$ as specified
[11, 19]. Likewise,
$\mathbf{m_{*}}=\\{m\left(\mathbf{x_{*}}[i]\right)\\}^{M}_{i=1}$ is the mean
vector for new observations [11].

Conditioning the joint distribution in Eq. (13) on the training data in
$\mathcal{D}$ yields the predictive distribution over $\mathbf{y_{*}}$, with
mean $\boldsymbol{\mu_{*}}$ and covariance $\boldsymbol{\varSigma_{*}}$ [19],
defined as,

$\begin{split}&p(\mathbf{y_{*}}\;|\;\mathbf{x_{*}},\mathcal{D})=\mathcal{N}(\boldsymbol{\mu_{*}},\boldsymbol{\varSigma_{*}})\\\
&\boldsymbol{\mu_{*}}\overset{\Delta}{=}\mathbf{m_{*}}+\mathbf{K_{x_{*}x}}\left(\mathbf{K_{xx}}+\mathbf{R}\right)^{-1}\left(\mathbf{y}-\mathbf{m}\right)\\\
&\boldsymbol{\varSigma_{*}}\overset{\Delta}{=}\mathbf{K_{x_{*}x_{*}}}-\mathbf{K_{x_{*}x}}\left(\mathbf{K_{xx}}+\mathbf{R}\right)^{-1}\mathbf{K_{xx_{*}}}+\mathbf{R_{*}}\end{split}$
(15)

This work applied a Type-II maximum likelihood approach to learn the
hyperparameters, which maximises the marginal likelihood of the model
$p(\mathbf{y}\;|\;\mathbf{x},\boldsymbol{\theta})$ [19]. This procedure
involves minimising the negative log-marginal-likelihood to improve numerical
stability [18]. Therefore, the hyperparameters were optimised via,

$\hat{\theta}=\underset{\theta}{\operatorname{argmin}}\\{-\log
p(\mathbf{y}\;|\;\mathbf{x},\boldsymbol{\theta})\\}$ (16)

where,

$\begin{split}&-\log
p(\mathbf{y}\;|\;\mathbf{x},\boldsymbol{\theta})=-\log\mathcal{N}(\mathbf{y}\;|\;\mathbf{m},\;\mathbf{K_{xx}}+\mathbf{R})\\\
&=\frac{N}{2}\log(2\pi)+\frac{1}{2}\log|\mathbf{K_{xx}}+\mathbf{R}|\\\
&+\frac{1}{2}[(\mathbf{y}-\mathbf{m})^{\textsf{T}}(\mathbf{K_{xx}}+\mathbf{R})^{-1}(\mathbf{y}-\mathbf{m})]\end{split}$
(17)

Constrained nonlinear optimisation was performed via the _fmincon_ function in
MATLAB, and appropriate bounds were applied to the mean-function
hyperparameters, given prior knowledge of the modal characteristics of the
helicopter blades.

###  2.3 An overlapping mixture of Gaussian processes (OMGP)

In one example, this work used the conventional GP equations to evaluate a
supervised mixture of GPs by assuming known which data belonged to each
helicopter blade, and fitting each data set separately. In a second example,
the data were fitted using an overlapping mixture of Gaussian processes
(OMGP), where the number of trajectory functions were assumed known _a priori_
, but it was not assumed which data belonged to a given class. In this case,
the unknown class was the particular helicopter blade that produced the
measurements. In practice, it is likely that this information would be known;
however, there are regularly other forms of missing data labels, such as
changes in operating conditions, or possible damage. The helicopter blades
case study is also representative of these types of situations, commonly
encountered in PBSHM, where certain information is unavailable. As it is an
unsupervised learning algorithm, the OMGP is a useful and flexible method for
developing the population form in the absence of complete information.

The OMGP [13, 20, 3, 11], describes a data set by evaluating each observation
using one of _K_ latent functions plus additive noise,

$y_{i}^{(k)}=\\{f^{(k)}(x_{i})+\epsilon_{i}\\}_{k=1}^{K}$ (18)

The OMGP is unsupervised, i.e., the labels that assign the data to a given
function are unknown. Therefore, another latent variable is introduced,
$\mathbf{Z}$, which is a binary indicator matrix. The matrix $\mathbf{Z}$,
with entries $\mathbf{Z}[i,k]$, is entirely populated with zeroes, except for
one non-zero entry per row. The non-zero entries, $\mathbf{Z}[i,k]\neq 0$,
indicate that observation $i$ was generated by function $k$. The likelihood of
the OMGP is therefore written as [13, 11],

$p(\mathbf{y}\;|\;\\{\mathbf{f}^{(k)}\\}_{k=1}^{K},\mathbf{Z},\mathbf{x})=\prod_{i,k=1}^{N,K}p{(y_{i}\;|\;{f}^{(k)}(x_{i}))}^{\mathbf{Z}[i,k]}$
(19)

Prior distributions are then placed over the latent functions and variables,
such that [11],

$P(\mathbf{Z})=\prod_{i,k=1}^{N,K}{\boldsymbol{\varPi}[i,k]}^{\mathbf{Z}[i,k]}$
(20)

$f^{(k)}(x_{i})\sim\mathcal{GP}(m^{(k)}(x_{i}),k^{(k)}(x_{i},x_{j})),\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$
(21)

where $P(\mathbf{Z})$ is the prior over the indicator matrix $\mathbf{Z}$, and
$\boldsymbol{\varPi}[i,:]$ is a histogram over the $K$ components for the
$i$th observation, and $\sum_{k=1}^{K}\boldsymbol{\varPi}[i,k]=1$ [11]. The
terms in Eq. (21) are independent GP priors over each latent function
$f^{(k)}$ with distinct mean $m^{(k)}(x_{i})$ and kernel
$k^{(k)}(x_{i},x_{j})$ functions [11]. A shared hyperparameter $\sigma$, is
used to define the prior over the noise variances, to reduce the number of
latent variables in the calculations [11]. As with the conventional GP, this
work used the modal damping-based FRF estimation from Eq. (8) as the mean
function for the GP regression. Modal damping, residue, and natural frequency
provided the mean-function hyperparameters.

####  2.3.1 Variational inference

Exact computation of the posterior distribution
$p(\\{\mathbf{f}^{(k)}\\}_{k=1}^{K},\mathbf{Z}\;|\;\mathbf{x},\mathbf{y})$ is
intractable; therefore, a variational inference [21] and Expectation
Maximisation (EM) scheme was employed to estimate the posterior, by specifying
an approximate density family $q(\mathbf{a})\in\mathcal{Q}$ over the target
conditional $p(\mathbf{a}\;|\;\mathbf{b})$ [21, 13, 3, 11]. The best
approximation $\hat{q}(\mathbf{a})$ of the true posterior distribution
$p(\mathbf{a}\;|\;\mathbf{b})$ can be expressed via the KL-divergence [21,
11],

$\hat{q}(\mathbf{a})=\underset{q(\mathbf{a})\in\mathcal{Q}}{\operatorname{argmin}}\\{\text{KL}(q(\mathbf{a})\;||\;p(\mathbf{a}\;|\;\mathbf{b}))\\}$
(22)

$\begin{split}\text{KL}(q(\mathbf{a})\;||\;p(\mathbf{a}\;|\;\mathbf{b}))&=\mathbb{E}_{q(\mathbf{a})}[\text{log
}q(\mathbf{a})]-\mathbb{E}_{q(\mathbf{a})}[\text{log
}p(\mathbf{a}\;|\;\mathbf{b})]\\\ &=\mathbb{E}_{q(\mathbf{a})}[\text{log
}q(\mathbf{a})]-\mathbb{E}_{q(\mathbf{a})}[\text{log
}p(\mathbf{a},\mathbf{b})]\\\ &+\text{log }p(\mathbf{b})\end{split}$ (23)

where, for this work,
$\mathbf{a}\overset{\Delta}{=}\\{\\{\mathbf{f}^{(k)}\\},\mathbf{Z}\\}$ and
$\mathbf{b}\overset{\Delta}{=}\\{\mathbf{y}\\}$. However, $p(\mathbf{b})$ is
intractable [21, 11]; therefore, rather than the KL-divergence, an alternative
object is defined, called the _evidence lower bound_ or _elbo_. The elbo,
$\mathcal{L}_{b}$, is equivalent to the negative KL-divergence in Eq. (23), up
to the term log $p(\mathbf{b})$, which is a constant with respect to
$q(\mathbf{z})$ [11],

$\begin{split}\mathcal{L}_{b}(\mathbf{a})&=\mathbb{E}_{q(\mathbf{a})}[\text{log
}p(\mathbf{a},\mathbf{b})]-\mathbb{E}_{q(\mathbf{a})}[\text{log
}q(\mathbf{a})]\\\ &=\int
q(\mathbf{a})\text{log}\frac{p(\mathbf{a},\mathbf{b})}{q(\mathbf{a})}\text{d}\mathbf{a}\end{split}$
(24)

Maximising $\mathcal{L}_{b}$ minimises the KL-divergence between
$q(\mathbf{a})$ and $p(\mathbf{a}\;|\;\mathbf{b})$. Rearranging Eq. (23) and
substituting in Eq. (24) yields,

$\text{log
}p(\mathbf{b})=\text{KL}(q(\mathbf{a})\;||\;p(\mathbf{a}\;|\;\mathbf{b}))+\mathcal{L}_{b}$
(25)

which shows that because $\text{KL}(\cdot)\geq 0$ [22], the evidence is lower-
bounded by the elbo $\mathcal{L}_{b}$, i.e., log
$p(\mathbf{b})\geq\mathcal{L}_{b}$ [11]. Substituting
$\mathbf{a}\overset{\Delta}{=}\\{\\{\mathbf{f}^{(k)}\\},\mathbf{Z}\\}$ and
$\mathbf{b}\overset{\Delta}{=}\\{\mathbf{y}\\}$ into Eq. (24) yields,

$\begin{split}&\text{ log }p(\mathbf{y}\;|\;\mathbf{x})=\text{ log}\iint
p(\\{\mathbf{f}^{(k)}\\},\mathbf{Z},\mathbf{y},\mathbf{x})p(\\{\mathbf{f}^{(k)}\\}),p(\mathbf{Z})\text{d}\\{\mathbf{f}^{(k)}\\}\text{d}\mathbf{Z}\\\
&\geq\mathcal{L}_{b}=\iint q(\\{\mathbf{f}^{(k)}\\},\mathbf{Z})\text{ log
}\frac{p(\\{\mathbf{f}^{(k)}\\},\mathbf{Z},\mathbf{y},\mathbf{x})}{q(\\{\mathbf{f}^{(k)}\\},\mathbf{Z})}\text{d}\\{\mathbf{f}^{(k)}\\}\text{d}\mathbf{Z}\\\
&=\iint q(\\{\mathbf{f}^{(k)}\\},\mathbf{Z})\text{ log
}\frac{p(\mathbf{y}\;|\;\\{\mathbf{f}^{(k)}\\},\mathbf{Z},\mathbf{x})p(\mathbf{Z})\prod_{k=1}^{K}p(\mathbf{f}^{(k)}\;|\;\mathbf{x})}{q(\\{\mathbf{f}^{(k)}\\},\mathbf{Z})}\text{d}\\{\mathbf{f}^{(k)}\\}\text{d}\mathbf{Z}\\\
\end{split}$ (26)

An improved bound, $\mathcal{L}_{bc}$, that lower-bounds the marginal
likelihood while also being an upper-bound on the standard variational bound
$\mathcal{L}_{b}$ [23, 11], was implemented in this work. The improved lower
bound can be expressed as [24, 11],

$\begin{split}&\text{ log
}p(\mathbf{y}\;|\;\mathbf{x})\geq\boldsymbol{\mathcal{L}}_{bc}\\\
&=\sum_{k=1}^{K}\left(\frac{1}{2}\;\Big{|}\\!\Big{|}\;\mathbf{R}^{(k)\textsf{T}}\backslash\left(\mathbf{B}^{(k)\frac{1}{2}}\left(\mathbf{y}-\mathbf{m}^{(k)}\right)\right)\;\Big{|}\\!\Big{|}\,\rule{-2.0pt}{15.0pt}^{2}-\sum_{i=1}^{N}\text{log
}\mathbf{R}^{(k)}[i,i]\right)\;...\;\\\
&-\text{KL}\left(q(\mathbf{Z})\;\Big{|}\\!\Big{|}\;p(\mathbf{Z})\right)-\frac{1}{2}\sum_{i=1,k=1}^{N,K}\boldsymbol{\hat{\varPi}}[i,k]\text{
log}\left(2\pi\sigma^{2}\right)\\\
&\mathbf{R}^{(k)}\overset{\Delta}{=}\text{chol}\left(\mathbf{I}+\mathbf{B}^{(k)\frac{1}{2}}\mathbf{K}_{\mathbf{xx}}^{(k)}\mathbf{B}^{(k)\frac{1}{2}}\right)\\\
&\text{KL}\left(q(\mathbf{Z})\;\Big{|}\\!\Big{|}\;p(\mathbf{Z})\right)\overset{\Delta}{=}\sum_{i=1,k=1}^{N,K}\boldsymbol{\hat{\varPi}}[i,k]\text{
log
}\frac{\boldsymbol{\varPi}[i,k]}{\boldsymbol{\hat{\varPi}}[i,k]}\end{split}$
(27)

where the backslash operator solves the system of linear equations
$\mathbf{Ax}=\mathbf{B}$ for $\mathbf{x}$ via $\mathbf{A}\backslash\mathbf{B}$
and $\text{chol}(\cdot)$ is the Cholesky decomposition. The quantity
$\mathbf{B}^{(k)}$ is an $N\times N$ diagonal matrix with elements,

$\mathbf{B}^{(k)}=\text{diag}\left(\left\\{\frac{\boldsymbol{\hat{\varPi}}[1,k]}{\sigma^{2}}\;,\;...\;,\;\frac{\boldsymbol{\hat{\varPi}}[N,k]}{\sigma^{2}}\right\\}\right)$
(28)

The improved lower bound $\mathcal{L}_{bc}$ depends only on
$p\left(\mathbf{Z}\right)$ and is independent of
$p\left(\\{\mathbf{f}^{(k)}\\}\right)$; as such, $\mathcal{L}_{bc}$ can be
referred to as the _marginalised variational bound_ [13, 11]. For a given
$p\left(\mathbf{Z}\right)$, when an optimal choice for
$p\left(\\{\mathbf{f}^{(k)}\\}\right)$ is made, $\mathcal{L}_{bc}$ is
equivalent to $\mathcal{L}_{b}$ [13, 11]. This idea means that the bound is
more stable over a variety of hyperparameters, because the improved bound is
independent of $p\left(\\{\mathbf{f}^{(k)}\\}\right)$ [13, 23, 11]. Using the
improved lower bound is therefore useful in the M-step of the Expectation
Maximisation (EM) scheme, when $\mathcal{L}_{bc}$ is optimised with respect to
the hyperparameters until convergence, with $p\left(\mathbf{Z}\right)$ held
fixed [13].

####  2.3.2 E-step: Mean-field updates

A mean-field assumption was employed, and a family of variational
distributions $q(\mathbf{a})\in\mathcal{Q}$ were defined such that each latent
variable could be optimised independently, i.e.,
$q(\\{\mathbf{f}^{(k)}\\},\mathbf{Z})=q(\\{\mathbf{f}^{(k)}\\})q(\mathbf{Z})$.
This assumption allowed for the optimisation of $\mathcal{L}_{b}$ (or
$\mathcal{L}_{bc}$) with respect to a given latent variable while keeping the
others fixed, and iteratively updating each variable until convergence of the
lower bound was reached (i.e., EM Maximisation) [11].

If $q(\\{\mathbf{f}^{(k)}\\})$ and the marginals for each component
$q(\mathbf{f}^{(k)})=\mathcal{N}(\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})$
are assumed known, $\mathcal{L}_{b}$ can be analytically maximised with
respect to $q(\mathbf{Z})$ by constraining $q$ to be a probability density and
setting the derivative of the bound to zero [13, 11],

$\begin{split}&q(\boldsymbol{Z})=\prod_{i,k=1}^{N,K}{\boldsymbol{\hat{\varPi}}[i,k]}^{\boldsymbol{Z}[i,k]}\text{,
s.t. }\boldsymbol{\hat{\varPi}}[i,k]\propto\boldsymbol{\varPi}[i,k]\text{
exp}(a_{ik})\\\
&a_{ik}\overset{\Delta}{=}-\frac{1}{2\sigma^{2}}\left(\left(y_{i}-\boldsymbol{\mu}_{i}^{(k)}\right)^{2}+\boldsymbol{\varSigma}^{(k)}[i,i]\right)-\frac{1}{2}\text{
log}\left(2\pi\sigma^{2}\right)\end{split}$ (29)

where Eq. (29) implies that $q(\mathbf{Z})$ is factorised for each sample [13,
11].

Now, assuming that $q(\mathbf{Z})$ is known, the lower bound can be
analytically maximised with respect to $q(\\{\mathbf{f}^{(k)}\\})$, as in,

$\begin{split}&q(\mathbf{f}^{(k)})=\mathcal{N}\left(\mathbf{f}^{(k)}\;|\;\boldsymbol{\mu}^{(k)},\boldsymbol{\varSigma}^{(k)}\right)\\\
&\boldsymbol{\varSigma}^{(k)}\overset{\Delta}{=}\left(\mathbf{K}_{xx}^{-1(k)}+\mathbf{B}^{(k)}\right)^{-1}\\\
&\boldsymbol{\mu}^{(k)}\overset{\Delta}{=}\mathbf{m}^{(k)}+\boldsymbol{\varSigma}^{(k)}\mathbf{B}^{(k)}\left(\mathbf{y}-\mathbf{m}^{(k)}\right)\\\
\end{split}$ (30)

The best candidate $\hat{q}(\mathbf{a})$ that most closely approximates the
true posterior distribution was found by initialising $q(\mathbf{Z})$ and
$q\left(\\{\mathbf{f}^{(k)}\\}\right)$ from their respective priors, and
iteratively updating them by alternating Eqs. (29) and (30). As both updates
are optimal with respect to the distribution that they compute, they are
guaranteed to increase the lower bound on the log-marginal-likelihood [13,
11], computed via Eq. (LABEL:eq:Lb3).

####  2.3.3 M-step: Hyperparameter optimisation

In the E-step, $\mathcal{L}_{bc}$ was optimised with respect to
$p\left(\mathbf{Z}\right)$ until convergence by iterating Eqs. (29) and (30),
with the hyperparameters held fixed. In the M-step, the lower bound is
optimised with respect to all model hyperparameters, as in,

$\Big{\\{}\\{\hat{\theta}_{k}\\}_{k=1}^{K},\hat{\boldsymbol{\varPi}}\Big{\\}}=\underset{\big{\\{}\\{\theta_{k}\\}_{k=1}^{K},\boldsymbol{\varPi}\big{\\}}}{\operatorname{argmax}}\\{\mathcal{L}_{bc}\\}$
(31)

while keeping the approximate posterior $q\left(\mathbf{Z}\right)$ fixed. The
E- and M-steps are then alternated until the lower bound converges. In this
work, constrained nonlinear optimisation was performed via the _fmincon_
function in MATLAB, and appropriate bounds were applied to the mean-function
hyperparameters given prior knowledge of the modal characteristics of the
helicopter blades.

####  2.3.4 Predictive equations

Once learnt, the OMGP fit (or form) can be used to estimate the latent
variables and functions and to inform damage detection. For training data
$\mathcal{D}$, the _maximum a posteriori_ (MAP) estimate can be used to
categorise observations according to the most likely component $k$, via [11],

$\hat{k}_{i}=\underset{k}{\operatorname{argmax}}\\{\hat{\boldsymbol{\Pi}}[i,k]\\}$
(32)

The MAP class component for new data $\hat{k}_{*}$ is then defined as [11],

$\hat{k}_{*}=\underset{k_{*}}{\operatorname{argmax}}\\{p({k}_{*}\;|\;\mathbb{x}_{*},\mathbb{y}_{*},\mathcal{D})\\}$
(33)

where, for a set of test data, the posterior predictive class component
$k_{*}$ is expressed using Bayes Rule (note that to classify new data
according to $k_{*}$, both $\mathbf{x_{*}}$ and $\mathbf{y_{*}}$ must be
observed [11]),

$p({k}_{*}\;|\;\mathbb{x}_{*},\mathbb{y}_{*},\mathcal{D})=\frac{{p(\mathbb{y}_{*}\;|\;\mathbb{x}_{*},k,\mathcal{D})p(k_{*})}}{{p(\mathbb{y}_{*}\;|\;\mathbb{x}_{*},\mathcal{D})}}$
(34)

The numerator of Eq. (34), or unnormalised posterior, is defined as [11],

$p(\mathbb{y}_{*}\;|\;\mathbb{x}_{*},k,\mathcal{D})p(k_{*})\overset{\Delta}{=}\mathcal{N}(\mathbb{y}_{*}\;|\;\boldsymbol{\mu}_{*}^{(k)},\boldsymbol{\varSigma}_{*}^{(k)})\boldsymbol{\varPi}[*,k]$
(35)

and the denominator of Eq. (34), which is the marginal likelihood (evidence),
is [11],

$\begin{split}p(\mathbb{y}_{*}\;|\;\mathbb{x}_{*},\mathcal{D})&\approx\sum_{k=1}^{K}\boldsymbol{\varPi}[*,k]\int
p\left(\mathbb{y}_{*}\;|\;\mathbf{f}^{(k)},\mathbb{x}_{*},\mathcal{D}\right)q\left(\mathbf{f}^{(k)}\;|\;\mathcal{D}\right)\text{d}\mathbf{f}^{(k)}\\\
&=\sum_{k=1}^{K}\mathcal{N}\left(\mathbb{y}_{*}\;|\;\boldsymbol{\mu}_{*}^{(k)},\boldsymbol{\varSigma}_{*}^{(k)}\right)\boldsymbol{\varPi}[*,k]\\\
\end{split}$ (36)

where [11],

$\begin{split}&\boldsymbol{\mu}_{*}^{(k)}\overset{\Delta}{=}\mathbf{m}_{*}^{(k)}+\mathbf{K}^{(k)}_{\mathbf{{x_{*}x}}}\left(\mathbf{K}^{(k)}_{\mathbf{{xx}}}+\mathbf{B}^{(k)-1}\right)^{-1}\left(\mathbf{y}-\mathbf{m}^{(k)}\right)\\\
&\boldsymbol{\Sigma}_{*}^{(k)}\overset{\Delta}{=}\mathbf{K}^{(k)}_{\mathbf{{x_{*}x_{*}}}}-\mathbf{K}^{(k)}_{\mathbf{{x_{*}x}}}\left(\mathbf{K}^{(k)}_{\mathbf{{xx}}}+\mathbf{B}^{(k)-1}\right)^{-1}\mathbf{K}^{(k)}_{\mathbf{{xx_{*}}}}+\mathbf{R}_{*}^{(k)}\\\
&\mathbf{R}_{*}^{(k)}\overset{\Delta}{=}\sigma^{2}\mathbf{I}_{\mathbf{M}}\\\
\end{split}$ (37)

The prior mixing proportions for new observations, $\boldsymbol{\varPi}[*,k]$,
weight each component equally _a priori_ , such that the sum of the weights is
equal to 1, or $\boldsymbol{\varPi}[*,k]=1/K$ [11]. The predictive equations
for the OMGP are quite similar to those for the conventional GP in Eq. (15),
with the exception of the noise component for the training data
$\mathbf{B}^{(k)-1}$ which is scaled for each sample according to the inverse
of the mixing proportion $\boldsymbol{\hat{\varPi}}[i,k]^{-1}$ [13, 11]. As
such, the contribution of each sample is weighted with respect to its
posterior predictive component in the mixture [11].

####  2.3.5 Practical implementation of the described technology

In this work, a population form was developed for the helicopter blades using
a mixture of Gaussian processes, and new (simulated) FRF data were used to
evaluate the sensitivity of the form for novelty detection. There were four
healthy helicopter blades and therefore four classes ($K=4$). Because the new
data had an equal chance of belonging to any of the four classes, the prior
mixing proportion weighted each class equally, with
$\boldsymbol{\varPi}[*,k]=0.25$. The marginal likelihood (evidence), defined
in Eq. (LABEL:eq:PredictiveEqs2) was used as a novelty index, and provided an
averaged likelihood of new data belonging to any of the $k$ classes. For
computational ease and to avoid numerical underflow, the calculations were
performed in log space. Specifically, negative log quantities were used, such
that the visual presentation of data novelty was consistent with the
Mahalanobis squared-distance (MSD) novelty index from [3, 4, 12]. As such,
data below a certain threshold were considered normal or _inyling_ , and data
above the threshold were considered novel or _outlying_.

A suitable threshold for novelty detection was determined using a similar
method as with the MSD in [3, 4, 12]. Bootstrap-sampling was used, where 1000
samples were randomly selected from the normal-condition data used to train
the form. The negative log-marginal-likelihood was then calculated for the
samples, for a large number of trials. The critical value was the threshold
with a certain percentage of the calculated values below it. This work used a
99% confidence interval, or $2.58\sigma$.

####  2.3.6 Model evaluation

The normalised mean squared-error (NMSE) and Mahalanobis squared-distance
(MSD) were calculated to compare the different models developed in this work
(e.g., supervised versus unsupervised mixture of GPs), particularly when it
was not feasible to use the model’s objective function as a performance metric
(e.g., when comparing two models learnt using different training data). The
NMSE was computed to assess each test observation against the mean of its most
likely class component $k_{*}$, via [11],

$\text{NMSE}=\frac{100}{M\sigma^{2}_{\mathbf{y_{*}}}}\left(\boldsymbol{\mu}^{(\hat{k}_{*})}_{*}-\mathbf{y_{*}}\right)^{\textsf{T}}\left(\boldsymbol{\mu}^{(\hat{k}_{*})}_{*}-\mathbf{y_{*}}\right)$
(38)

Likewise, the (normalised) MSD was computed via,

$\text{MSD}=\frac{1}{M}\sum_{i=1}^{M}\frac{\left(\boldsymbol{\mu}^{(\hat{k}_{*})}_{*}[i]-\mathbf{y_{*}}[i]\right)^{2}}{\boldsymbol{\varSigma}^{(\hat{k}_{*})}_{*}[i,i]}$
(39)

where $M$ is the number of test data points (not included in training set). In
this work, the MSD was used only to assess the predictive variance
$\boldsymbol{\varSigma}^{(k)}_{*}$. The MSD is less useful when assessing the
fit of the OMGP, as the error is scaled by the predictive variance
$\boldsymbol{\varSigma}^{(k)}_{*}$ [11].

##  3 Experimental campaign on a population of helicopter blades

Vibration data were collected over a series of tests from a set of four
healthy, full-scale composite helicopter blades, in fixed-free and free-free
configurations. This section discusses how differences among the blades and
changes in boundary conditions present as changes in the FRFs of the blades,
particularly with respect to the locations of the peaks.

###  3.1 Experiments

Vibration data were collected at ambient temperature from four healthy, full-
scale composite helicopter blades using Siemens PLM LMS SCADAS hardware and
software at the Laboratory for Verification and Validation (LVV) in Sheffield,
UK. Testing was performed in fixed-free and free-free configurations. To
approximate a fixed-free boundary condition, the root end of each blade was
placed in a substantiated strong-wall mount. To approximate a free-free
boundary condition, the blade was suspended from a gantry via springs and
cables.

Ten uniaxial 100 mV/g accelerometers were placed along the length of the
underside of each blade. Note that the same accelerometers were used on each
blade, and care was taken to ensure that they were attached to approximately
the same locations on each blade. For the fixed-free tests, an electrodynamic
shaker with force gauge was mounted to a fixture bolted to the laboratory
floor and attached to the blade 0.575 metres from the fixed end. The shaker
was attached to the underside of the blade in the flapwise direction. A
continuous random excitation was generated in LMS (note: LMS refers to Siemens
PLM LMS SCADAS hardware and software) and applied to excite the blade up to
400 Hz. Approximately 7.4 minutes of throughput force and acceleration data
were collected for each test, with a time step of 1.25e-03 seconds and a
sample rate of 800 Hz. The data were then divided into 20 blocks, each with a
dimension of 16384. Some data at the start and end of the acquisition were
discarded as these data corresponded to powering up and down of the shaker. A
Hanning window was applied to each data block and FRFs were computed in LMS.
The FRFs were then averaged in the frequency domain. For the free-free tests,
the shaker was connected at the same location on the blade as the fixed-free
tests, but was suspended from a gantry via springs. A continuous random
excitation was generated in LMS and applied to excite the blade up to 512 Hz.
Approximately 5.8 minutes of throughput force and acceleration data were
collected for each test, with a time step of 9.77e-04 seconds and a sample
rate of 1024 Hz. As with the fixed-free tests, the data were divided into 20
blocks, each with a dimension of 16384, and some data corresponding to
powering up/powering down of the shaker were discarded. A Hanning window was
applied to each data block and FRFs were computed for each measurement and
then averaged in the frequency domain. The measured bandwidths for the fixed-
free and free-free tests were selected to maximise the number and accuracy of
the captured modes (there was a large dip in the input spectrum just after 400
Hz for the fixed-free tests and just after 512 Hz for the free-free tests).
Data were collected over a relatively-long time period for each test and a
large number of averages were used to improve the spectra at the lower modes
near the limits of the operating range of the sensors. The experimental setup
for the fixed-free and free-free tests are shown in Figures 1a and 1b,
respectively. The sensor layout is shown in Figure 2. Acquisition parameters
are listed in Table 1.

An example drive-point coherence spectrum and FRF for Blade 3 in fixed-free
and free-free conditions (all computed from data obtained at the accelerometer
located near the force-input location, 0.575 metres from the blade root) are
shown in Figures 3a and 3b, respectively. These spectra are representative of
all measurements and show that coherence values were close to 1 at the
resonances, with the exception of some difficulty with collecting the first
two modes below 10 Hz. This difficulty in acquiring the first few modes
resulted in part from the operating frequency range of the accelerometers and
the extremely high flexibility of the blades, which limited the shaker
placement.

(a)

(b)

Figure 1: Helicopter blade in (a) fixed-free and (b) free-free configurations. Figure 2: Sensor locations on helicopter blades. Table 1: Acquisition parameters for vibration testing of helicopter blades. Acquisition Settings |  Fixed-Free |  Free-Free
---|---|---
Time step |  1.25e-03 s |  9.77e-04 s
Acquisition time |  20.48 s |  16.00 s
Bandwidth |  400 Hz |  512 Hz
Lines |  8192 |  8192
Frequency step |  4.88e-02 Hz |  6.25e-02 Hz
Window |  Hanning |  Hanning

(a)

(b)

Figure 3: Representative drive-point coherence and FRF from Blade 3 in (a)
fixed-free and (b) free-free configurations.

###  3.2 Discussion of experimental results

This work is concerned with demonstrating how differences among the nominally-
identical blades (with consideration for boundary condition effects), present
as changes in the peak positions as visible in the FRFs. Data were collected
from 10 accelerometers during each test. This discussion considers the results
from the accelerometer near the blade tip, as the feature changes visible in
these spectra are representative of all test points along the blade length.
The averaged FRFs for each blade in the fixed-free condition are shown in
Figures 4a and 4b. Figure 4a shows the full measured 400 Hz bandwidth and
Figure 4b shows only modes below 80 Hz. Likewise, the averaged FRF for each
blade in free-free is shown in Figures 5a and 5b, where Figure 5a shows the
full measured 512 Hz bandwidth and Figure 5b shows only modes below 80 Hz.

(a)

(b)

Figure 4: FRF magnitudes for the helicopter blades in fixed-free (a) full
bandwidth and (b) first 80 Hz.

(a)

(b)

Figure 5: FRF magnitudes for the helicopter blades in free-free (a) full
bandwidth and (b) first 80 Hz.

Figures 4a and 4b show increasing variability with respect to frequency, which
is an expected result, given that higher-frequency modes are more sensitive to
small physical changes than lower-frequency modes. For modes less than 80 Hz,
the maximum frequency difference among the blades was approximately 2.5 Hz,
while for modes greater than 80 Hz, the maximum frequency difference was
approximately 6.3 Hz. Figures 5a and 5b show a very different distribution of
peaks compared to the fixed-free tests (as expected, given such a dramatic
change in boundary condition), but also show similar variation among the
blades. For modes less than 80 Hz, the maximum frequency difference among the
blades was approximately 1.6 Hz. For modes greater than 80 Hz, the maximum
frequency difference among the blades was approximately 6.8 Hz.

Note that both the fixed-free and free-free results show grouping at several
of the peaks, where Blades 1 and 2 appear closely aligned in frequency while
Blades 3 and 4 appear closely aligned. These results are particularly
noticeable in Figure 4b, between 40 Hz and 60 Hz, and in Figure 5b, between 60
Hz and 80 Hz. These findings are quite relevant for PBSHM, which seeks to
transfer valuable information across similar structures [6]. If these data
were used to develop a normal condition for the blades against which to
evaluate new data for novelty, this variation could present some interesting
problems (e.g., negative transfer). For example, if Blade 4 experienced some
damage that reduced its natural frequency, such that it aligned with the same
peak from Blade 1, the damage may be missed. Another potential issue would be
if a new blade were added to the population and was slightly less stiff than
Blade 2, the new blade may be flagged as damaged when compared to the normal
condition developed using the original blades and test results. The
development of a population form to represent the behaviour of the group is
discussed in the next section. Only the fixed-free results were used in
developing the form, although significant changes in boundary conditions could
be addressed using techniques such as transfer learning [6]. Both the fixed-
free and free-free tests showed considerable variability among the blades and
a similar grouping pattern.

##  4 Development of a Gaussian process form

To differentiate between normal and damaged states, an important development
in PBSHM will involve establishing a normal condition for a population of
structures. In certain cases, the behaviour of the group can be represented
using a general model, called a _form_ , against which new data can be
evaluated for novelty [3, 4, 11]. This work used Gaussian process (GP)
regression to model the form for the FRFs of the helicopter blades in fixed-
free (recall that the free-free tests showed similar variability among the
blades as the fixed-free tests).

Two examples are discussed. In the first example, the form was developed using
a supervised mixture of GPs with normal-condition experimental data from the
blades. This example was supervised because it was assumed _a priori_ which
data belonged to each helicopter blade. In the second example, the normal-
condition form was developed using an (unsupervised) OMGP, where the number of
trajectory functions were assumed known _a priori_ , but it was not assumed
which data belonged to a given helicopter blade, representing a situation with
incomplete label information. This situation often occurs in SHM, when the
labels correspond to a range of operational, environmental, or damage
conditions for each structure.

###  4.1 A form using a supervised mixture of Gaussian processes

A supervised mixture of Gaussian processes assumes that the data labels are
known _a priori_. In this case, it was assumed to be known which FRF data were
associated with each helicopter blade. This assumption allowed for each GP in
the mixture to be treated independently using the equations for a conventional
GP (Eqs. (9-17)). While in a testing environment, FRFs are often described via
magnitude and phase; however, the uncertainty distribution of the magnitude is
non-Gaussian. As such, it is not proper to fit the magnitude directly using
GPs (doing so will result in support for negative values, but
$\text{mag}[\text{H}_{ij}]\geq 0$). If the time domain data are random
realisations of a stochastic process, the fast Fourier transforms (FFT) of the
signals are random realisations of the spectra, because the FFT is a linear
operator [25]. As quotients of FFTs (rather, quotients of cross-spectra and
input auto-spectra), the real and imaginary parts of the FRF are not strictly
Gaussian; however, for this paper it is assumed that each can be approximated
as GPs. In this work, the FRF was expressed as real and imaginary parts, and
the real and imaginary parts of the FRF were fit separately via Eqs. (9-17).
Note that several conference papers related to this research [26, 27],
provided a preliminary approach where magnitude was fitted directly. Although
the resulting form was useful for damage detection, the support for negative
magnitudes does not make sense and was therefore not employed in later work.

To develop the form, a narrow frequency band was selected between 48 and 56
Hz, containing the fifth bending mode of the blades. Closer inspection of the
data revealed a likely second mode in the same frequency range; however, the
modes were so closely-spaced that they could conceivably be a slightly
distorted single mode. An SDOF assumption was imposed, and an SDOF mean
function was applied to the GPs. In other words, the analysis was treated as a
grey box, in that some physics were imposed via the mean function. Simplifying
the analysis to include a single mode in the band of interest was a reasonable
assumption, as most of the physics (and variation), were accounted for with
the SDOF mean function, and the remaining distortion was resolved via the
black box component of the GP. (Note that the decision to focus on a narrow
frequency band was made to simplify the analysis while demonstrating the
proposed technology, although the form could be developed over a larger band
and with a multiple DOF assumption).

To fit the real part, the real part of the FRFs from the fixed-free tests from
each blade was copied 20 times, and random Gaussian noise with a magnitude
equal to 5% of the absolute peak value of each FRF was added to each copy. The
data copies were then concatenated into a vector, and 300 training points were
randomly selected from the vector. The modal damping FRF model from Eq. (8)
was used to generate a mean function, and the modal damping, residue, and
natural frequency were used as the mean-function hyperparameters. Upper and
lower bounds were applied to constrain the hyperparameter optimisation, based
on prior knowledge of the blades and visual inspection of the data (e.g., the
natural frequency was bounded between 40 and 60 Hz). In addition, the sign of
the data ($+/-$) was assumed to be known _a priori_ , and can easily be
determined via data inspection. This process was repeated using the imaginary
part of the FRFs, and was repeated for each blade. The posterior predictive
mean of the GP is plotted with the mean function, variance, and training data
for the real and imaginary FRFs in Figures 6b and 6a. For comparison, a single
GP was learnt for the same data (although using a different training set
comprised of 600 data points), as shown in Figures 7b and 7a.

(a)

(b)

Figure 6: Predictions for supervised mixture of GPs using (a) real and (b)
imaginary FRF data.

(a)

(b)

Figure 7: Predictions for single GP using (a) real and (b) imaginary FRF data.

Figures 6b and 6a show that the mixture provided a reasonable fit for the FRF
data, and the majority of the training data were enclosed by the variance
bounds. Furthermore, the NMSE was calculated using Eq. (38) using an
independent test data set ($M=12480$), and showed acceptable values (2.4 for
the real fit, and 6.9 for the imaginary fit). Conversely, Figures 7b and 7a
show that the differences in the natural frequencies of the blades resulted in
significant variability in the data along the horizontal axis, for which the
single GP is not well-suited. Because of this variation along the horizontal
axis, the GP was forced to apply excessive damping when fitting the data. The
NMSE was again calculated for the single GP ($M=13080$), and was
unsurprisingly, much higher than those for the mixture, at 55.8 for the real
fit and 65.1 for the imaginary fit. These results suggest that it is more
appropriate to use a mixture of GPs in cases where the data have significant
horizontal variability, as occurred here.

###  4.2 A form using an unsupervised mixture of Gaussian processes (OMGP)

An overlapping mixture of Gaussian processes (OMGP), assumes that the number
of trajectory functions are known _a priori_ , but does not assume which data
belong to a given class. In this case, because four helicopter blades
contributed to the data set, it was assumed that four trajectories (i.e.,
$K=4$ FRFs) describe the population form. The data were then classified and
fitted by approximating the intractable posterior distribution
$p(\\{\mathbf{f}^{(k)}\\}_{k=1}^{K},\mathbf{Z}\;|\;\mathbf{x},\mathbf{y})$ via
the OMGP approach described in [13], using the variational inference and EM
scheme laid out in Eqs. (18-LABEL:eq:PredictiveEqs3).

As with the previous example, a narrow frequency band containing the fifth
bending mode of the blades, located between 48 and 56 Hz, was isolated for
analysis and treated as an SDOF. Real and imaginary parts of the FRF were fit
separately via Eqs. (18-LABEL:eq:PredictiveEqs3). The same replicated data
with added noise were used as in the previous case, and 600 training points
were randomly selected from these data. Likewise, the modal damping FRF model
from Eq. (8) was used to generate a mean function, and modal damping, residue,
and natural frequency were used as the mean-function hyperparameters. Upper
and lower bounds were applied to constrain the hyperparameter optimisation,
based on prior knowledge of the blades and visual inspection of the data, and
the FRF direction was again assumed known. Because the real part of the FRF
was more separable than the imaginary part (indeed, visual inspection of the
FRFs showed that four distinct trajectories in the training data were harder
to decipher from the imaginary part versus the real part), the OMGP performed
much better for a variety of initialisations on the real data compared to the
imaginary data. As such, the real data were fitted first, and an optimal
solution was found by randomly initialising the process (i.e., via prior
guesses for the hyperparameters), ten times while keeping the training set
fixed, and choosing the solution with the maximum lower bound and therefore
the highest marginal likelihood. The initialisation parameters were chosen
randomly, but were restricted to reasonable intervals in the same manner as
the hyperparameters were bounded during the optimisation. Once learnt, the
optimised hyperparameters for the real data were used to initialise the OMGP
for the imaginary data. The posterior predictive mean of the OMGP is plotted
with the mean function, variance, and training data for the real and imaginary
FRFs in Figures 8a and 8b.

(a)

(b)

Figure 8: OMGP predictions using (a) real and (b) imaginary FRF data.

Visual inspection of Figures 8a and 8b show that the OMGP provided comparable
results to the supervised mixture from the previous example (shown in Figures
6a and 6b). Likewise, the NMSE were again calculated via Eq. (38) using an
independent test data set ($M=13080$). The calculated NMSE (2.3 for the real
fit and 6.5 for the imaginary fit) were similar to those obtained with the
supervised mixture (2.4 for the real fit, and 6.9 for the imaginary fit). The
NMSE and MSD (computed via Eq. (39)) for the OMGP are listed in Table 2, along
with those computed for the supervised mixture and single GP in the previous
example.

Table 2: Model performance metrics for (OM)GP forms. | Single GP | S. mixture | OMGP
---|---|---|---
| real | imag | real | imag | real | imag
NMSE | 55.76 | 65.05 | 2.40 | 6.92 | 2.32 | 6.49
MSD | 1.01 | 0.71 | 0.99 | 0.98 | 0.76 | 0.57

The OMGP provides a unified form that generalises the normal condition for a
population of structures, despite slight changes in boundary conditions and
differences among the nominally-identical members. Unlike the conventional GP,
that cannot accommodate large amounts of input-dependent noise, the OMGP (and
the supervised mixture of GPs) can account for the horizontal variance that
results from normal variations. The OMGP further generalises from the
supervised mixture of GPs by allowing for an unlabelled training data set.

##  5 Novelty detection via the form

One practical implementation of the OMGP form is for damage identification,
where unseen data can be compared against the posterior predictive
distributions of the OMGP and evaluated for novelty. This section demonstrates
how novelty detection could be accomplished, using a technique similar to that
involving the MSD in [3, 4, 12], but using the marginal likelihood from Eq.
(LABEL:eq:PredictiveEqs2) as a novelty index. Because the real and imaginary
parts of the FRFs were fitted separately, the sum of the (negative log)
marginal likelihood from each fit provided the novelty index. A Monte Carlo
approach (as in [12]) was used to determine an appropriate normal-condition
threshold.

To evaluate the sensitivity of the form to changes in natural frequency, real
and imaginary SDOF FRFs were synthesised using Eq. (8) to simulate stiffness
reduction resulting from damage, which may present as a downward shift in
natural frequency. FRFs were generated for each blade using modal information
from the fixed-free tests (i.e., residues and modal damping obtained via
hyperparameter optimisation from fitting the OMGP prior mean functions –
representing the ‘true’ values for these parameters) and natural frequency
incrementally decreased by 0.5% for a maximum reduction of 3.5% (i.e., the
natural frequency in the band of interest associated with each blade was
reduced such that it was 0.5%, 1%, 1.5%, 2%, 2.5%, 3%, and 3.5% lower in
frequency than the true, undamaged value). In other words, the ‘damaged’ FRFs
were very similar to the undamaged experimental FRFs, except that they were
slightly lower in frequency. The simulated FRFs were copied 1000 times and
random Gaussian noise with a magnitude equal to 5% of the absolute peak value
of the FRF was added to each copy. The negative log-marginal-likelihood was
then computed for each copy, with the novelty index equal to the negative log
sum of the marginal likelihoods for the real and imaginary parts. Likewise,
the experimental data from the fixed-free tests for each blade were copied
1000 times and random Gaussian noise with a magnitude equal to 5% of the
absolute peak value of the FRF was added to each copy. As with the simulated
data, the copied experimental data (with added noise), were tested against the
OMGP form, to evaluate the performance of the novelty detector for data known
to be inlying.

To identify a suitable normal-condition threshold, a similar technique to that
applied in [3, 4, 12] was used, where 1000 samples were randomly selected from
the normal-condition data (i.e., experimental FRFs, replicated with added
noise) used to train the form. The negative log-marginal-likelihood was then
calculated for the samples, for a large number of trials. The critical value
was the threshold with a certain percentage of the calculated values below it.
This work used a 99% confidence interval, or $2.58\sigma$.

An approximation of the FRF form magnitude was computed for visualisation
purposes. The posterior predictive distributions of the real and imaginary
OMGPs were sampled 10,000 times, and the FRF magnitudes were computed from
these samples. The mean and uncertainty bounds of the resulting magnitudes are
shown in Figures 9a-10b, along with the magnitudes of the simulated FRFs. Note
that in Figures 9a-10b, ‘Mag. mean’ refers to the mean of the FRF magnitudes
as computed from samples from the OMGP posterior while ‘Mag. bounds’ shows the
range within the magnitudes. The negative log-marginal-likelihood computed
using the experimental data from the fixed-free tests as test data and the
various synthesised FRFs as test data for Blades 1-4 are shown in Figures
11a-11d, respectively.

(a)

(b)

Figure 9: Synthesised FRFs with incrementally decreasing natural frequency for
(a) Blade 1 and (b) Blade 2, superimposed on the magnitude of the predicted
OMGP.

(a)

(b)

Figure 10: Synthesised FRFs with incrementally decreasing natural frequency
for (a) Blade 3 and (b) Blade 4, superimposed on the magnitude of the
predicted OMGP.

(a)

(b)

(c)

(d)

Figure 11: Negative log-marginal-likelihood computed using experimental and
simulated damaged results as test data for (a) Blade 1, (b) Blade 2, (c) Blade
3, and (d) Blade 4.

Examining the negative log-marginal-likelihood from right to left to
correspond with the downward frequency shift of the FRFs, Figure 11a shows
that -log$p(\mathbf{y}_{*}|\mathbf{x}_{*},\mathcal{D})$ computed using
replicated experimental FRFs for Blade 1 approached but did not surpass the
threshold (shown in Figures 11a-11d as a solid red line), which is expected,
given that the peak belonging to Blade 1 (as shown in Figure 9a) was located
at the edge of the training data with respect to frequency. As the natural
frequency of the synthesised FRFs decreased, the novelty index became
increasingly outlying. In Figure 11b, similar behaviour for Blade 2 is
visible, although the novelty index was further below the threshold for each
case as the FRF from Blade 2 had a lower magnitude than that of Blade 1, as
shown in Figures 9a and 9b. Figure 11c shows that for Blade 3, the novelty
index gradually increased until approaching the threshold at -1.5%, which
corresponded to the space between the peak groupings. As the natural frequency
of the simulated FRFs for Blade 3 decreased further (-1.5% to -2%), the
novelty index fell, because the simulated FRFs aligned with the natural
frequencies of Blades 1 and 2. As the natural frequency of the synthesised
FRFs continued to decrease (-2.5% to -3.5%), the novelty index became
increasingly outlying. Figure 11d shows similar results for Blade 4; however,
because Blade 4 had a larger FRF magnitude than Blade 3, the novelty index was
closer to the threshold. These results suggest that for Blades 1 and to a
lesser extent Blade 3, the technique would be quite sensitive to downward
frequency shifts. The technique would be less sensitive to frequency
reductions for Blades 3 and 4, as a decrease in frequency may result in
overlap with the normal-condition peaks for Blades 1 and 2. Note that this
metric was computed in a functional sense, meaning that the test data included
the full (band-limited) FRF, rather than individual test points, and the full
covariance matrices of the OMGP were used to compute the novelty index.
Considering the metric in a functional sense means that if the test data fit a
function significantly different than the form (e.g., the test data were not
an FRF), then the data could be flagged as outlying regardless of whether the
majority of individual test points were within the variance bounds.

##  6 Concluding remarks

This paper presented data obtained from a comprehensive testing campaign on
four healthy, nominally-identical helicopter blades, and used the measured
data to develop a generic model, called a _form_ , to represent the small
‘population.’ In the first example presented, a supervised mixture of GPs was
used to develop the form, where the labels were assumed known _a priori_. The
second example used an (unsupervised) OMGP to develop the form, where the data
classes were presumed unknown and needed to be inferred via the OMGP
(variational inference/EM scheme) approach. Comparable results were obtained
from both techniques. In addition, the posterior predictive distributions from
the OMGP were applied to a novelty detection approach to evaluate new
(simulated) data, providing insight into the generalisation capabilities of
the model.

Population-based SHM seeks to transfer valuable information across similar
structures, such as normal conditions and damaged states. The results
presented herein demonstrate preliminary techniques for establishing a normal
condition for a population of structures that are nominally-identical, but
distinct, with consideration for variations such as discrepancies in material
properties or fluctuations in boundary conditions. Although the proposed model
provides an accurate representation of the population data, it does not (in
its current form) consider correlation between the related functions. Further
work will consider methods of inferring correlated task parameters for
knowledge transfer within the group.

##  7 Acknowledgements

The authors gratefully acknowledge the support of the UK Engineering and
Physical Sciences Research Council (EPSRC), via grant references EP/R003645/1
and EP/R004900/1. L.A. Bull was supported by Wave 1 of the UKRI Strategic
Priorities Fund under the EPSRC grant EP/W006022/1, particularly the
Ecosystems of Digital Twins theme within that grant and the Alan Turing
Institute. For the purpose of open access, the author(s) has/have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising.

This research made use of The Laboratory for Verification and Validation
(LVV), which was funded by the EPSRC (grant numbers EP/J013714/1 and
EP/N010884/1), the European Regional Development Fund (ERDF) and the
University of Sheffield. The authors would like to extend special thanks to
Michael Dutchman, for assisting with experimental setup.

## References

  * Dowsett et al. [2018] A. Dowsett, D. O’Boy, S. Walsh, A. Abolfathi, S. Fisher,  The prediction of measurement variability in an automotive application by the use of a coherence formulation,  Proc. Inst. Mech. Eng. D: J. Automob. Eng. 232 (2018) 1694–1700. doi:10.1177/0954407017734768.
  * Chen et al. [2006] C. Chen, D. Duhamel, C. Soize,  Probabilistic approach for model and data uncertainties and its experimental identification in structural dynamics: Case of composite sandwich panels,  J. Sound Vib. 294 (2006) 64–81. doi:10.1016/j.jsv.2005.10.013.
  * L.A. Bull et al. [2021] L.A. Bull, P.A. Gardner, J. Gosliga, T.J. Rogers, N. Dervilis, E.J. Cross, E. Papatheou, A.E. Maquire, C. Campos, K. Worden,  Foundations of population-based SHM, Part I: Homogeneous populations and forms,  Mech. Syst. Signal Process. 148 (2021) 107141. doi:10.1016/j.ymssp.2020.107141.
  * L.A. Bull et al. [2019] L.A. Bull, T.J. Rogers, N. Dervilis, E.J. Cross, K. Worden,  A Gaussian process form for population-based structural health monitoring,  in: Proceedings of the 13th International Conference on Damage Assessment of Structures (DAMAS 2019), Porto, Portugal, Springer, July 9-10, 2019.
  * Gosliga et al. [2021] J. Gosliga, P.A. Gardner, L.A. Bull, N. Dervilis, K. Worden,  Foundations of population-based SHM, Part II: Heterogeneous populations–graphs, networks, and communities,  Mech. Syst. Signal Process. 148 (2021) 107144. doi:10.1016/j.ymssp.2020.107144.
  * P.A. Gardner et al. [2021] P.A. Gardner, L.A. Bull, J. Gosliga, N. Dervilis, K. Worden,  Foundations of population-based SHM, Part III: Heterogeneous populations–mapping and transfer,  Mech. Syst. Signal Process. 149 (2021) 107142. doi:10.1016/j.ymssp.2020.107142.
  * Worden et al. [2002] K. Worden, H. Sohn, C. Farrar,  Novelty detection in a changing environment: Regression and interpolation approaches,  J. Sound Vib. 258 (2002) 741–761. doi:10.1006/jsvi.2002.5148.
  * Alampalli [1998] S. Alampalli,  Influence of in-service environment on modal parameters,  in: Proceedings of the 16th International Modal Analysis Conference (IMAC), Society for Experimental Mechanics, Santa Barbara, CA, USA, 1998.
  * Cawley [1997] P. Cawley,  Long range inspection of structures using low frequency ultrasound,  in: Proceedings of the International Conference on Damage Assessment of Structures (DAMAS 97), Sheffield, UK, Sheffield Academic Press, June 30 - July 2, 1997, p. 1–17.
  * Sohn [2007] H. Sohn,  Effects of environmental and operational variability on structural health monitoring,  Philos. Trans. A: Math. Phys. Eng. Sci. 365 (2007) 539–560. doi:10.1098/rsta.2006.1935.
  * L.A. Bull et al. [2021] L.A. Bull, P.A. Gardner, T.J. Rogers, N. Dervilis, E.J. Cross, E. Papatheou, A.E. Maquire, C. Campos, K. Worden,  Bayesian modelling of multivalued power curves from an operational wind farm,  Mech. Syst. Signal Process. (2021) 108530. doi:10.1016/j.ymssp.2021.108530.
  * Worden et al. [2000] K. Worden, G. Manson, N. Fieller,  Damage detection using outlier analysis,  J. Sound Vib. 229 (2000) 647–667. doi:10.1006/JSVI.1999.2514.
  * Lázaro-Gredilla et al. [2012] M. Lázaro-Gredilla, S. Van Vaerenbergh, N.D. Lawrence,  Overlapping mixtures of Gaussian processes for the data association problem,  Pattern Recognit. 45 (2012) 1386–1395. doi:10.1016/j.patcog.2011.10.004.
  * Hong et al. [2011] S.-K. Hong, B.I. Epureanu, M.P. Castanier, D.J. Gorsich,  Parametric reduced-order models for predicting the vibration response of complex structures with component damage and uncertainties,  J. Sound Vib. 330 (2011) 1091–1110. doi:10.1016/j.jsv.2010.09.022.
  * Yang and Kessissoglou [2013] J. Yang, N. Kessissoglou,  Modal analysis of structures with uncertainties using polynomial chaos expansion,  in: Proceedings of Acoustics, Australian Acoustical Society, Victor Harbor, Australia, Nov. 17-20, 2013.
  * E.J. Cross et al. [2022] E.J. Cross, S.J. Gibson, M.R. Jones, D.J. Pitchforth, S. Zhang, T.J. Rogers, Physics-Informed Machine Learning for Structural Health Monitoring, Springer International Publishing, 2022, pp. 347–367. doi:{10.1007/978-3-030-81716-9}.
  * Worden and G.R. Tomlinson [2001] K. Worden, G.R. Tomlinson, Nonlinearity in Structural Dynamics: Detection, Identification and Modelling, Institute of Physics Publishing, Bristol, UK, 2001.
  * K.P. Murphy [2012] K.P. Murphy, Machine Learning: A Probabilistic Perspective, The MIT Press, Cambridge, MA, USA, 2012.
  * C.E. Rasmussen and C.K.I. Williams [2006] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, Adaptive Computation and Machine Learning, The MIT Press, Cambridge, MA, USA, 2006.
  * M.C.K. Tay and Laugier [2007] M.C.K. Tay, C. Laugier,  Modelling smooth paths using Gaussian processes,  in: Proceedings of the International Conference on Field and Service Robotics (FSR07), Springer, Chamonix, France, July 9-12, 2007, pp. 381–390.
  * D.M. Blei et al. [2017] D.M. Blei, A. Kucukelbir, J.D. McAuliffe,  Variational inference: A review for statisticians,  J. Am. Stat. Assoc. 112 (2017) 859–877.
  * D.J.C. MacKay [2003] D.J.C. MacKay, Information Theory, Inference and Learning Algorithms, Cambridge University Press, Cambridge, UK, 2003.
  * N.J. King and N.D. Lawrence [2006] N.J. King, N.D. Lawrence,  Fast variational inference for Gaussian process models through KL-correction,  in: Proceedings of the 17th European Conference on Machine Learning (EMCL 06), Springer-Verlag, Berlin, Germany, Sept. 18-22, 2006, pp. 270–281.
  * Lázaro-Gredilla and Titsias [2011] M. Lázaro-Gredilla, M. K. Titsias,  Variational heteroscedastic Gaussian process regression,  in: Proceedings of the 28th International Conference on Machine Learning (ICML 11), Omnipress, Bellevue, WA, USA, June 28 - July 2, 2011.
  * J.S. Bendat and A.G. Piersol [2011] J.S. Bendat, A.G. Piersol, Random Data: Analysis and Measurement Procedures (4th ed.), John Wiley & Sons, New York, NY, 2011\.
  * T.A. Dardeno et al. [2022a] T.A. Dardeno, L.A. Bull, N. Dervilis, K. Worden,  Investigating experimental repeatability and feature consistency in vibration-based SHM,  in: Proceedings of the 13th International Workshop on Structural Health Monitoring (IWSHM), 2022a.
  * T.A. Dardeno et al. [2022b] T.A. Dardeno, M. Haywood-Alexander, R.S. Mills, L.A. Bull, N. Dervilis, K. Worden,  Investigating the effects of ambient temperature on feature consistency in vibration-based SHM,  in: Proceedings of the 13th International Workshop on Structural Health Monitoring (IWSHM), 2022b.