# Some Elements of Operational Modal Analysis

Authors: Rune Brincker
Journal: Shock and Vibration
 (2014)

---

## Abstract

This paper gives an overview of the main components of operational modal analysis (OMA) and can serve as a tutorial for research oriented OMA applications. The paper gives a short introduction to the modeling of random responses and to the transforms often used in OMA such as the Fourier series, the Fourier integral, the Laplace transform, and the Z-transform. Then the paper introduces the spectral density matrix of the random responses and presents the theoretical solutions for correlation function and spectral density matrix under white noise loading. Some important guidelines for testing are mentioned and the most common techniques for signal processing of the operating signals are presented. The algorithms of some of the commonly used time domain and frequency domain identification techniques are presented and finally some issues are discussed such as mode shape scaling, and mode shape expansion. The different techniques are illustrated on the difficult case of identifying the three first closely spaced modes of the Heritage Court Tower building.


---

## Body

## 1. Introduction

While in traditional experimental modal analysis (EMA) the forces exciting the test specimen are controlled and normally the testing is carried out in the laboratory, in OMA the forces are just the ones that are naturally present during the operation of the structure and the test should be carried under the actual operating “in situ” conditions. For a civil engineering structure forces that might be ambient forces like wind and waves and for a mechanical structure that might be the operating forces on an engine or a gearbox, in both cases nothing is done to control temperature and other conditions that might influence the result.In OMA all modal parameters are to be determined without knowing the excitation forces. Therefore it is normally assumed that the excitation forces are Gaussian white noise, or at least that spectral densities of these forces are all flat. It is not necessary to satisfy that assumption for the actual physical forces because the physical forces acting on the structure can be thought of as created by a linear filter loaded by white noise, Figure1. In this case we maintain the assumption of white noise system input but add the properties of a linear filter to the system that is going to be identified. The properties of the filter do not change the properties of the structural system to be identified, Asmussem et al. [1], but of course we have to deal with the challenge of separating the “modes” of the loading system from the structural modes of interest.Figure 1
The commonly accepted assumption of white noise input should be thought of as loading an imaginary linear filter that produces the unknown forces. Thus the actual physical forces do not need to be white noise or have a flat spectrum.In this paper we shall focus on the theories behind the OMA technology. The paper can serve as a short overview of the present knowledge of these theories and methods. Details, derivations, and more references and information related to the treated subjects can be found in Brincker and Ventura [2] and for a more broad description of the subject, in the literature on OMA from the IMAC and IOMAC proceedings, we refer the reader to [3, 4].


## 2. Random Modeling

Since in OMA we assume the forces to be unknown we need to treat everything from a probabilistic point of view. Any parameter that we are going to observe is considered as a stochastic variable, sayX, and is in principle only known in terms of its probability density function p
(
x
). If we know the density function we can find, for instance, the mean μ
x and variance σ
x
2 as follows:
(1)
μ
x
=
∫
-
∞
∞
x
p
(
x
)
d
x
σ
x
2
=
E
[
(
X
-
μ
x
)
2
]
=
∫
-
∞
∞
(
x
-
μ
x
)
2
p
(
x
)
d
x
.
However, since we seldom will know the density function in OMA we use time averaging. If we have observed the signal x
(
t
), then we can calculate the mean and variance for that signal using time averaging as, Newland [5], follows:
(2)
μ
x
=
1
T
∫
0
T
x
(
t
)
d
t
σ
x
2
=
E
[
(
x
(
t
)
-
μ
x
)
2
]
=
1
T
∫
0
T
(
x
(
t
)
-
μ
x
)
2
d
t
.
As a main rule in OMA, we cannot use the mean values for much in practice due to large measurement errors in the low frequency region, and therefore normally we will remove the mean from the signals and calculate correlation based on the resulting zero mean signals as follows:
(3)
cor
[
x
(
t
)
,
y
(
t
)
]
=
E
[
x
(
t
)
y
(
t
)
]
=
1
T
∫
0
T
x
(
t
)
y
(
t
)
d
t
.
Since we are considering dynamic systems where the response is a linear combinations of the responses to many independent force impulses from the past, according to the central limit theorem, a random structural response is Gaussian or nearly Gaussian distributed. Therefore, since a Gaussian distribution is totally described by its second order properties (we discard the first order properties as explained above), we only need second order properties to describe random responses and in case of two signals x
(
t
) and y
(
t
) all information is contained in the correlation functions
(4)
R
x
y
(
τ
)
=
E
[
x
(
t
)
y
(
t
+
τ
)
]
=
E
[
x
(
t
-
τ
)
y
(
t
)
]
,
R
y
x
(
τ
)
=
E
[
y
(
t
)
x
(
t
+
τ
)
]
=
E
[
y
(
t
-
τ
)
x
(
t
)
]
.
The correlation functions possess the following symmetry properties that rely on the fact that we assume stationary signals and thus the time can be shifted arbitrarily so that R
x
y
(
τ
)
=
E
[
x
(
t
′
-
τ
)
y
(
t
′
)
]
=
R
y
x
(
-
τ
). In the general case of a vector response containing the individual response channels similarly we estimate the correlation function matrix as follows:
(5)
R
(
τ
)
=
E
[
y
(
t
)
y
T
(
t
+
τ
)
]
=
1
T
∫
0
T
y
(
t
)
y
T
(
t
+
τ
)
d
t
.
Again we assume stationary conditions and we can also calculate the correlation function matrix as E
[
y
(
t
-
τ
)
y
T
(
t
)
] and we get the symmetry relation for the correlation function matrix
(6)
R
(
-
τ
)
=
R
T
(
τ
)
.


## 3. Transforms

The many transforms used in signal processing and development of methods and theories in OMA constitute a problem for an easy introduction to the field. However, all the transforms are closely related. The classical Fourier series(7)
y
(
t
)
=
a
0
+
∑
k
=
1
∞
(
a
k
cos
⁡
2
π
k
t
T
+
b
k
sin
2
π
k
t
T
)
is used to describe signals with the period T. However, it is normal to write the series in the complex form
(8)
y
(
t
)
=
∑
k
=
-
∞
∞
Y
k
e
i
Δ
ω
k
t
that allows us to express the Fourier coefficients in the simple way
(9)
Y
k
=
1
T
∫
-
T
/
2
T
/
2
y
(
t
)
e
-
i
Δ
ω
k
t
d
t
.
The Fourier coefficients Y
k are discrete functions of the frequency due to the limited period T. If we extend the period to infinity we then get the Fourier integral and its inverse
(10)
y
(
t
)
=
∫
-
∞
∞
Y
(
ω
)
e
i
ω
t
d
ω
,
Y
(
ω
)
=
1
2
π
∫
-
∞
∞
y
(
t
)
e
-
i
ω
t
d
t
.
Just like periodic time functions result in discrete frequency, periodic frequency functions result in discrete time functions, thus assuming also periodic frequency functions introduce the discrete Fourier series and discrete Fourier transform
(11)
y
n
=
∑
k
=
1
N
Y
k
e
i
2
π
(
k
-
1
)
(
n
-
1
)
/
N
;
Y
k
=
1
N
∑
n
=
1
N
y
n
e
-
i
2
π
(
k
-
1
)
(
n
-
1
)
/
N
in this formulation with time and frequency shifted to comply with Matlab preferences. For the discrete time case it is important to note that according to the Shannon sampling theorem, we do not lose any information about the signal in between the sample points.The Laplace transform can be seen as a generalization of the Fourier integral where we limit our time functions to the positive axis and multiply all time functions with an exponential term so that all time functions are damped likey
′
(
t
)
=
e
-
σ
t
y
(
t
). We can then define the complex variable s
=
σ
+
i
ω and the Laplace transform and its inverse
(12)
Y
(
s
)
=
∫
0
∞
y
(
t
)
e
-
s
t
d
t
,
y
(
t
)
=
1
2
π
i
∫
σ
-
i
∞
σ
+
i
∞
Y
(
s
)
e
s
t
d
s
.
The Z-transform can be seen as the transform corresponding to the discrete series given by (8) and (9) but where we have swapped time and frequency, and thus we have continuous periodic frequency and discrete time
(13)
Y
(
z
)
=
∑
n
=
-
∞
∞
y
n
z
-
n
.
Here we have introduced the periodicity in the frequency domain by defining the complex variable z
=
e
s
Δ
t and as before for the Laplace transform we are dealing with damped version of time functions using s
=
σ
+
i
ω.The reason for using the transforms is their attractive properties that they all have in common, like, for instance, the convolution property (that convolution in the time domain corresponds to multiplication in the transform domain)(14)
h
(
t
)
*
g
(
t
)
⟷
H
(
s
)
G
(
s
)
here expressed for the Laplace transform.


## 4. Random Vibration

Random vibration is often characterized by the power spectral density (PSD) function that for a time seriesx
(
t
) is defined as the Fourier transform of the correlation function R
x
(
τ
)
(15)
G
x
y
(
ω
)
=
1
2
π
∫
-
∞
∞
R
x
y
(
τ
)
e
-
i
ω
τ
d
τ
.
The PSD is popular mainly because modes are clearly indicated by spectral peaks and according to the Parseval theorem; the area below a PSD for any frequency band is equal to the variance of the corresponding time signal (band-pass filtered to the same frequency band), and therefore the PSD has a simple physical interpretation of energy distribution (therefore the name “power” spectral density).One of the most important equations in random vibrations is the fundamental theorem that is relating the product of the PSD matrix of the input and the FRF matrix of the system to the PSD matrix of the response(16)
G
y
(
ω
)
=
H
~
*
(
i
ω
)
G
x
(
ω
)
H
~
T
(
i
ω
)
=
H
~
(
-
i
ω
)
G
x
(
ω
)
H
~
(
i
ω
)
.
The last equation follows from the identity H
~
*
(
i
ω
)
=
H
~
(
-
i
ω
) and from the fact that the transfer function is symmetric.Other central relations in random vibrations are the modal decompositions of the correlation function and PSD function matrices. The modal decomposition of the correlation function matrix is due to James et al. [6]. Expressing a general response by its modal decomposition and assuming white noise input where the correlations functions all degenerate to the Dirac delta functions it can be shown that the correlation function matrix for negative times R
y
-
(
τ
) and for positive times R
y
+
(
τ
) is given by
(17)
R
y
-
(
τ
)
=
2
π
∑
n
=
1
N
(
b
n
γ
n
T
e
-
λ
n
τ
+
b
n
*
γ
n
H
e
-
λ
n
*
τ
)
,
R
y
+
(
τ
)
=
2
π
∑
n
=
1
N
(
γ
n
b
n
T
e
λ
n
τ
+
γ
n
*
b
n
H
e
λ
n
*
τ
)
,
where b
n is the mode shape for mode n and γ
n is a vector describing the modal participation of the considered mode. It is important to note that the negative time part of the correlation function matrix R
y
-
(
τ
) is in fact a free decay because it is written as a linear combination of modal contributions (terms proportional to the mode shape times a complex exponential), whereas the positive time part R
y
+
(
τ
) is in fact only a free decay if it is used in its transposed form so the terms γ
n
b
n
T turn into the form b
n
γ
n
T and the response becomes proportional to the mode shapes. This means that whenever correlation functions are used as free decays, using the positive part of the correlation function matrix, the matrix must be used in its transposed form. It should also be noted that if we had not followed the definition given by (5) but instead defined the correlation matrix as R
(
τ
)
=
E
[
y
(
t
+
τ
)
y
T
(
t
)
] which is common in some presentations of random vibration theory, then because of stationarity this is equal to E
[
y
(
t
)
y
T
(
t
-
τ
)
]. So this swaps the time of the solutions for the correlation function matrix in (17), and thus in this case the positive part of the correlation function matrix can indeed be used as free decays without taking the transpose.The decomposition in the frequency domain can be found by taking the Fourier transform of (17) or by assuming a white noise input and using the fundamental theorem (16), Brincker et al. [7, 8],
(18)
G
y
(
ω
)
=
∑
n
=
1
N
(
b
n
γ
n
T
-
i
ω
-
λ
n
+
b
n
*
γ
n
H
-
i
ω
-
λ
n
*
+
γ
n
b
n
T
i
ω
-
λ
n
+
γ
n
*
b
n
H
i
ω
-
λ
n
*
)
.
As it appears, in the frequency domain we do not have a modal decomposition that can be considered as a linear combination of free decays from the system, because some terms are proportional to the mode shapes and some terms are proportional to the modal participation vectors. This is due to the fact that when taking the Fourier transform of (17) the terms from the negative and the positive time axes get mixed so that the first two terms in (18) are from the negative time axis and the last two ones are from the positive time axis.Including only the positive part of the correlation function matrix in the Fourier transform defines the so-called half spectrum matrix that consists only of the two last terms of (18), and, as we can see, the half spectrum matrix is a spectral representation of free decays, but again only in its transpose form.It is important to note that by using the transposed correlation function matrix as free decays in any identification technique—or using its counterpart in the frequency domain as the corresponding half spectrum—we see that we have as many free decays as we have sensors. Therefore OMA is so to speak born as a multiple input technique. From (16) we can see that a reduced rank of the input spectral density matrix G
x
(
ω
) will reduce the rank of the output spectral density accordingly, so normally it is a common assumption in OMA that the excitation of the structure is also multiple input, that is, using many independent excitation sources.


## 5. Testing

The most important concerning the testing part is to make a clear plan for the test, to secure that all measurements are carried out well, all data have the required quality, and the testing is well documented. For each data set that is to be used for OMA it can be argued, see Brincker and Ventura [2], that the total length of the time series should not be shorter than
(19)
T
tot
>
20
2
ς
f
min
⁡
=
10
ς
f
min
⁡
and the sampling frequency should not be smaller than
(20)
f
s
>
2.4
f
max
⁡
,
where f
min
⁡ is the smallest natural frequency that we are looking for and f
max
⁡ is the largest. It is important to secure a reasonable signal-to-noise ratio. This is done by making sure that the noise floor of the sensors (and the total measurement system) is well below the expected response level. Using two sensors to measure the same signal, the PSD of the measurement noise can be estimated as, Brincker and Larsen [9], follows:
(21)
G
n
(
f
)
=
(
1
-
γ
12
(
f
)
)
G
1
(
f
)
G
2
(
f
)
,
where y
1 and y
2 are the measured signals from the two sensors, G
1 and G
2 are the corresponding auto spectral densities, and γ
12 is the ordinary coherence between the two signals.


## 6. Signal Processing

As we have concluded in Section2, we extract the information from the random signals by calculating correlation functions. However before we do that we need to go through some initial preprocessing steps as follows.(i)
Validate data quality (check for clipping, dropouts, etc.).
(ii)
Calibrate signals to refer to physical units.
(iii)
Detrending (remove mean or slowly varying trend).After this initial step the user might want to evaluate and classify the operating condition during the test (for instance one or many cars on a bridge), judge the stationarity of the signals (make a time-frequency analysis), and finally evaluate the presence of harmonics (if possible remove them). Some optional preprocessing steps often used in OMA are (i)
adjustments of the sampling frequency (upsampling and downsampling, also denoted decimation),
(ii)
filtering to reduce the frequency band (low-pass, band-pass, or high-pass filters),
(iii)
integration/differentiation of signals.The different kinds of filtering can be carried out using digital FIR and IIR filters often used in electrical engineering. However, in OMA we do not need real-time filtering because common practice is to store the raw data during testing and to perform needed filtering afterwards. Because of this, FFT filters might be considered due to small phase and amplitude errors.When the preprocessing has been performed the correlation function matrix can be estimated by direct calculation according to (6) adjusting the integration to fit the total time length of the data T
=
N
Δ
t and to take the sampling into account which lead to the following simple and unbiased estimator
(22)
R
^
(
k
)
=
1
(
N
-
k
)
Δ
t
∑
n
=
1
N
-
k
y
(
n
)
y
T
(
n
+
k
)
Δ
t
=
1
N
-
k
∑
n
=
1
N
-
k
y
(
n
)
y
T
(
n
+
k
)
.
The tradition is to calculate the spectral density first by segmenting the data and using Welsh method, Brandt [10]
(23)
G
^
(
ω
)
=
1
S
∑
s
=
1
S
y
~
s
(
ω
)
*
y
~
s
T
(
ω
)
,
where y
~
s
(
ω
) is the Fourier transformed response of segment s. The correlation function can then be found by inverse Fourier transform. In OMA one should be careful with possible bias on spectral function and correlation function estimates that eventually might result in large errors in the damping values. It is worth making some comments about bias and the Welch formula for spectral estimation given by (23) because in many software implementations the Welch technique is the basis of both spectral function and correlation function estimation.It follows directly from (5) and the convolution property of the Fourier series, corresponding to (14), that (23) is only a spectral density estimate under the assumption that each data segment is periodic. Applying this technique without any windowing on the data segments corresponds to estimating the so-called circular correlation function in the time domain, and this clearly will introduce bias on the estimates. This bias is often denoted as “wrap-around” bias due to the wrong correlation so introduced between the ends of the data segment. In the frequency domain this leads to blunting of the spectral peaks and since this can be seen as energy “leaking” from the peaks to adjacent frequencies, the phenomenon is also denoted as leakage bias. The bias can be reduced giving the end points of the data segment a smaller weight by applying a windowing function on the data segments. This will reduce the leakage error but will not completely remove it.An alternative is to increase all data segments to double size by zero padding. This corresponds to assuming that the signal is zero outside of the data segment, which is also wrong and introduces a bias, but it can be shown that this bias is well defined and can be removed in the time domain by dividing the circular correlation function estimate with a triangular window. This correlation function estimate has been known since the seventies and is often denoted as the “unbiased” FFT estimate. The properties of this estimate are very close to the properties of the direct estimate given by (22).Another unbiased alternative to the direct estimation given by (22) or to zero padding is to use the random decrement technique that allows for unbiased estimating of both the correlation function and its derivative. It might also be of interest to apply the random decrement technique in cases where the user wants to apply only one single response signal to check if modal parameters depend on the excitation level of the structure.In the later years it has become popular to skip the negative part of the correlation functions so that according to (17) when taking the Fourier transform to get the spectral density the two first terms in (18) disappear and we obtain a so-called half spectrum that is a spectral representation of time domain free decays; see the discussion about this issue in Section 4. This allows for application of curve fitters known from traditional modal analysis.


## 7. Time Domain Identification

In time domain identification (TD-ID) it is normal to use parametric models obtained by least square (LS) fitting. In practice this is done by formulating an overdetermined set of equations that is solved using the pseudo inverse of the equation matrix. We will shortly summarize the ID recipes when using some popular ID algorithms like the polyreference (PR) technique, Ibrahim time domain (ITD), the eigensystem realization algorithm (ERA), and the stochastic subspace identification (SSI) technique.In PR the free decays are established (taking columns from the transposed correlation function matrix) by the correlation functions and the free decays are then arranged in a Hankel matrix, Vold et al. [11, 12],
(24)
H
1
=
[
y
(
1
)
y
(
2
)
⋯
y
(
n
p
-
n
a
)
y
(
2
)
y
(
3
)
y
(
n
p
-
(
n
a
-
1
)
)
⋮
⋮
⋱
⋮
y
(
n
a
)
y
(
n
a
+
1
)
y
(
n
p
-
1
)
]
and a “Hankel matrix” with only a single block row
(25)
H
2
=
[
y
(
n
a
+
1
)
y
(
n
a
+
2
)
⋯
y
(
n
p
)
]
.
Here the operating responses y
(
n
) are given in terms of the discrete time t
n
=
n
Δ
t. The matrix containing the AR matrices of the free decays
(26)
A
=
[
A
n
a
,
A
n
a
-
1
,
…
,
A
1
]
is then found by the LS solution
(27)
A
^
=
H
2
H
1
+
,
where H
1
+ is the pseudo inverse of H
1. The modal parameters can then be found by forming the companion matrix and performing an eigenvalue decomposition. Thus PR is an AR model-based technique. The order n
a of the AR model determines the number of modes in the model. If the number of measurement channels is n
c, then the number of rows and column of the companion matrix is n
a
×
n
c and the number of eigenvalues is then also n
a
×
n
c corresponding to n
a
×
n
c
/
2 modes.ARMA models where the response data can be modeled directly have never become popular in OMA due to the large convergence problems when several modes and channels of data are present.In ITD (in a modern formulation) a Hankel matrix is formed with four block rows, Ibrahim [13–15],
(28)
H
=
[
y
(
1
)
y
(
2
)
⋯
y
(
n
p
-
3
)
y
(
2
)
y
(
3
)
⋯
y
(
n
p
-
2
)
y
(
3
)
y
(
4
)
⋯
y
(
n
p
-
1
)
y
(
4
)
y
(
5
)
⋯
y
(
n
p
)
]
=
[
H
1
H
2
]
that is split in the middle defining H
1 and H
2. The system matrix is then simply found by the LS solution
(29)
A
^
=
H
2
H
1
T
(
H
1
H
1
T
)
-
1
,
where H
1
T
(
H
1
H
1
T
)
-
1 can be considered as the pseudo inverse of H
1. The modal parameters are found performing the eigenvalue decomposition of the system matrix A
^. This matrix defines the model order. With the previous defined variables, we see that the number of eigenvalues is 2
×
n
c, and the model has n
c number of modes. This means that the ITD model has a fixed model order corresponding to n
a
=
2 for the AR model.In ERA two Hankel matrices are formed, Juang and Pappa [16] and Pappa et al. [17, 18],
(30)
H
(
0
)
=
[
Y
(
0
)
Y
(
1
)
⋯
Y
(
1
)
Y
(
2
)
⋯
⋮
⋮
Y
(
s
-
1
)
Y
(
s
)
⋯
]
;
H
(
1
)
=
[
Y
(
1
)
Y
(
2
)
⋯
Y
(
2
)
Y
(
3
)
⋯
⋮
⋮
Y
(
s
)
Y
(
s
+
1
)
⋯
]
and an SVD is performed on the first matrix
(31)
H
(
0
)
=
U
S
V
T
.
We can then estimate the observability and controllability matrices as follows:
(32)
Γ
^
=
U
S
Λ
^
=
S
V
T
and finally the discrete time system matrix is estimated as follows:
(33)
D
^
=
Γ
^
+
H
(
1
)
Λ
^
+
.
The modal parameters are found performing the eigenvalue decomposition of D
^, but in this case the eigenvectors must be brought back to physical coordinates by the observation matrix. It should be noted that using the number of block rows of the block Hankel matrices equal to s
=
n
a, the ERA has the same number of modes as an AR model.The above mentioned techniques are all based on using the correlation functions as free decays. In SSI we use a different approach and we use the responses to construct the block Hankel matrix with2
s block rows, Overschee and de Moor [19], Peeters [20], and Peeters and de Roeck [21]
(34)
H
=
[
y
(
1
)
y
(
2
)
⋯
y
(
n
p
-
2
s
+
1
)
y
(
2
)
y
(
3
)
⋯
y
(
n
p
-
2
s
+
2
)
⋮
⋮
⋮
y
(
2
s
)
y
(
2
s
+
1
)
y
(
n
p
)
]
=
[
H
1
H
2
]
that is split in the middle defining H
1 and H
2. A projection matrix is then formed by the LS solution
(35)
O
=
H
2
H
1
T
(
H
1
H
1
T
)
+
H
1
.
Parallel to the solution idea in the ERA we now take the SVD of the projection matrix
(36)
O
=
U
S
V
T
and we estimate the observability and Kalman state matrix
(37)
Γ
^
=
U
n
S
n
,
X
^
=
S
n
V
n
T
.
The last matrix can be thought of as containing the initial conditions of the free decays in the projection matrix. Finally the discrete time system matrix and the observation matrix can be found by solving a least squares problem and the modal parameters are found performing the eigenvalue decomposition of the system matrix. The number of eigenvalues of the SSI model is equal to s
×
n
c, and the number of modes in the model is then equal to s
×
n
c
/
2.


## 8. Frequency Domain Identification

Frequency domain (FD) methods are mainly popular due to their ability to appeal to our intuition by the nice plots where we can inspect spectral peaks and have an idea about modal participation by evaluating the height of each peak and the damping by evaluating its width. But they tend to suffer from bias problems due to leakage because even though the spectral density can be estimated so that it is asymptotically unbiased (bias is zero when information approaches infinity), in practice we need to deal with a limited amount of information, and thus the leakage bias will tend to lead to overestimation of the damping in the frequency domain. Also the missing modal (free decay) decomposition of the spectral density as given by (18) is a problem because no commonly known identification techniques are designed to deal with this form, but this problem can be solved by working with the half spectral density functions. Similar to what we did for time domain we shortly summarize the ID recipes when using some popular ID algorithms in FD like classical FD (also called basic FD), frequency domain decomposition (FDD), and frequency domain polyreference (FD-PR, also denoted Polymax).In classical FD we have the simplest possible recipe, Bendat and Piersol [22] and Felber [23]:(i)
natural frequency is estimated from the location of the peak in the considered PSD;
(ii)
damping is estimated from the width of the peak;
(iii)
mode shape is estimated from any column or row in the PSD matrix.This works well only in case of well-separated modes. Also it is a problem that the user has to deal with the large number of PSD plots. In case of closely spaced modes an alternative is the FDD recipe where the PSD matrix is decomposed using SVD, Brincker et al. [7],
(38)
G
y
(
f
)
=
U
S
U
H
=
U
[
s
n
2
]
U
H
.
In this case the application is even simpler for the user because he is only inspecting a single PSD plot, that is, a plot of the singular values taken from the diagonal matrix S. These singular values can be considered as a combination of estimates of the modal coordinate auto-PSDs and the noise in the operational data. Just like in classical FD, the natural frequency can be estimated from the location of a peak in the plot, but in the FDD the mode shape is estimated as the first singular vector (first column in U) at the same frequency line. A better frequency estimate of the natural frequency (and damping) can be found by modal filtering of the PSD matrix isolating each modal coordinate in FD, taking the modal coordinate PSD to time domain and finally finding modal parameters from the 1DOF free decay, Brincker et al. [24].In FD-PR where the idea is to take the polyreference to the frequency domain, we have the problem that we cannot consider “just free decays” like in the TD, because in the FD the whole time axis is transformed into every point in the FD. Therefore in principle we have to deal with a full ARMA model in the FD; that is, the free decay in the TD resulting in a homogeneous equation of motion becomes a nonhomogenous equation of motion in the FD, Parloo [25] and Peters et al. [26, 27]. The simplest possible case is achieved assuming that the right hand side is a constant matrix (this is of course only a reasonable assumption for a narrow band estimator), Brincker and Ventura [2], and in this case the corresponding recipe is quite simple. Taking the response to be equal to the half spectrum transposed Y
(
f
)
=
G
y
T
(
f
) we form the two Hankel matrices
(39)
H
1
=
[
Y
n
a
(
f
(
k
1
)
)
,
Y
n
a
(
f
(
k
1
+
1
)
)
,
…
,
Y
n
a
(
f
(
k
2
)
)
]
H
2
=
[
I
,
I
,
…
,
I
]
.
The matrix containing the autoregressive matrices given by (25) is found by the LS solution
(40)
A
^
=
H
2
H
1
+
.
And finally the modal parameters are then found by forming the companion matrix based on the autoregressive coefficient matrices and performing an eigenvalue decomposition.


## 9. Example: The Heritage Court Tower Data

We will illustrate the OMA techniques on the well-known case of operational data from the Heritage Court Tower building. The case is described in detail in Dyck and Ventura, [28].The operational data was obtained using four datasets measuring only horizontal acceleration; two sensors close to the top of the building were used as references and the remaining sensors were then roved down the building. The first dataset applies 6 sensors placed close to the top of the building, the second dataset applies 8 sensors, the six roving sensors now moved downwards, the third dataset is similar but the roving sensors again moved downwards the building, and finally the fourth dataset includes 8 sensors where the six roving sensors are now close to the bottom of the building. All datasets have a total measurement time of 328 s and a sampling frequency of 40 Hz. Three sensors were used on each measured floor.The simplest way to graphically illustrate the operational data is to make an FDD plot. The results for the first dataset are shown in Figure2 showing the lowest quarter of the frequency band. The figure is also showing the results of the FDD performing a modal filtering taking mode shapes as the first singular vector of (38) in the points indicated in Figure 2. As it appears the first three modes are in the frequency band from 1.3 to 1.5 Hz.Figure 2
Results of using FDD to identify the first three modes on the first dataset of the HCT case. The singular values of the SD matrix are shown in dotted line. In the plot the frequency lines where the mode shape vectors are estimated are indicated by an asterisk and the corresponding modal decomposition is shown in solid line.The expected excitation of the building is a combination of wind, traffic, and excitation from people moving around in the building, so it seems reasonable to assume that the assumption of multiple input loading is fulfilled. The multiple input assumption is also supported by the fact that the SVD plot in Figure2 shows a relatively good modal separation (sufficient rank of the spectral density matrix).Since the lowest natural frequency is around 1.3 Hz, and if we assume the damping ratio to be around 1%, then the total length of each record according to (19) should be approximately 770 s. Thus the actual measurement time of 328 s is lower than the half of what is recommended by (19), and since the three first modes of this example are relatively closely spaced, it should be expected that we have some difficulties identifying the three first modes of this structure consistently. It, especially, should be expected that we have difficulties identifying the modes for the datasets where the roving sensors get close to the base of the building where the response is low and we have a decreasing signal-to-noise ratio in the measurements.The results of the identification of the first three modes of dataset 1 are shown in Table1 and the similar results for dataset 4 are shown in Table 2.Table 1
Modal identification on dataset 1 of the Heritage Court Tower case.
Modal quantity
Modal identification method
AR
ITD
ERA
SSI
FDD
FD-PR
Mode 1, f
n, (Hz)
1.227
1.227
1.227
1.227
1.225
1.228
Mode 1,ς
n, (%)
0.96
0.97
0.96
0.53
1.29
1.32
Mode 1,π
n, (%)
30.1
30.2
30.1
32.6
27.9
29.7
Mode 2, f
n, (Hz)
1.285
1.285
1.285
1.284
1.287
1.290
Mode 2,ς
n, (%)
1.19
1.20
1.19
0.58
1.41
1.46
Mode 2,π
n, (%)
19.2
19.3
19.3
18.5
17.4
19.0
Mode 3, f
n, (Hz)
1.452
1.451
1.452
1.450
1.450
1.454
Mode 3,ς
n, (%)
1.11
1.11
1.11
0.63
1.59
1.36
Mode 3,π
n, (%)
50.6
50.5
50.6
48.9
55.8
51.3Table 2
Modal identification on dataset 4 of the Heritage Court Tower case.
Modal quantity
Modal identification method
AR
ITD
ERA
SSI
FDD
FD-PR
Mode 1, f
n, (Hz)
1.246
1.246
1.246
1.215
1.235
1.213
Mode 1,ς
n, (%)
1.12
1.05
1.12
1.92
2.72
2.52
Mode 1,π
n, (%)
23.4
26.0
23.4
16.3
35.5
31.8
Mode 2, f
n, (Hz)
1.299
1.300
1.299
1.279
1.293
1.301
Mode 2,ς
n, (%)
1.57
1.40
1.57
7.20
2.45
2.62
Mode 2,π
n, (%)
3.4
3.6
3.4
10.0
0.7
0.2
Mode 3, f
n, (Hz)
1.441
1.442
1.441
1.443
1.450
1.446
Mode 3,ς
n, (%)
0.33
0.31
0.33
0.89
1.27
1.08
Mode 3,π
n, (%)
58.6
34.1
58.6
69.0
63.8
63.6For all the time domain identifications the modal participation vectorsγ
n in (17) are used to find the relative modal participation factor π
n as described in Brincker and Ventura [2]. Similarly in the frequency domain the modal participation vectors γ
n in (18) are obtained using the half spectral density matrix including only the two last terms in (18).In the time domain the first three modes were isolated using a band-pass filter with a center frequency of 1.35 Hz, a flat characteristic in a band around the center frequency with a band width of 0.4 Hz, and a roll-off band on each side with a width of 0.4 Hz.For the AR, ITD, and ERA techniques the correlation function matrices were estimated using the direct technique according to (22); the full correlation function matrix was transposed as described in Section 4 and all columns (all free decays) in the transposed correlation function matrix were then used for the identification. For the first dataset 500 discrete time lags were used in the correlation function matrix; however, for the last dataset where the modes are somewhat more difficult to identify, 950 time lags were used. For the AR, ITD, and ERA techniques a low model (corresponding to n
a
=
2) order was used including 6 modes for dataset 1 and 8 modes for dataset 4.The PC algorithm using the stochastic algorithm 1 from Overshcee and de Moor [19] was used for the SSI identification. In SSI the need for an oversized model is larger than for the previously mentioned techniques and therefore in this example a large model with s
=
80 block rows was used for the estimation. This corresponds to a model with 240 modes for the first dataset and to a model with 320 modes for dataset 4. The SVD matrices in (37) were reduced to the first 6 singular values for the first dataset (corresponding to a reduced model with only three modes) and to 18 singular values for the last dataset (corresponding to a reduced model with nine modes). It was needed to use nine modes in the last dataset in order to obtain a model including the three modes with a reasonable participation factor.In the frequency domain the modes were identified using the FDD and the FD-PR technique. In both cases the spectral density matrix was estimated taking the discrete Fourier transform of the directly estimated correlation function matrix.As it is indicated earlier, it is good idea to start any OMA with the simple FDD analysis, just looking at the plot of the singular values; see Figure2. As it appears three singular values are peaking inside the interval 1-2 Hz, while the fourth singular value is flat. Therefore the fourth singular value defines the noise floor; the first three singular values describe the physics of the system and it can be concluded that three modes are present in the considered frequency band. As mentioned earlier, the singular vectors at the three indicated points close to the three spectral peaks were chosen as mode shape estimates, and the modal coordinates were then found solving the equation y
(
t
)
=
A
q
(
t
), where the matrix A contains the estimated mode shapes and q
(
t
) are the modal coordinates. We can find the solution to this heavily overdetermined problem as the LS solution q
^
(
t
)
=
A
+
y
(
t
). This can also be done in the frequency domain, and the modal coordinate estimates are indicated by the solid lines in Figure 2. The modal coordinates estimates in the frequency domain can then be taken back to time domain by inverse FFT and the frequency and damping can be found from the modal coordinate correlation function by simple means (in this case using simple ITD for one single channel of data).The FD-PR identification was carried out based on the half spectrum matrix estimated with 1025 frequency lines in the whole frequency band up to 20 Hz. The autoregressive matrices were found using (40) over the frequency band centered at 1.35 Hz and with a band width of 0.6 Hz. The first dataset was identified using a model order of n
a
=
2 (corresponding to 6 modes for the first dataset); however, for the last dataset in order to identify all three modes with a reasonable modal participation, a model order of n
a
=
8 was used (corresponding to 32 modes).Using that each floor in the building is moving like a ridged body in the horizontal plane the movements of all points of a given floor can be estimated using only the measured three horizontal components and the different parts of the mode shape can be merged using the common reference DOFs. The mode shapes of the three first modes obtained by the different techniques are quite similar. The mode shapes of the HCT building from the AR identification are shown in Figure3.Mode shapes of the HCT building found by AR identification and merging the mode shapes from the 4 data using the reference DOFs. The units of the coordinate system are meters.
(a)
(b)
(c)The differences in the identification can be discussed on the basis of the results in Tables1 and 2. For dataset 1 it is clear that all techniques identify nearly the same natural frequencies for all three modes. Concerning the damping we can see a somewhat higher estimation uncertainty (as expected), and we can see that frequency domain techniques have a tendency to provide higher damping values than the time domain techniques. The higher damping is most probably due to leakage errors introduced by the discrete Fourier transform of the correlation functions to frequency domain. All techniques also agree that the relative modal participation is around 30% for mode 1, it is around 19% for mode 2, and it is around 50% for mode 3. The high modal participation for all three modes secures the relatively consistent identification results.Considering the results of Table2 we clearly see that this case is somewhat more difficult. We have already mentioned that in order to identify all three modes we needed to adjust the identification; for instance, the time domain techniques needed more time lags in the correlation functions, and the SSI needed to include more singular values to have more modes in the model. Even doing these adjustments in order to improve identification accuracy, we clearly now see some deviations of the natural frequencies. For the first mode we see that SSI and FD-PR provide the value f
1
=
1.215 Hz, whereas the AR, ITD, and ERA agree on f
1
=
1.246 Hz, but we know from the results of dataset 1 that the right natural frequency is f
1
=
1.227. These deviations are quite large to what would normally be expected. Modes 2 and 3 show smaller, but similar, deviations that are also larger than normally expected. Larger deviations on the damping value are observed and, for instance, the SSI provides an unrealistic estimate of more than 7% damping. These difficulties are most probably due to a small amount of data combined with a smaller signal-to-noise ratio caused by having most of the sensors close to the ground where the response is small. Another important reason to the difficulties identifying dataset 4 is that we see from Table 2 that the modal participation of mode 2 is relatively weak.This example stresses the need for good testing practice making sure that an appropriate amount of data is taken and that a good signal-to-noise ratio is present in all datasets. The analysis of this example was carried out using the MATLAB toolbox that comes with Brincker and Ventura [2].


## 10. Mode Shape Scaling and Expansion

The most commonly used method in mode shape scaling is to make a mass and/or stiffness perturbation of the test specimen and to use the corresponding change of natural frequencies and mode shape to estimate the scaling factorα that is defined as α
=
1
/
m, where m is the modal mass. Many formulas exist, but the most general one is the following, López-Aenlle et al. [29]:
(41)
α
j
2
=
(
ω
0
j
2
-
ω
I
i
2
)
t
j
i
a
0
j
T
(
ω
I
i
2
Δ
M
+
Δ
K
)
a
I
i
,
where the natural frequency ω
0
j of mode j in its unperturbed state and the corresponding mode shape a
0
j and the natural frequency ω
I
i of mode i in its perturbed state and the corresponding mode shape a
I
i are used to find an estimate of the scaling factor so that the corresponding mass scaled mode shape α
j
=
α
j
a
j can be defined. The terms t
j
i are found from the LS solution
(42)
T
=
[
t
j
i
]
=
A
0
+
A
I
,
where A
0 contains the unperturbed mode shapes and A
I contains the perturbed mode shapes. In principle this equation is exact and the only approximation is due to the estimation of the projection terms in (42). The first formulation of this kind of equation is due to Bernal [30]. However, because of the small changes of mass and stiffness that are often used due to practical reasons and because of the problems of detecting these small changes due to measurement and ID noise, the uncertainty on the scaling factor is often large and therefore it might in many cases be more accurate to expand the mode shapes and perform the scaling using the expanded mode shapes on the FE mass matrix, López-Aenlle and Brincker [31].Mode shape expansion is based on the idea of fitting a measured mode shapea with a limited number of mode shapes (subspace) from an FE model. Given a fixed subspace from the FE model with mode shapes arranged in mode shape matrix
(43)
B
=
[
B
a
B
d
]
,
where B
a contains the DOFs corresponding to the experiment (active DOF’s) and B
d contains the remaining DOFs in the FE model (deleted DOFs), we then have the classical fitting problem
(44)
a
≅
B
a
p
that can only be approximately fulfilled since we are dealing with overdetermined problem. We find the classical LS estimate for the parameter vector
(45)
p
^
=
B
a
+
a
and we have now a smoothed version of the experimental mode shape
(46)
a
^
=
B
a
p
^
.
The smoothed version can be expanded to all DOFs in the FE model just by including all DOFs in the FE mode shapes
(47)
a
^
=
B
p
^
.
The expanded version can be used for scaling as mentioned above and in damage detection and updating. If the expanded mode shapes are used for scaling together with the mass matrix of the finite element model as mentioned above, this provides a simple procedure for the scaling of OMA mode shapes. Assuming that the mode shapes from the finite element model are mass scaled, the modal mass is obtained as the inner product of the expanded experimental mode shape over the mass matrix of the finite element model
(48)
m
=
a
^
T
M
a
^
=
p
^
T
B
T
M
B
p
^
=
p
^
T
p
^
.
The result follows from (47), the well-known orthogonality principle, and the assumption of mass scaled FE modes so that the inner product B
T
M
B is equal to the identify matrix, Aenlle and Brincker, [32].Using the above mentioned expansion assuming a fixed subspace is equivalent to expansion using SEREP, O’Callahan et al. [33]. One of the problems of the expansion as outlined above is to know which modes should be included in the subspace matrix B. It is obvious that the subspace should be chosen minimal in order to obtain the best solution, (45), and thus that the optimal choice of subspace must change from mode to mode. This problem can be solved by the local correspondence (LC) principle, Brincker et al. [34], that states that any perturbed mode shape can be written as a linear combination of modes of the unperturbed system including only a few mode shapes around the corresponding unperturbed mode.


## 11. Conclusions

Some of the main elements of operational modal analysis (OMA) have been considered. It is argued that it is not necessary to assume a white noise input in order to use OMA; however, it is a central assumption in OMA that only second order information is considered (correlation and spectral density functions) and that the excitation is multiple input. The theoretical solutions for the correlation function matrix and the spectral density matrix are discussed and it is pointed out that care should be taken in order to use the second order information in a form that in fact represents free decays of the system. The identification recipes for some commonly used time domain and frequency domain techniques are presented and their ability to identify the first three closely spaced modes of the Heritage Court Tower building is illustrated. Finally the important issues of mode shape scaling and mode shape expansion are presented.



---

*Source: 325839-2014-07-06.xml*