Shock and Vibration

Volume 2014 (2014), Article ID 325839, 11 pages

http://dx.doi.org/10.1155/2014/325839

Research Article

Some Elements of Operational Modal Analysis

Rune Brincker

Department of Engineering, Aarhus University, Edison Building, Finlandsgade 22, 8200 Aarhus N, Denmark

Received 27 June 2013; Accepted 26 February 2014; Published 6 July 2014

Academic Editor: Nuno Maia

Copyright © 2014 Rune Brincker. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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, Figure 1. 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, say , and is in principle only known in terms of its probability density function . If we know the density function we can find, for instance, the mean and variance as follows:

However, since we seldom will know the density function in OMA we use time averaging. If we have observed the signal , then we can calculate the mean and variance for that signal using time averaging as, Newland [5], follows:

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:

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 and all information is contained in the correlation functions

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 . In the general case of a vector response containing the individual response channels similarly we estimate the correlation function matrix as follows:

Again we assume stationary conditions and we can also calculate the correlation function matrix as and we get the symmetry relation for the correlation function matrix

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

is used to describe signals with the period . However, it is normal to write the series in the complex form

that allows us to express the Fourier coefficients in the simple way

The Fourier coefficients are discrete functions of the frequency due to the limited period . If we extend the period to infinity we then get the Fourier integral and its inverse

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

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 like . We can then define the complex variable and the Laplace transform and its inverse

The -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

Here we have introduced the periodicity in the frequency domain by defining the complex variable and as before for the Laplace transform we are dealing with damped version of time functions using .

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)

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 series is defined as the Fourier transform of the correlation function

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

The last equation follows from the identity 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 and for positive times is given by

where is the mode shape for mode and 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 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 is in fact only a free decay if it is used in its transposed form so the terms turn into the form 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 which is common in some presentations of random vibration theory, then because of stationarity this is equal to . 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],

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 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

and the sampling frequency should not be smaller than

where is the smallest natural frequency that we are looking for and 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:

where and are the measured signals from the two sensors, and are the corresponding auto spectral densities, and is the ordinary coherence between the two signals.

6. Signal Processing

As we have concluded in Section 2, 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 and to take the sampling into account which lead to the following simple and unbiased estimator

The tradition is to calculate the spectral density first by segmenting the data and using Welsh method, Brandt [10]

where is the Fourier transformed response of segment . 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],

and a “Hankel matrix” with only a single block row

Here the operating responses are given in terms of the discrete time . The matrix containing the AR matrices of the free decays

is then found by the LS solution

where is the pseudo inverse of . 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 of the AR model determines the number of modes in the model. If the number of measurement channels is , then the number of rows and column of the companion matrix is and the number of eigenvalues is then also corresponding to 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],

that is split in the middle defining and . The system matrix is then simply found by the LS solution

where can be considered as the pseudo inverse of . The modal parameters are found performing the eigenvalue decomposition of the system matrix . This matrix defines the model order. With the previous defined variables, we see that the number of eigenvalues is , and the model has number of modes. This means that the ITD model has a fixed model order corresponding to for the AR model.

In ERA two Hankel matrices are formed, Juang and Pappa [16] and Pappa et al. [17, 18],

and an SVD is performed on the first matrix

We can then estimate the observability and controllability matrices as follows:

and finally the discrete time system matrix is estimated as follows:

The modal parameters are found performing the eigenvalue decomposition of , 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 , 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 with block rows, Overschee and de Moor [19], Peeters [20], and Peeters and de Roeck [21]

that is split in the middle defining and . A projection matrix is then formed by the LS solution

Parallel to the solution idea in the ERA we now take the SVD of the projection matrix

and we estimate the observability and Kalman state matrix

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 , and the number of modes in the model is then equal to .

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],

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 . 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 ) 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 we form the two Hankel matrices

The matrix containing the autoregressive matrices given by (25) is found by the LS solution

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 Figure 2 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 Figure 2 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 Table 1 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, , (Hz)
1.227
1.227
1.227
1.227
1.225
1.228

Mode 1, , (%)
0.96
0.97
0.96
0.53
1.29
1.32

Mode 1, , (%)
30.1
30.2
30.1
32.6
27.9
29.7

Mode 2, , (Hz)
1.285
1.285
1.285
1.284
1.287
1.290

Mode 2, , (%)
1.19
1.20
1.19
0.58
1.41
1.46

Mode 2, , (%)
19.2
19.3
19.3
18.5
17.4
19.0

Mode 3, , (Hz)
1.452
1.451
1.452
1.450
1.450
1.454

Mode 3, , (%)
1.11
1.11
1.11
0.63
1.59
1.36

Mode 3, , (%)
50.6
50.5
50.6
48.9
55.8
51.3

Table 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, , (Hz)
1.246
1.246
1.246
1.215
1.235
1.213

Mode 1, , (%)
1.12
1.05
1.12
1.92
2.72
2.52

Mode 1, , (%)
23.4
26.0
23.4
16.3
35.5
31.8

Mode 2, , (Hz)
1.299
1.300
1.299
1.279
1.293
1.301

Mode 2, , (%)
1.57
1.40
1.57
7.20
2.45
2.62

Mode 2, , (%)
3.4
3.6
3.4
10.0
0.7
0.2

Mode 3, , (Hz)
1.441
1.442
1.441
1.443
1.450
1.446

Mode 3, , (%)
0.33
0.31
0.33
0.89
1.27
1.08

Mode 3, , (%)
58.6
34.1
58.6
69.0
63.8
63.6

For all the time domain identifications the modal participation vectors in (17) are used to find the relative modal participation factor as described in Brincker and Ventura [2]. Similarly in the frequency domain the modal participation vectors 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 ) 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 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 Figure 2. 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 , where the matrix contains the estimated mode shapes and are the modal coordinates. We can find the solution to this heavily overdetermined problem as the LS solution . 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 (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 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 Figure 3.

(a)

(b)

(c)

Figure 3: 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.

The differences in the identification can be discussed on the basis of the results in Tables 1 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 Table 2 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  Hz, whereas the AR, ITD, and ERA agree on  Hz, but we know from the results of dataset 1 that the right natural frequency is . 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 , where is the modal mass. Many formulas exist, but the most general one is the following, López-Aenlle et al. [29]:

where the natural frequency of mode in its unperturbed state and the corresponding mode shape and the natural frequency of mode in its perturbed state and the corresponding mode shape are used to find an estimate of the scaling factor so that the corresponding mass scaled mode shape can be defined. The terms are found from the LS solution

where contains the unperturbed mode shapes and 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 shape 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

where contains the DOFs corresponding to the experiment (active DOF’s) and contains the remaining DOFs in the FE model (deleted DOFs), we then have the classical fitting problem

that can only be approximately fulfilled since we are dealing with overdetermined problem. We find the classical LS estimate for the parameter vector

and we have now a smoothed version of the experimental mode shape

The smoothed version can be expanded to all DOFs in the FE model just by including all DOFs in the FE mode shapes

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

The result follows from (47), the well-known orthogonality principle, and the assumption of mass scaled FE modes so that the inner product 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 . 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.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

References

J. C. Asmussen, S. R. Ibrahim, and R. Brincker, “Modal parameter identification from responses of general unknown random inputs,” in Proceedings of the 14th International Modal Analysis Conference (IMAC '96), 1996.

R. Brincker and C. E. Ventura, Introduction to Operational Modal Analysis, Wiley, 2014.

https://www.sem.org/Proceedings/ConferencePapers-Home.cfm

http://www.iomac.dk/iomac-proceedings

D. E. Newland, An Introduction to Random Vibrations, Spectral & Wavelet Analysis, Longman, 3rd edition, 1997.

G. H. James, T. G. Carne, and J. P. Lauffer, “The natural excitation technique (NExT) for modal parameter extraction from operating structures,” The International Journal of Analytical and Experimental Modal Analysis, vol. 10, no. 4, pp. 260–277, 1995.

R. Brincker, L. Zhang, and P. Andersen, “Modal identification of output-only systems using frequency domain decomposition,” Smart Materials and Structures, vol. 10, no. 3, pp. 441–445, 2000.

R. Brincker and L. Zhang, “Frequency domain decomposition revisited,” in Proceedings of the International Operational Modal Analysis Conference (IOMAC '09), pp. 615–626, 2009.

R. Brincker and J. A. Larsen, “Obtaining and estimating low noise floors in vibration sensors,” in Proceedings of the 25th Conference and Exposition on Structural Dynamics (IMAC '07), Bethel, Conn, USA, February 2007.

A. Brandt, Noise and Vibration Analysis. Signal Analysis and Experimental Procedures, Wiley, 2011.

H. Vold, J. Kundrat, G. T. Rocklin, and R. Russell, “A multi-input modal estimation algorithm for mini-computers,” SAE Paper 820194, SAE, 1982.

H. Vold and G. T. Rocklin, “The numerical implementation of a multi-input modal estimation method for mini-computers,” in Proceedings of the International Modal Analysis Conference (IMAC '82), pp. 542–548, 1982.

S. R. Ibrahim and E. C. Mikulcik, “A method for direct identification of vibration parameters from the free response,” The Shock and Vibration Bulletin, vol. 47, pp. 183–196, 1977.

S. R. Ibrahim, “Random decrement technique for modal identification of structures,” Journal of Spacecraft and Rockets, vol. 14, no. 11, pp. 696–700, 1977.

S. R. Ibrahim, “Modal confidence factor in vibration testing,” Journal of Spacecraft and Rockets, vol. 15, no. 5, pp. 313–316, 1978.

J. Juang and R. S. Pappa, “An eigen system realization algorithm for modal parameter identification and modal reduction,” Journal of Guidance, Control, and Dynamics, vol. 8, no. 5, pp. 620–627, 1985.

R. S. Pappa, K. B. Elliott, and A. Schenk, “Consistent-mode indicator for the eigensystem realization algorithm,” Journal of Guidance, Control, and Dynamics, vol. 16, no. 5, pp. 852–858, 1993.

R. S. Pappa, “Eigensystem realization algorithm, user’s guide for VAX/VMS computers,” NASA Technical Memorandum 109066, NASA, 1994.

P. Overschee and B. de Moor, Subspace Identication for Linear Systems, Theory, Implementation, Application, Kluwer Academic Publishers, 1996.

B. Peeters, System identification and damage detection in civil engineering [Ph.D. thesis], Faculteit Toegepaste Wetenschapen, Katholieke Universiteir Leuven, 2000.

B. Peeters and G. de Roeck, “Reference-based stochastic subspace identification for output-only modal analysis,” Mechanical Systems and Signal Processing, vol. 13, no. 6, pp. 855–878, 1999.

J. S. Bendat and A. G. Piersol, Engineering Applications of Correlation and Spectral Analysis, John Wiley & Sons, 2nd edition, 1993.

A. J. Felber, Development of a hybrid bridge evaluation system [Ph.D. thesis], University of British Colombia, 1993.

R. Brincker, C. E. Ventura, and P. Andersen, “Damping estimation by frequency domain decomposition,” in Proceedings of the 19th International Modal Analysis Conference (IMAC '01), pp. 698–703, Kissimmee, Fla, USA, February 2001.

E. Parloo, Application of frequency-domain system identification techniques in the field of operational modal analysis [Ph.D. thesis], Department of Mechanical Engineering, Vrije Universiteit Brussel, Brussel, Belgium, 2003.

B. Peeters, H. van der Auweraer, P. Guillaume, and J. Leuridan, “The PolyMAX frequency-domain method: a new standard for modal parameter estimation?” Shock and Vibration, vol. 11, no. 3-4, pp. 395–409, 2004.

B. Peeters, H. van der Auweraer, F. Vanhollebeke, and P. Guillaume, “Operational modal analysis for estimating the dynamic properties of a stadium structure during a football game,” Shock and Vibration, vol. 14, no. 4, pp. 283–303, 2007.

C. Dyck and C. E. Ventura, Ambient Vibration Measurements of the Heritage Court Tower, EQ LAB, University of British Columbia, Earthquake Engineering Reserach, 1998.

M. López-Aenlle, R. Brincker, F. Pelayo, and A. F. Canteli, “On exact and approximated formulations for scaling-mode shapes in operational modal analysis by mass and stiffness change,” Journal of Sound and Vibration, vol. 331, no. 3, pp. 622–637, 2012.

D. Bernal, “Modal scaling from known mass perturbations,” Journal of Engineering Mechanics, vol. 130, no. 9, pp. 1083–1088, 2004.

M. López-Aenlle and R. Brincker, “Scaling mode shapes in operational modal analysis using the finite element model mass matrix,” International Journal of Mechanical Sciences, vol. 76, pp. 86–101, 2013.

M. L. Aenlle and R. Brincker, “Modal scaling in OMA using the mass matrix of a finite element model,” in Proceedings of the 32nd International Modal Analysis Conference (IMAC '14), Orlando, Fla, USA, February 2014.

J. O'Callahan, P. Avitabile, and R. Riemer, “System equivalent reduction expansion process,” in Proceedings of the 7th International Modal Analysis Conference (IMAC '89), SEM, Las Vegas, Nev, USA, February 1989.

R. Brincker, A. Skafte, M. López-Aenlle, A. Sestieri, W. D'Ambrogio, and A. Canteli, “A local correspondence principle for mode shapes in structural dynamics,” Mechanical Systems and Signal Processing, vol. 45, no. 1–3, pp. 91–104, 2014.