Updating singular value decomposition for modal analysis of slow-varying non-stationary vibration structures

Abstract

This paper proposes a novel method for the modal analysis of slow-varying vibration structures based on vector autoregressive models. The basic idea of this method consists of using a short-time sliding window (STSW) to identify modal parameters for non-stationary vibrations. This method uses the recursive least-squares estimation for multivariable systems with the singular value decomposition (SVD) method to find the solutions within a segment of the data from each time window. Model identification is conducted by updating the SVD of the data matrix using the order and time from the previous computational window to monitor the modal parameters of a slow-varying system. Finally, this work was validated first by numerically simulating a system's gradual changes submitted to an exciting force and further by an experiment on a hydraulic turbine blade.

1. Introduction

System identification methods are generally divided into parametric and non-parametric methods (Petsounis and Fassois, 2001a and 2001b). Non-parametric methods are based on non-parameterized representations, which may be the impulse response function, the autocovariance function, and describing the signal’s power spectral density. Many studies have used the Hilbert-Huang transform (Shi et al., 2009), some focusing on the Cohen class of distribution (Lee et al., 2001; Meltzer and Ivanov, 2003; Roshan-Ghias et al., 2007), and others on the wavelet-based representation (Wang et al., 2013).

The parametric methods are the functional time-dependent autoregressive (TAR) series (Spiridonakos and Fassois, 2013), time-dependent autoregressive moving average (TARMA) series (Petsounis and Fassois, 2000), and specific functional subspaces. They have drawn much attention because of their broad application to many fields. Author Ma et al. (2018) presented the parametric output-only identification of time-varying structures using a kernel recursive extended least-squares TARMA approach. More specifically, the study used the TARMA model in kernel Hilbert space to track the time-varying dynamics. Other author Yang et al. (2015) proposed a moving kriging shape function modeling of vector TARMA models for modal identification and then validated the identification algorithm with a moving cantilever beam experiment. From another point of view, Li et al. (2019) presented a Bayesian estimation of operational modal parameters for linear time-varying mechanical systems based on the functional series vector TAR model. This built the analytical expression conjugate prior to the unknown parameters, the spanning AR coefficients, and showed the excellent performance of TAR models based on the Bayesian estimation for the time-varying vibration. Another method, employed by Spiridonakos and Fassois (2013), was to apply the stochastic functional series time-dependent autoregressive (FS-TAR) method in each state for effective fault diagnosis in inherently non-stationary structures, after which, the AR coefficients of the projection parameter vector are extracted and utilized as the characteristic quantity representing the structural state in each case. Based on parameterized representations of the time-dependent series models, parametric methods are advantageous in terms of their improved accuracy, resolution, and tracking of time-varying dynamics.

In real-life structural systems, the dynamic properties of these systems change under working conditions and are known as non-stationary systems (Chen et al., 2020). Typical structural systems include traffic-excited bridges, earthquake-excited structures, surfaces of any kind, sea vehicles, robotic devices, and rotating machinery (Au et al., 2004; Verboven et al., 2004; Vu et al., 2016). The extraction of modal parameters for non-stationary systems is more complex than for stationary systems whose dynamic properties remain constant over time.

The methods for non-stationary system identification may be classified as fast and slow non-stationary (fast and slow time-varying). In the fast non-stationary methods, the parameters are explicit functions of time. The slow non-stationary methods, on the other hand, are based on conventional stationary frequencies or time domain system identification and signal segmentation techniques. There are several difficulties in developing algorithms for identifying fast non-stationary systems. The presence of time-dependent coefficients results in more computational complexity and matrix singularity. In addition, it is very challenging to choose the functions for time-varying parameters in the methods for non-stationary vibration systems.

Therefore, the slow non-stationary methods have drawn much attention thanks to their potential application in various practical systems. Under the assumptions of short time-invariance and the theories of “time‐freezing” (Zadeh, 1950), the time-varying systems can be regarded as time-invariant systems over short periods. The vibration response signals of slow-varying systems are entirely obtained simultaneously, but they can be received by continuous sampling over time. Many researchers have developed methods to extract modal parameters for the slow-varying non-stationary vibration systems (Verboven et al., 2004; Cheng et al., 2020; Zong et al., 2020). Although many problems have been solved, the methods for slow-varying non-stationary vibration remain limited as computational complexity and matrix singularity. Hence, developing the algorithms for slow-varying non-stationary structures to overcome these disadvantages is an open problem. Furthermore, updating the computational model to track and monitor the modal parameters can be developed in conjunction with the identification method to provide a more efficient online modal analysis technique.

The time-varying autoregressive model (TVAR) is similar to the conventional autoregressive (AR) model. However, TVAR has more time-varying coefficients that could lead to certain disadvantages, such as computational complexity and matrix singularity for the identification. In addition, the time-varying system’s ambient excitation is usually difficult to measure under operating conditions. A powerful technique to mitigate these disadvantages is the locally stationary method based on the conventional stationary frequency domain or time domain system identification and signal segmentation technics. Many endeavors have been made in this direction. Author Ma and Ding (2019) assumed that the system parameters vary linearly with time in each window. A linear function describes the temporal variation of the parameter in the shifting window. Hence, the time-varying parameters are identified in the different time windows. In another identification approach, the short-time autoregressive (STAR) modeling was used for the operational modal analysis of a non-stationary mechanical system (Vu et al., 2011a). Based on the stationary state of each data segment, the modal parameter variations are monitored by autoregressive models for the emerging steel plate.

With respect to the mathematical component, the singular value decomposition method is a widely used technique to decompose a matrix into several component matrices (Gandhi and Rajgor, 2017; Stewart, 2006). It has been used in system identification (Brincker et al., 2001; Shen and Wai, 2021; Sun et al., 2021) to monitor the modal parameters for both cases: stationary and non-stationary vibrations. Other author Jiang et al. (2021) have presented a method of damage detection using singular value decomposition (SVD) for beam structures. SVD is applied to decompose the trajectory matrix of the attractor reconstructed from shape data to localize the damage to detect the defects for beam-like systems, simplifying the measurement method and reducing testing work. Another technique, found in Lobos et al. (2001), is to use singular value decomposition (SVD) to estimate the harmonics in signals in the presence of high noise. The method was developed to locate the frequencies in closely spaced sinusoidal signals. The study also presented the superiority of SVD with the standard FFT technique for signals buried in the noise. It concluded that the SVD method is especially suitable for offline analysis of recorded waveforms.

The singular value decomposition of a matrix is a valuable and important method used in the least-squares fitting of data. Many applications, for instance, signal processing, mechanical engineering, or statistics, employ SVD. In many cases, the computing procedure of SVD is repeated. This repetition could lead to high computational costs. Some authors have shared that updating the SVD is a further development that overcomes this drawback. Such an update algorithm, due to Businger, is described in Bunch and Nielsen (1978). It is reliable and efficient for a matrix with SVD and is applied when adding or deleting a row or column. In Brand (2006), author recently developed a fast algorithm to update a few dominant singular values of an augmented matrix used to perform background elimination in multiple analysis systems. A thin SVD is calculated through a matrix’s column updates and downdates.

This paper proposes a new algorithm for the online monitoring of slow-varying modal parameters in vibrating structures subjected to unknown excitation. The proposed method applies a vector autoregressive model (VAR) in a short-time sliding window (STSW) on measured signals. The model parameters are determined and updated through the order and time from the previous computational window. The recursive least-squares estimation for multivariable systems is used to find the solutions by the singular value decomposition (SVD). This work aims to avoid the computational complexity of identifying and monitoring modal parameter variations for non-stationary vibrations. In Bui et al. (2022), the Schur complement was used to update the parameters of the VAR model and monitor the varying modal parameters for a submerged plate. In terms of computational time, the Schur method is fast because it solves the standard equations of the least-squares. However, the main obstacle to vibration signal analysis is that the collected non-stationary signals are usually mixed with heavy noise caused by variable operating or environmental conditions. As a result, the rank deficient in least-squares estimation must be coped with to overcome this problem. This problem can be better handled by the generalized Schur complement (Ando, 1979). However, in this paper we are using the singular value decomposition method, which is a highly reliable, computationally stable mathematical tool that could obtain more accurate results and help to resolve these problems.

This study is structured as follows. Section 2 briefly introduces the vector autoregressive models and singular value decomposition. The updated approximation for the singular value decomposition of the matrix will be discussed in Section 3. Section 4 presents the proposed method for updating the modal parameters of the VAR model. The identification of the mechanical and operational systems will be presented in Section 5. The conclusion is summarized in the final Section.

2. Vector autoregressive models and singular value decomposition

Considering the general time-invariant recursive process for signal $\mathbf{y}\left[t\right]\in {\mathbf{R}}^{1\times n}$ , referred to as a multivariate autoregressive model at $p$ , dimension $n$ and sampling period $T_s$ , that is given by the following equation (Vu et al., 2011b):

$\mathbf{y}\left[t\right]+{\displaystyle \sum }_{i=1}^p{\mathbf{A}}_i\mathbf{y}\left[t-i\right]=\mathbf{e}\left[t\right]$

(1)

where $t$ designates the normalized discrete time, $\mathbf{e} \left[t\right]\in {\mathbf{R}}^{1\times n}$ a residual vector with zero means, and ${\mathbf{A}}_i\in {\mathbf{R}}^{n\times 1}$ the AR parameter matrix. Equation (1) is rewritten into the following linear regression form:

$\mathbf{y}\left[t\right]=\mathbf{z}\left[t\right]{\mathbf{P}}_{np\times n}+\mathbf{e}\left[t\right]$

(2)

where $\mathbf{z}\left[t\right]={\left(\mathbf{y}\left[t-1\right],\mathbf{y}\left[t-2\right],\dots ,\mathbf{y}\left[t-p\right]\right)}^T\in {\mathbf{R}}^{1\times np}$ is the corresponding regression vector, and

${\mathbf{P}}_{np\times n}=\left[-{\mathbf{A}}_1-{\mathbf{A}}_2\dots \hspace{0.17em}-{\mathbf{A}}_p\right]$ is the AR parameter matrix.

A least-squares estimation can be applied if the data are assumed to be measured in a white-noise environment. Considering N successive vectors of the output responses from $\mathbf{y} \left[t\right]$ to $\mathbf{y} \left[t+N-1\right]$ , the modal parameters matrix ${\mathbf{P}}_{np\times n}$ can be found in the least-squares method by minimizing the summed squared error between the left and right-hand sides of the equation. The objective function to be minimized may be expressed in the norm-2 vector notation form as follows (Lobos et al., 2001):

$\mathrm{E}=\frac{1}{2}{\Vert {\mathbf{K}}_{N\times np}{\mathbf{P}}_{np\times n}-\mathbf{Y}{\left[t\right]}_{N\times n}\Vert }_2^2$

(3)

where

$\mathbf{Y}{\left[t\right]}_{N\times n}=\left[\begin{array}{l}\mathbf{y}\left[t\right]\\ \mathbf{y}\left[t+1\right]\\ \dots \\ \mathbf{y}\left[t+N-1\right]\end{array}\right],{\mathbf{K}}_{N\times np}=\left[\begin{array}{l}\mathbf{z}\left[t\right]\\ \mathbf{z}\left[t+1\right]\\ \dots \\ \mathbf{z}\left[t+N-1\right]\end{array}\right]$

(4)

The singular value decomposition of the ${\mathbf{K}}_{N\times np}$ matrix is used to compute the solution of the least-squares method. There are orthogonal matrices ${\mathbf{U}}_{N\times np}$ , ${\mathbf{V}}_{np\times np}$ , and a diagonal matrix ${\mathbf{D}}_{np\times np}$ such that ${\mathbf{K}}_{N\times np}=\text{UD}{\mathbf{V}}^T$ . Here, ${\mathbf{U}}_{N\times np}$ and ${\mathbf{V}}_{np\times np}$ are the left singular vectors and the suitable right singular vectors of ${\mathbf{K}}_{N\times np}$ , respectively, and the diagonal entries $\mathbf{D}=diag\left(d_1,d_2,\dots ,d_{np}\right)$ are the singular values of ${\mathbf{K}}_{N\times np}$ . The model parameters of the AR model are estimated as follows (Lobos et al., 2001):

${\mathbf{P}}_{np\times n}=\mathbf{V}{\mathbf{D}}^{-1}{\mathbf{U}}^T\mathbf{Y}{\left[t\right]}_{N\times n}$

(5)

A state matrix of the AR model at order $p$ is constructed from the AR coefficient matrix, that is,

${\mathbf{A}}_{np\times np}=\left[\begin{array}{llllll}-{\mathbf{A}}_1& -{\mathbf{A}}_2& -{\mathbf{A}}_3& \dots & \dots & -{\mathbf{A}}_p\\ \mathbf{I}& \mathbf{0}& \mathbf{0}& \dots & \dots & \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \mathbf{0}& \dots & \dots & \mathbf{0}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \dots & \mathbf{I}& \mathbf{0}\end{array}\right]$

(6)

where $\mathbf{I}\in {\mathbf{R}}^{n\times n}$ is the identity matrix. The eigenvalue decomposition of the state matrix to determine modal parameters of a mechanic system is presented as:

${\mathbf{A}}_{np\times np}=\mathbf{L}\mathrm{\Lambda }{\mathbf{L}}^{-1}=\mathbf{L}\left|\begin{array}{llll}{g}_1& 0& 0& 0\\ 0& g_2& 0& 0\\ 0& 0& \ddots & ⋮\\ 0& 0& \dots & g_{np}\end{array}\right|{\mathbf{L}}^{-1}$

(7)

where $g_i,i=1,2,\dots ,np$ are discrete eigenvalues and $\mathbf{L}\in {\mathbf{R}}^{np\times np}$ are eigenvectors of the state matrix. Considering that, each complex eigenvalue $g_i$ of the discrete system corresponds to one natural frequency of the mechanical system, then: ${\lambda }_i=\frac{\ln \left(g_i\right)}{T_s}$ . Therefore, the natural frequencies $f_i$ and ${\xi }_i$ damping ratio are computed from complex conjugate pairs ${\lambda }_i$ as follows:

$f_i=\frac{\sqrt{{\text{Re}}^2\left({\lambda }_i\right)+{\text{Im}}^2\left({\lambda }_i\right)}}{2\pi },{\xi }_i=-\frac{\text{Re}\left({\lambda }_i\right)}{2\pi f_i}$

(8)

3. Updating the singular value decomposition of a matrix

In many least-squares and signal processing applications, one updates a matrix ${\mathbf{K}}_{N\times p}$ by appending or deleting a row or a column. After each update or downdate, the computing process of the SVD must be repeated for the resulting matrix. This section presents the updating formulation SVD for the matrix when appending and deleting a row or a column (Bunch and Nielsen, 1978).

Consider the singular value decomposition of a given matrix ${\mathbf{K}}_{N\times p}$ as:

${\mathbf{K}}_{N\times p}=\mathbf{U}\mathbf{D}{\mathbf{V}}^T$

(9)

where $\mathbf{U}\in {\mathbf{R}}^{N\times p}$ and $\mathbf{V}\in {\mathbf{R}}^{p\times p}$ are orthogonal, and $\mathbf{D}\in {\mathbf{R}}^{p\times p}$ is zero except on the main diagonal $\mathbf{D}=diag\left(d_1,d_2,\dots ,d_p\right)$ .

3.1. Updating the SVD of a matrix when appending a row

Define a new matrix ${\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}$ that is based on the given matrix ${\mathbf{K}}_{N\times p}$ when appending a row ${\mathbf{a}}^T$ as follows:

${\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}=\left(\begin{array}{l}{\mathbf{K}}_{N\times p}\\ {\mathbf{a}}^T\end{array}\right)=\stackrel{\sim }{\mathbf{U}}\stackrel{\sim }{\mathbf{D}}{\stackrel{\sim }{\mathbf{V}}}^T$

(10)

where $\stackrel{\sim }{\mathbf{U}}\in {\mathbf{R}}^{\left(N+1\right)\times p},\stackrel{\sim }{\mathbf{V}}\in {\mathbf{R}}^{p\times p},\stackrel{\sim }{\mathbf{D}}\in {\mathbf{R}}^{p\times p},\stackrel{\sim }{\mathbf{D}}=diag\left({\stackrel{\sim }{d}}_1,{\stackrel{\sim }{d}}_2,\dots ,{\stackrel{\sim }{d}}_p\right)$ .

One can compute $\stackrel{\sim }{\mathbf{U}},\stackrel{\sim }{\mathbf{D}},\stackrel{\sim }{\mathbf{V}}$ matrices of ${\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}$ by using the provided information of the matrices $\mathbf{U},\mathbf{D},\mathbf{V}$ of ${\mathbf{K}}_{N\times p}$ .

With $\mathbf{z}={\mathbf{V}}^T\mathbf{a}={\left[z_1\hspace{0.17em}{z}_2\dots \hspace{0.17em}{z}_p\right]}^T$ , a factorized representation, is defined by:

${\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}^T{\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}={\mathbf{K}}_{N\times p}^T{\mathbf{K}}_{N\times p}+\mathbf{a}{\mathbf{a}}^T=\mathbf{U}\left({\mathbf{D}}^2+\mathbf{z}{\mathbf{z}}^T\right){\mathbf{V}}^T$

(11)

From equation (11), the singular values of ${\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}$ are computed by the eigen decomposition of the matrix ${\mathbf{D}}^2+\mathbf{z}{\mathbf{z}}^T$ . The SVD of ${\mathbf{D}}^2+\mathbf{z}{\mathbf{z}}^T$ is expressed by: ${\mathbf{D}}^2+\mathbf{z}{\mathbf{z}}^T=\mathbf{Q}{\mathbf{\Omega }}^2{\mathbf{Q}}^T$ . Here, $\mathbf{Q}\in {\mathbf{R}}^{p\times p}$ is the orthogonal and $\stackrel{\sim }{\mathbf{D}}=\mathbf{\Omega }\in {\mathbf{R}}^{p\times p}$ .

The singular values of the new matrix are updated through the singular values of the matrix ${\mathbf{K}}_{N\times p}$ . The singular values of the new matrix can be updated as by the method developed in Bunch and Nielsen (1978) as follows:

${\stackrel{\sim }{d}}_i=d_i+{\mu }_i,1\le i\le p$

(12)

where ${\mu }_i$ satisfy the secular equation, more details on ${\mu }_i$ can be found in Bunch and Nielsen (1978):

$1+{\displaystyle \sum }_{i=1}^p\frac{z_i^2}{\left(d_i+d_j+\mu \right)\left(d_j-d_i-\mu \right)}=0,1\le j\le p$

(13)

Instead of computing the singular values directly of the matrix ${\stackrel{\sim }{\mathbf{K}}}_{\left(N+1\right)\times p}$ , one can update them through the secular equation (13). Once the singular values of the new matrix have been updated, the right singular vector $\stackrel{\sim }{\mathbf{V}}=\left[{\stackrel{\sim }{\mathbf{v}}}_1\hspace{0.17em}{\stackrel{\sim }{\mathbf{v}}}_2\dots \hspace{0.17em}{\stackrel{\sim }{\mathbf{v}}}_p\right]$ is obtained by:

${\stackrel{\sim }{\mathbf{v}}}_i={\frac{\mathbf{V}{\mathbf{T}}_i^{-1}{\mathbf{z}}_1}{\Vert {\mathbf{T}}_i^{-1}{\mathbf{z}}_1\Vert }}_2,{\mathbf{T}}_i={\mathbf{D}}^2-{\stackrel{\sim }{d}}_i\mathbf{I},{\mathbf{z}}_1=\frac{\mathbf{z}}{{\Vert \mathbf{a}\Vert }_2}$

(14)

The updated left singular vectors $\stackrel{\sim }{\mathbf{U}}=\left[{\stackrel{\sim }{\mathbf{u}}}_1\hspace{0.17em}{\stackrel{\sim }{\mathbf{u}}}_2\dots \hspace{0.17em}{\stackrel{\sim }{\mathbf{u}}}_p\right]$ are related to the updated right singular vectors as:

${\stackrel{\sim }{\mathbf{u}}}_i=\frac{1}{{\stackrel{\sim }{d}}_i}\left(\begin{array}{l}{\mathbf{K}}_{N\times p}\\ {\mathbf{a}}^T\end{array}\right){\stackrel{\sim }{\mathbf{v}}}_i,1\le i\le p$

(15)

Obviously, the need to repeat SVD on equation (10)’s matrix is avoided by using the updated formulas from the four equations (12)–(15), which further reduces the computational complexity of computing SVD. Thus, solving the least-squares problem by SVD is beneficial for the identification process.

3.2. Updating the SVD of a matrix when deleting a row

A new matrix based on the given matrix ${\mathbf{K}}_{N\times p}$ when deleting a row ${\mathbf{a}}^T$ is defined as equation (16):

${\mathbf{K}}_{N\times p}=\left(\begin{array}{l}{\mathbf{K}}_{\left(N-1\right)\times p}^{\%}\\ {\mathbf{a}}^T\end{array}\right)=\mathbf{U}\mathbf{D}{\mathbf{V}}^T$

(16)

where ${\stackrel{\sim }{\mathbf{K}}}_{\left(N-1\right)\times p}=\stackrel{\sim }{\mathbf{U}}\stackrel{\sim }{\mathbf{V}}\stackrel{\sim }{\mathbf{D}},\stackrel{\sim }{\mathbf{U}}\in {\mathbf{R}}^{\left(N-1\right)\times p},\stackrel{\sim }{\mathbf{V}}\in {\mathbf{R}}^{p\times p},\stackrel{\sim }{\mathbf{D}}\in {\mathbf{R}}^{p\times p},\stackrel{\sim }{\mathbf{D}}=diag\left({\stackrel{\sim }{d}}_1,{\stackrel{\sim }{d}}_2,\dots ,{\stackrel{\sim }{d}}_p\right)$

Equation (16) implies that ${\stackrel{\sim }{\mathbf{K}}}_{\left(N-1\right)\times p}^T{\stackrel{\sim }{\mathbf{K}}}_{\left(N-1\right)\times p}={\mathbf{K}}_{N\times p}^T{\mathbf{K}}_{N\times p}-\mathbf{a}{\mathbf{a}}^T=\mathbf{U}\left({\mathbf{D}}^2-\mathbf{z}{\mathbf{z}}^T\right)\mathbf{V}$ , where $={\mathbf{V}}^T\mathbf{a}={\left[z_1\hspace{0.17em}{z}_2\dots \hspace{0.17em}{z}_p\right]}^T$ .

Thus, the singular values of ${\stackrel{\sim }{\mathbf{K}}}_{\left(N-1\right)\times p}$ can be found by computing the eigen decomposition of ${\mathbf{D}}^2-\mathbf{z}{\mathbf{z}}^T=\mathbf{Q}{\mathbf{\Omega }}^2{\mathbf{Q}}^T$ , where $\mathbf{Q}\in {\mathbf{R}}^{p\times p}$ is the orthogonal matrix, and $\stackrel{\sim }{\mathbf{D}}=\mathbf{\Omega }\in {\mathbf{R}}^{p\times p}$ is a non-negative and diagonal matrix.

The singular values decomposition of the matrix ${\stackrel{\sim }{\mathbf{K}}}_{\left(N-1\right)\times p}$ is defined through the SVD of the matrix ${\mathbf{K}}_{N\times p}$ as follows (Bunch and Nielsen, 1978):

${\stackrel{\sim }{d}}_i=d_i+{\mu }_i,1\le i\le p$

(17)

where ${\mu }_i$ satisfy the secular equation:

$-1+{\displaystyle \sum }_{i=1}^p\frac{z_i^2}{\left(d_i+d_j+\mu \right)\left(d_j-d_i-\mu \right)}=0,1\le j\le p$

(18)

The updated right singular vectors $\stackrel{\sim }{\mathbf{V}}=\left[{\stackrel{\sim }{\mathbf{v}}}_1\hspace{0.17em}{\stackrel{\sim }{\mathbf{v}}}_2\dots \hspace{0.17em}{\stackrel{\sim }{\mathbf{v}}}_p\right]$ and the updated left singular vectors $\stackrel{\sim }{\mathbf{U}}=\left[{\stackrel{\sim }{\mathbf{u}}}_1\hspace{0.17em}{\stackrel{\sim }{\mathbf{u}}}_2\dots \hspace{0.17em}{\stackrel{\sim }{\mathbf{u}}}_p\right]$ are estimated as forms:

${\stackrel{\sim }{\mathbf{v}}}_i={\frac{\mathbf{V}{\mathbf{T}}_i^{-1}\mathbf{z}}{\Vert {\mathbf{T}}_i^{-1}\mathbf{z}\Vert }}_2,{\mathbf{T}}_i={\mathbf{D}}^2-{\stackrel{\sim }{d}}_i\mathbf{I},1\le i\le p$

(19)

${\stackrel{\sim }{\mathbf{u}}}_i=\frac{1}{{\stackrel{\sim }{d}}_i}{\stackrel{\sim }{\mathbf{K}}}_{\left(N-1\right)\times p}{\stackrel{\sim }{\mathbf{v}}}_i,1\le i\le p$

(20)

Thus, the set of four equations (17)–(20) is an effective procedure for updating the SVD of the matrix when deleting a row. The aim of this technique is to reduce the computational complexity when computing the SVD of a matrix.

3.3. Updating the SVD of a matrix when appending a column

In this case, a matrix based on the given matrix ${\mathbf{K}}_{N\times p}$ when adding a column is defined as:

${\stackrel{\sim }{\mathbf{K}}}_{N\times \left(p+1\right)}=\left({\mathbf{K}}_{N\times p}\mathbf{b}\right)=\stackrel{\sim }{\mathbf{U}}\stackrel{\sim }{\mathbf{D}}{\stackrel{\sim }{\mathbf{V}}}^T$

(21)

where $\stackrel{\sim }{\mathbf{U}}\in {\mathbf{R}}^{N\times \left(p+1\right)},\stackrel{\sim }{\mathbf{V}}\in {\mathbf{R}}^{\left(p+1\right)\times \left(p+1\right)},\stackrel{\sim }{\mathbf{D}}\in {\mathbf{R}}^{\left(p+1\right)\times \left(p+1\right)},\stackrel{\sim }{\mathbf{D}}=diag\left({\stackrel{\sim }{d}}_1,{\stackrel{\sim }{d}}_2,\dots ,{\stackrel{\sim }{d}}_{p+1}\right)$

Letting ${\mathbf{z}}^{\prime }={\mathbf{U}}^T\mathbf{b}={\left[z_1^{\prime }\hspace{0.17em}{z}_2^{\prime }\dots z^{\prime }\right]}^T,\alpha =\frac{1}{{\Vert \mathbf{b}\Vert }_2}$ , the factorized representation can be rewritten as:

${\stackrel{\sim }{\mathbf{K}}}_{N\times \left(p+1\right)}^T{\stackrel{\sim }{\mathbf{K}}}_{N\times \left(p+1\right)}={\mathbf{K}}_{N\times p}^T{\mathbf{K}}_{N\times p}+\mathbf{b}{\mathbf{b}}^T=\mathbf{U}\left({\mathbf{D}}^2+{\mathbf{z}}^{\prime }{{\mathbf{z}}^{\prime }}^T\right){\mathbf{V}}^T$

(22)

The eigen decomposition of ${\mathbf{D}}^2+{\mathbf{z}}^{\prime }{{\mathbf{z}}^{\prime }}^T$ can be written as $\mathbf{Q}{\mathbf{\Omega }}^2{\mathbf{Q}}^T$ , where $\mathbf{Q}\in {\mathbf{R}}^{\left(p+1\right)\times \left(p+1\right)}$ , $\stackrel{\sim }{\mathbf{D}}=\mathbf{\Omega }\in {\mathbf{R}}^{\left(p+1\right)\times \left(p+1\right)}$ .

The singular values of the matrix ${\stackrel{\sim }{\mathbf{K}}}_{N\times \left(p+1\right)}$ are updated via the eigen decomposition of ${\mathbf{D}}^2+{\mathbf{z}}^{\prime }{{\mathbf{z}}^{\prime }}^T$ as follows (Bunch and Nielsen, 1978):

${\stackrel{\sim }{d}}_i=d_i+{\mu }_i,1\le i\le p+1$

(23)

where ${\mu }_i$ satisfy the secular equation:

$1+\frac{1}{{\alpha }^2}{\displaystyle \sum }_{i=1}^p\frac{{z^{\prime }}_i^2}{\left(d_i+d_j+\mu \right)\left(d_j-d_i-\mu \right)}-\frac{{\left(1-{\Vert {\mathbf{z}}^{\prime }\Vert }_2^2\right)}^2}{{\alpha }^2\mu }=0,1\le j\le p+1$

(24)

The updated right singular vector $\stackrel{\sim }{\mathbf{V}}=\left[{\stackrel{\sim }{\mathbf{v}}}_1\hspace{0.17em}{\stackrel{\sim }{\mathbf{v}}}_2\dots \hspace{0.17em}{\stackrel{\sim }{\mathbf{v}}}_{p+1}\right]$ is determined through the left singular vector $\mathbf{V}\in {\mathbf{R}}^{p\times p}$ and the regular values $\stackrel{\sim }{\mathbf{D}}\in {\mathbf{R}}^{\left(p+1\right)\times \left(p+1\right)}$ form:

${\stackrel{\sim }{\mathbf{v}}}_i={\eta }_i\left[\begin{array}{l}\mathbf{V}{\mathbf{T}}_i^{-1}\mathbf{D}{\mathbf{z}}^{\prime }\\ -1\end{array}\right],{\mathbf{T}}_i={\mathbf{D}}^2-{\stackrel{\sim }{d}}_i\mathbf{I},{\eta }_i={\Vert \left[\begin{array}{l}\mathbf{V}{\mathbf{T}}_i^{-1}\mathbf{D}{\mathbf{z}}^{\prime }\\ -1\end{array}\right]\Vert }_2^{-1}$

(25)

The updated left singular vectors $\stackrel{\sim }{\mathbf{U}}=\left[{\stackrel{\sim }{\mathbf{u}}}_1\hspace{0.17em}{\stackrel{\sim }{\mathbf{u}}}_2\dots \hspace{0.17em}{\stackrel{\sim }{\mathbf{u}}}_{p+1}\right]$ are updated via the right singular vectors and the given matrix by:

${\stackrel{\sim }{\mathbf{u}}}_i=\frac{1}{{\stackrel{\sim }{d}}_i}\left({\mathbf{K}}_{N\times p}\hspace{0.17em}\mathbf{b}\right){\stackrel{\sim }{\mathbf{v}}}_i,1\le i\le \left(p+1\right)$

(26)

Instead of repeating the computing SVD of the matrix, the SVD updating procedure of the matrix when appending a column is presented in the set of four equations (23)–(26).

4. Modal analysis of non-stationary vibrations

As presented in the previous section, many studies have used the matrix's singular value decomposition (SVD) in system identifications. Updating the SVD of the matrices by appending or deleting a column and a row has been presented in many algorithms. This section develops a formula for updating model parameters for AR models by updating the SVD of a data matrix through the identification procedure with respect to time and model order. The proposed process for the modal analysis is also presented in this section.

4.1. Updating in model order for AR model parameters

Consider that, the data matrix ${\mathbf{K}}_{N\times np}$ and the output vector $\mathbf{Y}{\left[t\right]}_{N\times n}$ of the AR at order $p$ are formed from N successive samples by:

$\mathbf{Y}{\left[t\right]}_{N\times n}=\left[\begin{array}{l}\mathbf{y}\left[t\right]\\ \mathbf{y}\left[t+1\right]\\ \dots \\ \mathbf{y}\left[t+N-1\right]\end{array}\right],{\mathbf{K}}_{N\times np}^{\left(p\right)}=\left[\begin{array}{l}\mathbf{z}\left[t\right]\\ \mathbf{z}\left[t+1\right]\\ \dots \\ \mathbf{z}\left[t+N-1\right]\end{array}\right]$

(27)

The goal of the model identification is to determine all the model parameters. Based on equation (5), solving the least-squares problem for the model parameters of AR models by SVD is defined by:

${\mathbf{P}}_{np\times n}={\mathbf{V}}_{\left(p\right)}{\mathbf{D}}_{\left(p\right)}^{-1}{\mathbf{U}}_{\left(p\right)}^T\mathbf{Y}{\left[t\right]}_{N\times n}$

(28)

where ${\mathbf{K}}_{N\times np}^{\left(p\right)}={\mathbf{U}}_{\left(p\right)}{\mathbf{D}}_{\left(p\right)}{\mathbf{V}}_{\left(p\right)}^T$ with ${\mathbf{U}}_{\left(p\right)}\in {\mathbf{R}}^{N\times np},{\mathbf{V}}_{\left(p\right)}\in {\mathbf{R}}^{np\times np}$ are orthogonal matrices, and ${\mathbf{D}}_{\left(p\right)}\in {\mathbf{R}}^{np\times np}$ is a square diagonal matrix.

The data matrix at order $p+1$ can be built by appending a data sub-matrix into the data matrix ${\mathbf{K}}_{N\times np}^{\left(p\right)}$ as follows:

${\mathbf{K}}_{N\times n\left(p+1\right)}^{\left(p+1\right)}=\left({\mathbf{K}}_{N\times np}^{\left(p\right)}{\mathbf{K}}_{N\times n}^{\prime }\right),\text{where}{\mathbf{K}}_{N\times n}^{\prime }=\left[\begin{array}{l}\mathbf{y}\left[t-\left(p+1\right)\right]\\ \mathbf{y}\left[t+1-\left(p+1\right)\right]\\ \dots \\ \mathbf{y}\left[t+N-1-\left(p+1\right)\right]\end{array}\right]$

The model parameter estimation of AR models at order $p+1$ is obtained by updating the SVD of the data matrix. The flow diagram of updating order for model parameters is presented in Figure 1. The proposed algorithm is summarized as follows:

Figure 1.

Flow diagram of updating order for model parameters.

Input: ${\mathbf{K}}_{N\times np}^{\left(p\right)},{\mathbf{K}}_{N\times n}^{\prime },\mathbf{Y}{\left[t\right]}_{N\times n},{\mathbf{U}}_{\left(p\right)},{\mathbf{D}}_{\left(p\right)},{\mathbf{V}}_{\left(p\right)},{\mathbf{z}}^{\prime }={\mathbf{U}}^T{\mathbf{K}}_{N\times n}^{\prime }={\left[z_1^{\prime }\hspace{0.17em}{z}_2^{\prime }\dots z_{np}^{\prime }\right]}^T,\alpha ={\Vert {\mathbf{K}}_{N\times n}^{\prime }\Vert }_2^{-1}$

Step 1. Computation of the model parameters at order $p$ .

Using the SVD in equation (9) of the data matrix to estimate the model parameters at order $p$ in equation (28).

${\mathbf{P}}_{np\times n}={\mathbf{V}}_{\left(p\right)}{\mathbf{D}}_{\left(p\right)}^{-1}{\mathbf{U}}_{\left(p\right)}^T\mathbf{Y}{\left[t\right]}_{N\times n}$

(29)

Step 2. Update of the singular values of the data matrix at order $p+1$ .

Solving the secular equation in equation (24) to calculate the singular values of ${\mathbf{D}}_{\left(p+1\right)}$ via ${\mathbf{z}}^{\prime },{\mathbf{D}}_{\left(p\right)}$ at order $p$ .

Step 3. Update of the right vector ${\mathbf{V}}_{\left(p+1\right)}$ of the data matrix at order $\left(p+1\right)$ .

Based on ${\mathbf{V}}_{\left(p\right)},{\mathbf{z}}^{\prime },{\mathbf{D}}_{\left(p\right)},{\mathbf{D}}_{\left(p+1\right)}$ , the right vectors ${\mathbf{V}}_{\left(p+1\right)}$ at order $\left(p+1\right)$ are updated in equation (25).

Step 4. Update of the left vectors ${\mathbf{U}}_{\left(p+1\right)}$ of the data matrix at order $p+1$ .

Once the singular values and the right vectors at order $p+1$ are obtained, the left vectors ${\mathbf{U}}_{\left(p+1\right)}$ are updated through ${\mathbf{V}}_{\left(p\right)},{\mathbf{z}}^{\prime },{\mathbf{D}}_{\left(p\right)},{\mathbf{D}}_{\left(p+1\right)},{\mathbf{K}}_{N\times np}^{\left(p\right)},{\mathbf{K}}_{N\times n}^{\prime }$ in equation (26).

Step 5. Computation of the model parameters of AR models at order $\left(p+1\right)$ .

Because the singular values, the left and right vectors at order $p+1$ , are calculated via ${\mathbf{U}}_{\left(p\right)},{\mathbf{D}}_{\left(p\right)},{\mathbf{V}}_{\left(p\right)}^T$ at order $p$ , the solution of the least-square problem ${\mathbf{P}}_{n\left(p+1\right)\times n}={\mathbf{V}}_{\left(p+1\right)}{\mathbf{D}}_{\left(p+1\right)}^{-1}{\mathbf{U}}_{\left(p+1\right)}^T\mathbf{Y}{\left[t\right]}_{N\times n}$ can be directly updated.

The above algorithm identifies the model parameters at higher model orders through the SVD of the data matrix at the previous order. The model parameters are defined with low computational complexity and computational time effectiveness. This technique is preferable to the repetitive approach of equation (28) for each order value.

4.2. Updating in time for AR model parameters

Observation matrix $\mathbf{K}{\left[k\right]}_{N\times np}$ and output vector $\mathbf{Y}{\left[k\right]}_{N\times n}$ at time $t=k$ of the AR at order $p$ are extracted from the measured portion of the global response as follows:

$\mathbf{K}{\left[k\right]}_{N\times np}=\left[\begin{array}{l}\mathbf{z}\left[k\right]\\ \mathbf{z}\left[k+1\right]\\ \dots \\ \mathbf{z}\left[k+N-1\right]\end{array}\right],\mathbf{Y}{\left[k\right]}_{N\times n}=\left[\begin{array}{l}\mathbf{y}\left[k\right]\\ \mathbf{y}\left[k+1\right]\\ \dots \\ \mathbf{y}\left[k+N-1\right]\end{array}\right]$

(30)

Based on equation (5), solving the least-squares problem for the model parameters of AR models by SVD is defined as:

$\mathbf{P}{\left[k\right]}_{np\times n}=\mathbf{V}\left[k\right]\mathbf{D}{\left[k\right]}^{-1}\mathbf{U}{\left[k\right]}^T\mathbf{Y}{\left[k\right]}_{N\times n}$

(31)

where $\mathbf{K}{\left[k\right]}_{N\times np}=\mathbf{U}\left[k\right]\mathbf{D}\left[k\right]\mathbf{V}{\left[k\right]}^T$ . At time $t=k+m$ , the data matrix of AR models can be represented by:

$\mathbf{K}{\left[k+m\right]}_{N\times np}=\left[\begin{array}{l}\mathbf{z}\left[k+m\right]\\ \dots \\ \mathbf{z}\left[k+N+1\right]\\ \dots \\ \mathbf{z}\left[k+m+N-1\right]\end{array}\right],\mathbf{Y}{\left[k+m\right]}_{N\times n}=\left[\begin{array}{l}\mathbf{y}\left[k+m\right]\\ \mathbf{y}\left[k+m+1\right]\\ \dots \\ \mathbf{y}\left[k+m+N-1\right]\end{array}\right]$

(32)

From the data matrices of the model at time $t=k$ and $t=k+m$ , one can rewrite equations (30) and (32), respectively, as:

$\mathbf{K}{\left[k\right]}_{N\times np}=\left[\begin{array}{l}\mathbf{z}\left[k\right]\\ \mathbf{z}\left[k+1\right]\\ \dots \\ \mathbf{z}\left[k+N-1\right]\end{array}\right]=\left[\begin{array}{l}{\mathbf{R}}_{1_{m\times np}}\\ {\mathbf{R}}_{\left(N-m\right)\times np}\end{array}\right],\mathbf{K}{\left[k+m\right]}_{N\times np}=\left[\begin{array}{l}\mathbf{z}\left[k+m\right]\\ \dots \\ \mathbf{z}\left[k+N+1\right]\\ \dots \\ \mathbf{z}\left[k+m+N-1\right]\end{array}\right]=\left[\begin{array}{l}{\mathbf{R}}_{\left(N-m\right)\times np}\\ {\mathbf{R}}_{2_{m\times np}}\end{array}\right]$

(33)

where

${\mathbf{R}}_{1_{m\times np}}=\left[\begin{array}{l}\mathbf{z}\left[k\right]\\ \dots \\ \mathbf{z}\left[k+m-1\right]\end{array}\right],{\mathbf{R}}_{\left(N-m\right)\times np}=\left[\begin{array}{l}\mathbf{z}\left[k+m\right]\\ \dots \\ \mathbf{z}\left[k+N-1\right]\end{array}\right],{\mathbf{R}}_{2_{m\times np}}=\left[\begin{array}{l}\mathbf{z}\left[k+N\right]\\ \dots \\ \mathbf{z}\left[k+m+N-1\right]\end{array}\right]$

(34)

Instead of the iterative solution procedure required to calculate the model parameters at time $t=k+m$ , the model parameters of the AR model at time $t=k+m$ can be obtained by updating the SVD of the matrix. Figure 2 presents the flow diagram of updating in time for model parameters. The algorithm is implemented in five steps:

Figure 2.

Flow diagram of updating in time for model parameters.

Input: $\mathbf{K}{\left[k\right]}_{N\times np},\mathbf{Y}{\left[k\right]}_{N\times n},{{\mathbf{R}}_1}_{m\times np},{{\mathbf{R}}_2}_{m\times np},{\mathbf{R}}_{\left(N-m\right)\times np},\mathbf{Y}{\left[k+m\right]}_{N\times n},\mathbf{U}\left[k\right],\mathbf{D}\left[k\right],\mathbf{V}\left[k\right]$

Step 1Computation the model parameters at time $t=k$ .

The model parameters of the AR model at time $t=k$ are calculated using equation (31).

Step 2Update of the singular values of the data matrix at time $t=k+m$ .

The singular values of the data matrix $\mathbf{K}{\left[k+m\right]}_{N\times np}$ are identified by solving the secular equation in equations (13) and (18).

Step 3Update of the right vectors $\mathbf{V}\left[k+m\right]$ of the data matrix at time $t=k+m$ .

Once, the singular values of the data matrix $\mathbf{K}{\left[k+m\right]}_{N\times np}$ at time $t=k+m$ are obtained in Step 2. The right vector are identified by equations (19) and (14) through ${{\mathbf{R}}_1}_{m\times np},{{\mathbf{R}}_2}_{m\times np},{\mathbf{R}}_{\left(N-m\right)\times np},\mathbf{U}\left[k\right],\mathbf{D}\left[k\right],\mathbf{V}\left[k\right]$ .

Step 4Update of the left vectors $\mathbf{U}\left[k+m\right]$ of the data matrix at time $t=k+m$ .

The left vectors $\mathbf{U}\left[k+m\right]$ are computed using equations (15) and (20) through

${{\mathbf{R}}_1}_{m\times np},{{\mathbf{R}}_2}_{m\times np},{\mathbf{R}}_{\left(N-m\right)\times np},\mathbf{V}\left[k+m\right],\mathbf{D}\left[k\right],\mathbf{D}\left[k+m\right]$

Step 5Computation of the model parameters of the AR models at time $t=k+m$ .

The AR model’s model parameters are obtained by using the updated singular values, right vectors, and left vectors in steps 2, 3, and 4 as: $\mathbf{P}{\left[k+m\right]}_{np\times n}=\mathbf{V}\left[k+m\right]\mathbf{D}{\left[k+m\right]}^{-1}\mathbf{U}{\left[k+m\right]}^T\mathbf{Y}{\left[k+m\right]}_{N\times n}$ .

From this algorithm, the parameters of models $t=k+m$ are updated through the SVD of the data matrix at the time $t=k$ .

5. Results and discussion

In this section, the proposed method will be applied to a lumped-mass dynamic system and to a hydraulic turbine blade experimental setup to extract the modal parameters and monitor the systems.

5.1. Lumped-mass mechanical model

A numerical simulation was carried out to produce numerical system input-output data. A mechanical model of the simulated system used for this study is shown in Figure 3. The dynamical model of the system was derived in equation (35).

Figure 3.

Schematic model of a two-degree-of-freedom (2-DOF) time-varying mechanical system illustrating its key components and the dynamic interactions between them.

The motion equation of the system can be derived from either a Newton or Lagrange formulation, as follows:

$\left[\begin{array}{l}{m}_1\left(t\right)0\\ 0m_2\left(t\right)\end{array}\right]\left\{\begin{array}{l}{\ddot{x}}_1\left(t\right)\\ {\ddot{x}}_2\left(t\right)\end{array}\right\}+\left[\begin{array}{l}{c}_1+c_2+{\dot{m}}_1\left(t\right)\hspace{1em}-c_2\\ -c_2{c}_2+{\dot{m}}_2\left(t\right)\end{array}\right]\left\{\begin{array}{l}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\end{array}\right\}+\left[\begin{array}{l}{k}_1+k_2-k_2\\ -k_2{k}_2\end{array}\right]\left\{\begin{array}{l}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right\}=\left\{\begin{array}{l}f\left(t\right)\\ 0\end{array}\right\}$

(35)

where ${\dot{m}}_1\left(t\right)=\frac{d{m}_1\left(t\right)}{dt},{\dot{m}}_2\left(t\right)=\frac{d{m}_2\left(t\right)}{dt}$ .

Because the variation of the masses is small in this example, the derivative of the masses with respect to time is negligible in the damping matrix. Therefore, equation (35) can be rewritten as follows:

$\left[\begin{array}{l}{m}_1\left(t\right)0\\ 0m_2\left(t\right)\end{array}\right]\left\{\begin{array}{l}{\ddot{x}}_1\left(t\right)\\ {\ddot{x}}_2\left(t\right)\end{array}\right\}+\left[\begin{array}{l}{c}_1+c_2\hspace{1em}-c_2\\ -c_2\hspace{1em}{c}_2\end{array}\right]\left\{\begin{array}{l}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\end{array}\right\}+\left[\begin{array}{l}{k}_1+k_2-k_2\\ -k_2{k}_2\end{array}\right]\left\{\begin{array}{l}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right\}=\left\{\begin{array}{l}f\left(t\right)\\ 0\end{array}\right\}$

(36)

An excitation signal with a white-noise shape is applied in this sub-section with the different rate change of system’s masses. Figure 4 shows the gradual changing masses for two different rates. The numerical values of the system parameters are given as follows:

$c_1=10\left(Ns/m\right),c_2=20\left(Ns/m\right),k_1=10000\left(N/m\right),k_2=22000\left(N/m\right)$

Figure 4.

Mass changing function.

The system’s varying modal parameters are simulated in Figure 5 on a theoretical basis for the two different rates, 2.5 (%/s) and 5 (%/s). At the masse’s change rate 2.5 (%/s), the first natural frequency varies within the range of (2.32–2.47) (Hz), and the second varies within the range of (7.95–8.51) (Hz). At the masse’s change rate 5 (%/s), the two natural frequencies vary within the range of (2.19–2.47) (Hz) and (7.49–8.51) (Hz), respectively.

Figure 5.

Modal parameters of the system.

The displacement responses of the system under a particular excitation were obtained through the ODE45 method in MATLAB using a fixed-integration step and recorded at a sampling frequency of 100 (Hz). At the mass change rate 2.5 (%/s), Figure 6(a) plots the non-stationary vibration displacement signal. The signal’s spectrogram under a random excitation is presented Figure 6(b). Figure 7(a) and (b) show the non-stationary vibration displacement signal and the signal’s spectrogram under a random excitation at the mass change rate 5 (%/s).

Figure 6.

(a) Realization of the non-stationary vibration displacement and (b) Short-time Fourier transform of the signal with the gradual change at 2.5 (%/s).

Figure 7.

(a) Realization of the non-stationary vibration displacement and (b) short-time Fourier transform of the signal with the gradual change at 5 (%/s).

The proposed method uses the AR model to identify the mechanical system’s modal parameters in each data segment. The model’s parameters at the next segment are updated using the results of previous segments. The length of each window is sufficient to contain the modal parameters, and the length must be at least four times the highest period in order to track the changes of all modes. The model order selection procedure is based on the minimum description length (MDL) (Bui et al., 2022).

Figure 8 presents the identification of the modal parameters with the proposed method at the mass change rate 2.5 (%/s). Obviously, the natural frequencies are accurately determined and tracked. The results of damping ratios for the first and second modes are also in very good agreement with the calculated values. The proposed method performs better in terms of accuracy resolution than STFT. It can be observed from these figures that the proposed method matches well with theoretical variations. At the higher change rate of the masses, the modal parameters of the system using the proposed method are shown in Figure 9.

Figure 8.

Modal parameters of the system with the gradual change at 2.5 (%/s) using the proposed method.

Figure 9.

Modal parameters of the system with the gradual change at 5 (%/s) using the proposed method.

Table 1 presents a comparison of the computation time for the modal parameter diagrams of the proposed method with the recently published updated method by Schur complement (Bui et al., 2022). Generally, the identified results by the proposed method match with the simulated results presented in Figure 5. It can be seen that the first mode is still accurately identified and tracked. However, the second mode is largely dispersed, especially its damping ratio, which is known to associate with a greater uncertainty (Vu et al., 2011a). Conservatively, the proposed method can identify and monitor the mass change rate at 5 (%/s). This change rate is quite high for real systems. As can be seen from Table 1, the proposed method performs better than the Schur method.

Table 1.

Comparison of the computational time of two methods: method by Schur complement (Bui et al., 2022) and proposed method for modal parameter diagram (data from two-channels).

	Method by Schur complement (Bui et al., 2022)		Proposed method
Mass change	2.5 (%/s)	5 (%/s)	2.5 (%/s)	5 (%/s)
Time (s)	114	114	103	103

5.2. Experimental data

In this section, the proposed method is applied to a hydraulic turbine blade to monitor the modal parameters. The turbine blade is made of bronze alloy “M”-C92300 corresponding to the standard designation 87Cu-8Sn-IPb-4Zn.F. Figure 10 shows the configuration of the test in which four accelerometers have been mounted to record the accelerations of this blade.

Figure 10.

Blade with four accelerometers.

5.2.1. The turbine blade in air

The test was carried out using an LMS system. A PCB impact hammer with a sampling frequency of 6400 (Hz) acted on the structure in the static test (Figure 11). The natural frequencies of this blade shown in Table 2 are obtained by different methods: Ansys, the power spectral density (PSD) of responses in MATLAB, short-time autoregressive model (STAR) (Vu et al., 2011a), method by Schur complement (Bui et al., 2022), and the proposed method. As seen, the proposed method reveals excellent matching results with previous identification methods.

Figure 11.

Modal test of the blade in the air.

Table 2.

Modal identification of the blade in the air (Hz).

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
Ansys	Non-identify	332.9	643.16	896.9	963.6	1157.6
PSD in MATLAB	204.7	359.4	582.8	843.8	1050	1295.0
STAR (Vu et al., 2011a)	208.0	368.0	576.0	848.0	1040	1294.5
Method by Schur complement (Bui et al., 2022)	204.0	357.2	579.4	845.9	1047	1293.0
Proposed method	204.4	360.3	582.5	846.4	1055	1297.0

The first frequency (204.7 Hz) was not identified using the finite element method because it was the mounting-structure mode. The system’s PSD plot obtained from the spectrogram function in MATLAB is presented in Figure 12. The stabilization diagram obtained from the proposed method is shown in Figure 13. As shown in Table 2, the identifications of frequencies by the proposed method are pretty accurate for all modal parameters. The proposed method with order updating was applied to the systems to extract the natural frequencies.

Figure 12.

Spectra of the blade in the air by spectrogram functions in MATLAB.

Figure 13.

Order-updating stabilization diagram of the blade in the air using the proposed method.

5.2.2. Turbine blade raised from water

5.2.2.1. Experiment setup

A hydraulic test bench was constructed in the research laboratory to evaluate fluid-structure interactions. When the blade was placed inside a perforated tank, the vibration of the blade in air and water could be measured. The valves and various outlet nozzles are used to control the flow speed. The flow rates and therefore the depth of the blade in the water can be adjusted according to the desired flow velocities. A fastening system is a multi-tasking tool for mounting a blade in variable boundary conditions.

The PCB330A sensors are utilized for vibration measurements. The Vishay System 6000 was used to record pressure sensors during different tests. The data is exported and saved in different formats, such as.xls and .txt, by the software “Strain Smart” designed to work with the Vishay acquisition box. A PCB impact hammer is acted on the structure tests. The hammer is equipped with a steel extension to hit the plate at different depths. The configuration of the experiment is depicted in Figure 14. The blade is submerged at different depths, and the depth/length ratio varies from 0.4 (totally submerged) to 0 (in air).

Figure 14.

Modal test of the blade in water.

The natural frequencies of the blade change with respect to the submerged depth due to the effect of the fluid. Before the blade rises, its modal parameters were determined by analytical and experimental methods at different depth length ratios (D/L). The result is shown in Table 3. The temporal response data at sampling frequency 8192 (Hz) are depicted in Figure 15(a) and (b) plots the short-time Fourier transform of the signals.

Table 3.

Modal identification of submerged blade by impact test.

Mode	Frequencies (Hz) of submerged condition (depth/blade length)
	0.4 (totally submerged)	0.2	0.1	0 (totally in air)
1st	160.8	166.85	170.8	Non-identify
2nd	288.1	288.1	290	332.9
3rd	352.7	354	358	643.2
4th	496.9	498	500	896.9
5th	640.5	642.8	645	963.6

Figure 15.

(a) Response to the blade's acceleration and (b) short-time Fourier transform of the signal.

The modal signal-to-noise ratio (Vu et al., 2007) is utilized to separate spurious and real modes. The smallest frequency is then obtained for the next window computation. In this section, the proposed method was applied to track modal parameters for the submerged blade. The length of each data segment was chosen to track all the modes of the signals. The size length of the sliding window varies and is chosen to be at least 4 times that of the longest natural period of the previous block (Bui et al., 2022).

Before the blade rises, its modal parameters was determined by analytical and experimental methods at different depth length ratios (D/L). The result is shown in Table 3.

In this section, the proposed method was applied to the measured signals. The length of each data segment was chosen to track all the modes of the signals. In each window, the minimum description length (MDL) was used to obtain the order of the models. The system's natural frequencies were extracted using the proposed method as shown in Figure 16.

Figure 16.

Natural frequencies of the blade using the proposed method.

Compared with the modal testing in Table 3 and STTF of the signals shown in Figure 15, the identifications of frequencies (147–151.6, 191.2–199, 292.2–303, 435.3–451.4, 579.6–605.2) (Hz) are quite accurate for all modal parameters. Figure 16 shows that the natural frequencies increase slightly when the blade rises from the water to the surface under turbulence. This increase is relevant to the depth ratio of the blade submerged in water. This result agrees well with the conclusion in other authors (Vu et al., 2007).

During the rise-up, the blade is in a moving condition, and the boundary condition is affected by the movement. The clamp and the fixing systems are also in vibration. These lead to a little discrepancy to results in Table 3 which are from the modal testing (static), Figure 16 is dynamic turbulent testing, so some other effects, such as the boundary condition, change of mass, and stiffness, can explain the discrepancy. However, the two representations agree on the frequencies and the variation of the frequencies.

7. Conclusion

In this work, a new method of the modal parameter identification and monitoring using a vector autoregressive model has been introduced for non-stationary vibrations. The technique uses SVD and a short-time sliding window as autoregressive models to extract the model parameters. The method aims to reduce the computational complexity for the online slow-varying monitoring of non-stationary cases by updating the SVD to estimate AR model parameters through the order and time of the previous computational window. The proposed method was validated first through a numerical simulation of a mechanical system at the different rates of masses and then through experiments on a submerged hydraulic turbine blade. Results show that the proposed method is a powerful technique for analysis and monitoring the modal parameters in non-stationary vibration systems under a reasonable varying rate at 5% per second. Natural frequencies can be accurately obtained and monitored. The identification and monitoring of damping ratios are sufficiently convincing.