# Efficient Subspace-Based Operational Modal Analysis Using Video-Based Vibration Measurements

Zhilei Luo, Boualem Merainani, Michael Döhler, Vincent Baltazart, Qinghua Zhang

► To cite this version:

Zhilei Luo, Boualem Merainani, Michael Döhler, Vincent Baltazart, Qinghua Zhang. Efficient Subspace-Based Operational Modal Analysis Using Video-Based Vibration Measurements. EVACES 2023 - 10th International Conference on Experimental Vibration Analysis for Civil Engineering Structures, Aug 2023, Milan, Italy. pp.32-42, 10.1007/978-3-031-39117-0_4. hal-04249311

HAL Id: hal-04249311

https://inria.hal.science/hal-04249311v1

Submitted on 19 Oct 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Efficient Subspace-based Operational Modal
Analysis Using Video-based Vibration
Measurements

Zhilei Luo¹, Boualem Merainani², Michael Döhler¹, Vincent Baltazart², and Qinghua Zhang¹

¹ Université Gustave Eiffel, Inria, COSYS-SII, I4S Team, F-35042 Rennes, France

² Université Gustave Eiffel, COSYS-SII, F-44344 Bouguenais, France

**Abstract.** Computer vision-based vibration measurement methods are contactless and offer advantages over traditional sensor measurements like accelerometers that have to be installed on the investigated structure. In particular, measurements with a high spatial resolution are obtained at relatively low cost. When processing such measurements for modal analysis with system identification methods, the high dimensional data corresponding to thousands of traditional sensors pose a challenge regarding the computational complexity and the memory requirements of the identification algorithm. In this paper, strategies for dimension reduction in subspace-based modal analysis are implemented and evaluated with regards to the obtained modal parameter uncertainties. In particular, the high spatial resolution of the mode shapes is preserved, while computation time and memory requirements are drastically reduced. The proposed method is applied to numerical and experimental data of a beam.

**Keywords:** Vibration measurements · Operational modal analysis · Dimension reduction · Vision measurements · Estimation uncertainty

# 1 Introduction

In operational modal analysis (OMA), the monitored structures are subject to natural ambient excitation that is unmeasured. Based on the measured vibration response, a linear time-invariant (LTI) state space model is identified, from where the modal parameters (natural frequencies, damping ratios and mode shapes) are obtained. Stochastic subspace methods [8] are widely used for this task, with typical methods being the covariance-driven subspace method and the Unweighted Principal Component algorithm (UPC, "data-driven") [2, 8, 9]. The popularity of subspace methods for operational modal analysis is, on one side, due to their favorable statistical properties, like consistency even when the noise driving the system is non-stationary [1]. On the other side, they are direct methods that use straightforward numerical operations like the SVD or linear least squares to obtain the estimates, which makes them numerically efficient [2].
With recently developed measurement technologies, the number of sensors that can be used for vibration analysis are rapidly increasing, and even (near) full field measurements become possible with vision-based technologies [3]. Using image processing methods, the displacement or velocity outputs are computed from the recorded videos at possibly thousands or more positions. This leads to a great spatial resolution of the identified mode shapes. The explosion of the number of outputs is a numerical challenge for the subspace methods. One way to handle a large number of sensors, so far, is to select a subset of the outputs as projection channels [4,8]. However, the information on the other outputs may not be optimally exploited. Furthermore, only one dimension of the involved covariances is reduced.

In this paper, we propose the use of projection strategies based on the PCA (considering information at all sensors) instead of the individual selection of projection channels [6]. Furthermore, both dimensions of the involved covariances are reduced in the subspace identification algorithm, in contrast to only one dimension in previous works like [4,8]. This reduces the computational burden as well as the memory usage significantly. Despite this data size reduction, the full size mode shapes are retrieved. The novel method is compared to the individual projection channel selection, and the modal parameter uncertainties are compared in Monte Carlo simulations.

The paper is organized as follows. In Section 2, the stochastic subspace identification method is briefly explained. The data reduction schemes are then introduced, including the conventional method of selecting projection channels and PCA based data reduction. In Section 3, the reduction scheme is introduced for reducing both dimensions of the involved output covariances. In Section 4, the proposed method is applied on numerical data and compared with the traditional method of individual channel selection.

# 2 Stochastic Subspace Identification with Dimension Reduction

## 2.1 Stochastic Subspace Identification and Projection Channels

It is assumed that the vibration behavior of the mechanical system is described by the discrete-time state-space model [5]

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

where $x_k \in \mathbb{R}^n$ is the state vector, $y_k \in \mathbb{R}^r$ is the output vector (displacement, velocity, or acceleration measurements), $w_k \in \mathbb{R}^n$ and $v_k \in \mathbb{R}^r$ are assumed to be white noises, corresponding to ambient excitation and measurement noise. $A \in \mathbb{R}^{n \times n}$ is the state transition matrix and $C \in \mathbb{R}^{r \times n}$ is the output matrix.

The system matrices $A$ and $C$ can be estimated with covariance-based stochastic subspace identification [8]. When dealing with high dimensional data (e.g.
large number of image pixel data), usually $r_0$ (with $r_0 \le r$) projection channels are selected [4,8]

$$y_k^{\text{red}} = L^T y_k, \quad L \in \mathbb{R}^{r \times r_0}, \qquad (2)$$

where the columns in the binary selection matrix L indicate the projection channels. Then, instead of computing the full output covariances of the measurement data, the right side reduced output covariances

$$R_i^{rr} = \mathbf{E}[y_{k+i}(y_k^{\text{red}})^T], \quad i = 1, 2, \dots \qquad (3)$$

are computed. These covariances are assembled into the block Hankel matrix $\mathcal{H}^{rr} \in \mathbb{R}^{(p+1)r \times qr_0}$, and the factorization property $R_i^{rr} = CA^{i-1}G^{\text{red}}$, $G^{\text{red}} = \mathbf{E}[x_{k+1}(y_k^{\text{red}})^T]$ [9] is utilized, then

$$\mathcal{H}^{rr} = \begin{bmatrix} R_1^{rr} & R_2^{rr} & \cdots & R_q^{rr} \\ R_2^{rr} & R_3^{rr} & \cdots & R_{q+1}^{rr} \\ \vdots & \vdots & \ddots & \vdots \\ R_{p+1}^{rr} & R_{p+2}^{rr} & \cdots & R_{q+p}^{rr} \end{bmatrix} = \mathcal{O}\mathcal{C}^{\text{red}}, \text{ with } \mathcal{O} = \begin{bmatrix} C \\ CA \\ \vdots \\ CA^p \end{bmatrix}, \qquad (4)$$

where $\mathcal{C}^{\text{red}} = [G^{\text{red}} AG^{\text{red}} \dots A^{q-1}G^{\text{red}}]$. Matrices $\mathcal{O}$ and $\mathcal{C}$ can be retrieved using the singular value decomposition, then $\mathcal{C}$ is directly obtained from the first block row of $\mathcal{O}$, and $A$ is estimated by

$$A = \mathcal{O}^{\dagger}\bar{\mathcal{O}}, \qquad (5)$$

where $\bar{\mathcal{O}}$ and $\mathcal{O} \in \mathbb{R}^{pr \times n}$ are respectively the parts of $\mathcal{O}$ after removing the first and the last block rows, and $(\cdot)^{\dagger}$ denotes the Moore-Penrose pseudoinverse.

Jacobsen et al. [4] proposed a scheme to select the projection channels based on their contribution to the data correlation. Denote the $i$th element of output vector $y_k$ by $y_k^i$, $i = 1, 2, \dots, r$. Then the output correlations yield

$$D_{ij} = \frac{\mathbf{E}(y_k^i y_k^j)}{\sqrt{\mathbf{E}(y_k^i y_k^j)\mathbf{E}(y_k^j y_k^i)}} \approx \frac{\sum_k y_k^i y_k^j}{\sqrt{\sum_k y_k^i y_k^i \sum_k y_k^j y_k^j}}. \qquad (6)$$

The first projection channel $s_1 = \arg\max_j \sum_{i=1}^r |D_{ij}|$ is selected corresponding to the column of matrix $D = [D_{ij}]_{i,j=1,\dots,r}$ with maximal absolute correlation, implying that it carries most information correlated to other sensors. The remaining projection channels ($s_m$, $m = 2,3,\dots,r_0$) are selected corresponding to the column of $D$ having the least correlation with respect to the previously selected ones

$$s_m = \arg\min_j \sum_{i \in \{s_1, \dots, s_{m-1}\}} |D_{ij}|, \quad j \neq s_1, s_2, \dots, s_{m-1}. \qquad (7)$$

In this way, most new information can be involved. The resulting projection channels compose the selection matrix

$$L_{ij} = \begin{cases} 1, & i = s_n, j = n, n = 1, 2, \dots, r_0, \\ 0, & \text{others.} \end{cases} \qquad (8)$$
## 2.2 PCA Based Dimension Reduction

Instead of choosing individual projection channels, linear combinations of the sensors can offer a more informative choice since all the sensors can actually be used and not just a subset. Optimal linear combinations can be achieved by means of Principal Component Analysis (PCA) of the correlation matrix D, i.e., based on the SVD

$$D = [U^{\text{red}} U_0] \begin{bmatrix} \Sigma^{\text{red}} & 0 \\ 0 & \Sigma_0 \end{bmatrix} \begin{bmatrix} (U^{\text{red}})^T \\ U_0^T \end{bmatrix}. \quad (9)$$

The first $r_0$ singular vectors $U^{\text{red}} \in \mathbb{R}^{r \times r_0}$ are selected for data reduction using

$$L = U^{\text{red}} \quad (10)$$

in (2), instead of binary selection matrix L in (8).

## 3 Both Sides Reduction of Hankel Matrix

The right side reduction discussed in the previous section preserves the original A and C matrices. However, only the column size of the Hankel matrix is reduced leaving the rows unchanged. Similarly, both sides (left and right) reduction can be defined on the output covariances [6]

$$R_i^{\text{br}} = \mathbf{E}[y_{k+i}^{\text{red}} (y_k^{\text{red}})^T] = C^{\text{red}} A^{i-1} G^{\text{red}} \in \mathbb{R}^{r_0 \times r_0}, \quad (11)$$

where $C^{\text{red}} = L^T C \in \mathbb{R}^{r_0 \times n}$. The corresponding Hankel matrix $\mathcal{H}^{\text{br}} \in \mathbb{R}^{(p+1)r_0 \times q r_0}$ factorizes as

$$\mathcal{H}^{\text{br}} = \begin{bmatrix} R_1^{\text{br}} & R_2^{\text{br}} & \cdots & R_q^{\text{br}} \\ R_1^{\text{br}} & R_2^{\text{br}} & \cdots & R_{q+1}^{\text{br}} \\ \vdots & \vdots & \ddots & \vdots \\ R_{p+1}^{\text{br}} & R_{p+2}^{\text{br}} & \cdots & R_{q+p}^{\text{br}} \end{bmatrix} = \mathcal{O}^{\text{red}} \mathcal{C}^{\text{red}}, \text{ where } \mathcal{O}^{\text{red}} = \begin{bmatrix} C^{\text{red}} \\ C^{\text{red}} A \\ \vdots \\ C^{\text{red}} A^p \end{bmatrix}. \quad (12)$$

Recall (5), the matrix A can still be retrieved analogously by $A = (\mathcal{O}^{\text{red}})^{\dagger} \bar{\mathcal{O}}^{\text{red}}$. However, only $\mathcal{C}^{\text{red}}$ instead of C can be obtained here. The both sides reduced Hankel matrix (12) does not contain enough information to recover a full size C. For this, only the first block row of right side reduced Hankel matrix (4) can be taken

$$\begin{aligned} \mathcal{H}_{\text{1st row}}^{\text{rr}} &= [R_1^{\text{rr}} R_2^{\text{rr}} \cdots R_q^{\text{rr}}] \\ &= C [G^{\text{red}} A G^{\text{red}} \cdots A^{q-1} G^{\text{red}}] \\ &= CC^{\text{red}}. \end{aligned} \quad (13)$$

Since $\mathcal{C}^{\text{red}}$ is already obtained in (12) when estimating A from the both sides reduced Hankel matrix, the full size C can then be estimated by

$$C = \mathcal{H}_{\text{1st row}}^{\text{rr}} (\mathcal{C}^{\text{red}})^{\dagger}. \quad (14)$$
# 4 Application

The different data size reduction schemes discussed in the previous sections are evaluated on simulated data from a finite element model of a cantilever beam:

- **FULL**. Full size identification without dimension reduction.

- **C-rr** Traditional scheme of individual projection channel selection with right side reduction on Hankel matrix (Section 2.1).

- **C-br** Individual projection channel selection with both sides reduction on Hankel matrix (Section 2.1 and 3).

- **PCACOR-rr**. PCA-based right side reduction scheme on Hankel matrix (Section 2.2).

- **PCACOR-br**. PCA-based both sides reduction scheme on Hankel matrix (Section 2.2 and 3).

The steel cantilever beam of size 30 mm × 2000 mm × 5 mm shown in Fig. 1 (left). A 3D finite element (FE) Timoshenko beam model has been made for simulation and comparison, as illustrated in Fig. 1 (right). The Young's modulus of beam is 2.1 × 10¹¹ Pa, mass density is 7850 kg/m³, and Poisson's ratio is 0.21. It is composed of 160 beam elements. Each node has three translational and three rotational degrees of freedom (DOFs), i.e. $d_x, d_y, d_z$ and $r_{xy}, r_{yz}, r_{xz}$. The DOFs at node 1 are constrained.

In a Monte Carlo simulation, $N_{MC} = 50$ displacement vibration data sets of the beam are generated from Gaussian white noise excitation at all 160 DOFs of the beam in the (horizontal) x-direction, with time step $\tau = 1/500$ for a data length of $N = 10^5$ using model (1). The datasets are processed with the proposed methods, using parameters $p + 1 = q = 61$ for the block Hankel matrices.

To test the performance of the methods, three projection channels or principal components, respectively, are selected. The computational times are compared between both reduction techniques and the full data size identification in

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Experimental setup (left) and finite element model (right) of cantilever beam.</figcaption></figure>
Table 1. The results show that the computation times are reduced by factor 40 from the full data to right side reduction schemes (C-rr and PCACOR-rr), and by another factor of 15 when using the both sides reduction schemes (C-br and PCACOR-br), demonstrating a good computational efficiency.

To examine the accuracy and precision of modal parameter estimation with the reduction techniques, three different measurement noise cases are simulated:

- **Noise case 1.** All 160 outputs are imposed with measurement noise having each 5% of the standard deviation of respective output signal.

- **Noise case 2.** Outputs at the first and last 50 nodes are imposed with noise of 5%, and the middle 60 nodes with 100% of the respective output signal standard deviation.

- **Noise case 3.** The 160 outputs are imposed with measurement noise having the same absolute amplitude for all outputs, namely 5% of the averaged standard deviation of all outputs.

The first eight structural modes are identified in all datasets. Relative bias and standard deviation of frequencies are evaluated from the Monte Carlo simulation

$$E_i = \frac{|\bar{f}_i - f_i|}{f_i}, \quad D_i = \frac{\sigma_{f_i}}{f_i}, \quad i = 1, 2, \dots, 8. \qquad (15)$$

The evaluation results are shown in Figs. 2, 3 and 4. For all noise cases, the PCA-based size reductions showed standard deviations that are either lower or in the same order as the reductions based on individual projection channels. The standard deviations with the PCA-based techniques are up to 10 times lower than with the individual projection channels for some of the modes. For all modes, standard deviations are in the same order or lower than without any size reduction (FULL), while the reductions based on the individual projection channels sometimes led to significantly higher standard deviations than FULL. While the both sides reduction with PCACOR-br yields similar or even lower standard deviations than the right side reduction only, the contrary happens for C-br. Concerning the bias, the PCA-based techniques perform for most modes similarly as FULL and better than using individual projection channels.

The individual projection channel selection during the Monte Carlo simulation and the linear combination of the sensor outputs based on the PCA for one dataset are shown in Figs. 5 and 6. In Noise case 1 and 2 the three projection channels are chosen at three distinct locations on the structure for C-rr

<table><caption>**Table 1.** Computational time of identification algorithms for one dataset.</caption><thead><tr><th>Data size reduction schemes</th><th>Computational time (seconds)</th></tr></thead><tbody><tr><td>FULL</td><td>170.4</td></tr><tr><td>C-rr / PCACOR-rr</td><td>4.8 / 4.3</td></tr><tr><td>C-br / PCACOR-br</td><td>0.3 / 0.3</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Relative bias (left) and standard deviation (right) of identified natural frequencies using 3 projection channels in **Noise case 1**.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Relative bias (left) and standard deviation (right) of identified natural frequencies using 3 projection channels in **Noise case 2**.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Relative bias (left) and standard deviation (right) of identified natural frequencies using 3 projection channels in **Noise case 3**.</figcaption></figure>

and **C-br**. However, in **Noise case 3**, always two adjacent channels are chosen which are moreover at the bottom of the beam, exhibiting the least vibration amplitudes and lowest SNR. For the PCA-based method (**PCACOR-rr** and **PCACOR-br**), the coefficients of the linear combinations are similar in the three noise cases. It can be noted that the sensors close to the bottom of the beam are hardly considered in **Noise case 3**, having the smallest SNR.

The influence of the number of chosen channels is tested in **Noise case 3** using right side reduction schemes and shown in Figs. 7 and 8. The performance of the traditional method **C-rr** increases with additional channels until around 10 channels are chosen, where it is comparable to the full size identification and becomes stable. **PCACOR-rr** has already a much better performance before ten channels selected, and its general performance is almost comparable to **FULL**.

In the PCA, the percentage of the total variance that is explained by the chosen principal components (sum of selected singular values divided by sum of all singular values) is illustrated for one dataset for each of the different noise cases in Fig. 9. **Noise case 1** has the highest initial variance percentage (using
<figure><img src="image_6.png" /><figcaption>Fig. 5. Channel selection within traditional method (C-rr and C-br) in Noise cases 1 (left), 2 (middle) and 3 (right).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Channel selection within PCA based method (PCACOR-rr and PCACOR-br) in Noise cases 1 (left), 2 (middle) and 3 (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Bias (upper) and standard deviation (lower) of identified natural frequencies with C-rr using different number of projection channels in Noise case 3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Bias (upper) and standard deviation (lower) of identified natural frequencies by PCACOR-rr using different number of principal components in Noise case 3.</figcaption></figure>
only one principal component) up to 88%, and it quickly grows to nearly 100% with only five components selected. **Noise cases 2** and **3** have a lower initial variance percentage for the first component (< 70%), which reaches nearly 100% for 63 components. It should be noted that the performance of modal parameter estimation with three principal components is already close to the FULL case (see Figs. 2–4) in all three noise cases, whereas the respective variance percentage that is explained by these three components is quite different in Fig. 9.

Finally, the algorithms have been applied to experimental data from the beam in the laboratory, where details can be found in [7, 10]. From the recorded video flows, displacements were extracted in 160 regions of interests (ROIs) as described in [7, 10], coinciding with the nodes of the numerical model. From the obtained dataset, the modal parameters were estimated with the presented methods. In Fig. 10, the identified mode shapes are plotted, which are of good accuracy and show no degradation with the reductions. The MAC values between the mode shape estimates and the numerical FE mode shapes have been evaluated, which are all above 0.99 except for mode 7, where the MAC for the

<figure><img src="image_10.png" /><figcaption>Fig. 9. Variance percentage for different numbers of principle components.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Mode shapes of beam from experimental data.</figcaption></figure>
estimates from FULL, PCACOR-rr and PCACOR-br was at 0.964 while being
at only 0.83 and 0.88 for C-rr and C-br, respectively. For the other modes, there
was no significant difference with or without the size reductions, indicating good
accuracy of the size reduction methods also on experimental data.

5 Conclusion

In this paper, a novel data reduction scheme based on principal component anal-
ysis is evaluated for subspace-based modal analysis with high dimensional output
data. The proposed method is validated not only in Monte Carlo simulations but
also on experimental vibration data from vision-based motion estimation, and
compared to traditional individual projection channel selection. The PCA-based
approaches showed better performance than the traditional method in the tests.
Even though dealing with size reduced data, the high spatial resolution and
accurate estimates of mode shapes are preserved.

Acknowledgements

The support from the ANR France Relance program is gratefully acknowledged.

References

1. Benveniste, A., Mevel, L.: Nonstationary consistency of subspace methods. IEEE Transactions on Automatic Control **AC-52**(6), 974–984 (2007)

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

3. Feng, D., Feng, M.Q.: Computer vision for SHM of civil infrastructure: From dynamic response measurement to damage detection–a review. Engineering Structures **156**, 105–117 (2018)

4. Jacobsen, N.J., Andersen, P., Brincker, R.: Applications of frequency domain curve-fitting in the EFDD technique. In: Proc. 26th IMAC (2008)

5. Juang, J.N.: Applied system identification. Prentice Hall, Englewood Cliffs, NJ, USA (1994)

6. Luo, Z., Merainani, B., Döhler, M., Baltazart, V., Zhang, Q.: High dimensional data reduction in modal analysis with stochastic subspace identification. In: Proc. 22nd IFAC World Congress (2023)

7. Merainani, B., Xiong, B., Baltazart, V., Döhler, M., Dumoulin, J., Zhang, Q.: Subspace-based modal identification and uncertainty quantification from video image flows. Journal of Sound and Vibration, submitted

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

9. Van Overschee, P., De Moor, B.: Subspace Identification for Linear Systems: Theory, Implementation, Applications. Kluwer, Dordrecht, The Netherlands (1996)

10. Xiong, B.: Video-based vibration analysis for structural health monitoring in civil engineering. Ph.D. thesis, Université Rennes 1 (Dec 2021)