Three-dimensional operational modal analysis based on self-iteration principal component extraction and direct matrix assembly

3.1. Self-iteration principal component extraction algorithm

The SIPCE algorithm is the improvement of (Nor-linear Iterative Partial Least Squares) NIPALS algorithm [17]. NIPALS gets mutually orthogonal eigenvectors by projecting the high-dimensional data samples to the low-dimensional feature space, then establishing a linear regression relation between eigenvectors [18]. NIPALS model has less Variables, it’s more robust to noise. Since the birth of the NIPALS algorithm, it’s wildly used to calculate PCs, the calculation process is actually equivalent to EVD.

The following are algorithmic steps, a column from the physical responses matrix $\mathbf{X}\left(t\right)\in {\mathbb{R}}^{m\times T}$, $t=$ 1, 2, …, $T$ is chosen. The first principal component vectors ${\overrightarrow{y}}_1\left(t\right)\in {\mathbb{R}}^{1\times T}$ is randomly chosen from signal matrix $\mathbf{X}\left(t\right)\in {\mathbb{R}}^{m\times T}$, $t=$ 1, 2, ..., $T$, and then the first transformation matrix vector ${\overrightarrow{w}}_1\in {\mathbb{R}}^{m\times 1}$ is iteratively calculated. As Eq. (17) shows, $\mathbf{E}\left(t\right)$ denotes error matrix, error matrix is updated to calculate the second principal component vectors ${\overrightarrow{y}}_2\left(t\right)\in {\mathbb{R}}^{1\times T}$ and the second transformation matrix vector ${\overrightarrow{w}}_2\in {\mathbb{R}}^{m\times 1}$ and so on:

Initially setting max iterations times $J_{\mathrm{m}\mathrm{a}\mathrm{x}}$, accuracy threshold $\alpha$, contribution threshold $\eta$. Assume index $i$ denotes the $i$th principal moment, and $i=$ 1. Setting original error matrix ${\mathrm{{\rm E}}}_i\left(t\right)=\mathbf{X}\left(t\right)$.

Chose the arbitrary column data vector from signal matrix, assume the chosen column is the $l$th column, so vector ${\overrightarrow{x}}_l\left(t\right)\in {\mathbb{R}}^{1\times T}$ is chosen from signal matrix $\mathbf{X}\left(t\right)\in {\mathbb{R}}^{m\times T}$, initializing $j=$ 1.

Calculate ${\overrightarrow{w}}_i^{\left(j\right)}\in {\mathbb{R}}^{m\times 1}$: ${\overrightarrow{w}}_i^{\left(j\right)}=\left(\mathbf{X}\right(t\left){{\overrightarrow{y}}_i^{\left(j\right)}}^T\right(t\left)\right)/\left({\overrightarrow{y}}_i^{\left(j\right)}\right(t\left){{\overrightarrow{y}}_i^{\left(j\right)}}^T\right(t\left)\right)$;

Normalize the vector ${\overrightarrow{w}}_i^{\left(j\right)}$: ${\overrightarrow{w}}_i^{\left(j\right)}={\overrightarrow{w}}_i^{\left(j\right)}/‖{\overrightarrow{w}}_i^{\left(j\right)}‖$;

Update ${\overrightarrow{y}}_i^{(j+1)}\left(t\right)$: ${\overrightarrow{y}}_i^{(j+1)}\left(t\right)={{\overrightarrow{w}}_i^{\left(j\right)}}^T\mathbf{X}\left(t\right)/\left({{\overrightarrow{w}}_i^{\left(j\right)}}^T{\overrightarrow{w}}_i^{\left(j\right)}\right)$;

Compare the ${\overrightarrow{y}}_i^{(j+1)}\left(t\right)$ in the step 4) with the ${\overrightarrow{y}}_i^{\left(j\right)}$ in the step 2), define variable ${\beta }_i^{\left(j\right)}=\frac{‖{\overrightarrow{y}}_i^{(j+1)}-{\overrightarrow{y}}_i^{\left(j\right)}‖}{‖{\overrightarrow{y}}_i^{\left(j\right)}‖}$, if ${\beta }_i^{\left(j\right)}\le \alpha$, go to the step 7); else, $j=j+$ 1 and return to the step 2).

Calculate ${\lambda }_i$: ${\lambda }_i=‖\frac{1}{T}\mathbf{X}\left(t\right){\mathbf{X}}^{Tr}\left(t\right){\overrightarrow{w}}_i‖/‖{\overrightarrow{w}}_i‖$.

Update the error matrix ${\mathbf{E}}_i\left(t\right)$:

Define the contribution of current mode${\epsilon }_i=\frac{{\lambda }_i}{\sum _{s=1}^i{\lambda }_s}$, if ${\epsilon }_i>\eta$ and $j<J_{\mathrm{m}\mathrm{a}\mathrm{x}}$, $i=i+$ 1 and return to the step 2), else go to step 10)

Get the PCs matrix $\mathbf{Y}\left(t\right)=\left[{\overrightarrow{y}}_1\right(t),{\overrightarrow{y}}_2(t),\dots ,{\overrightarrow{y}}_i(t),\dots ,{\overrightarrow{y}}_n(t){]}^{Tr}\in {\mathbb{R}}^{n\times T}$ and the linear transformation matrix $\mathbf{W}=[{\overrightarrow{w}}_1,{\overrightarrow{w}}_2,\cdots ,{\overrightarrow{w}}_n]\in {\mathbb{R}}^{m\times n}$.

The step 9) is to calculate the oversized contribution rates ${\epsilon }_i$ correspond to the $i$th principal components. So, ${\lambda }_i$ is the approximate contribution rate, which corresponds to the $i$th principal component. Principal components are extracted by the order of cumulative variance proportion. Contribution rates ${\epsilon }_i$ is the proportional size of the current principal component in all calculated principal components. $\eta$ is contribution threshold, if ${\epsilon }_i<\eta$, stopping iteration, which means the current principal component is not important enough. And the principal components extraction process is complete. Fig. (3) is the flowchart of self-iteration principal component extraction algorithm based on OMA.

Fig. 3 shows the flowchart of self-iteration principal component extraction based OMA algorithm. Variable ${\lambda }_i(i=$ 1, 2, …, $m)$ represents the contribution of its corresponding principal components, ${\epsilon }_i$ is the importance index of next principal components, and $\eta$ is a setting threshold value which is used to control the iterations. Therefore, it realizes the principal components extraction.

Difference from calculating all modal orders at one time in traditional batch-processing PCA method, SIPCE just extracts several modal orders according to practical need. When there are many response measurement points and sample times, the SIPCE method has obvious advantages in time and space overhead. Besides that, SIPCE calculate modal parameters by iteration, which avoids singular values and ill-posed problem in traditional PCA method, so it always has higher identification accuracy. Table 1 is the comparison between SIPCE and traditional PCA.

Fig. 3The flowchart of self-iteration principal component extraction algorithm based OMA

3.2. Three-dimensional OMA and its differences

Actual engineering structures are three-dimensional. From one-dimensional to three-dimensional, it takes a big step from scientific researches to engineering applications of modal analysis based on a series of blind source separation methods.

The matrix assembles: the vibration matrix of one-dimensional structures is a one-dimensional vector, the proportion of the value of each modal makes sense, but the amplitude of the modal does not make any sense. And the vibration matrix of the three-dimensional structure is combined by three directions $X$, $Y$ and $Z$. Each direction of the scale factor and the modulus ratio should stay the same, otherwise it becomes a one-dimensional vibration mode of three directions, rather than a three-dimensional vibration mode of structures. How to carry out this process is a big challenge.

For continuous three-dimensional structures, in finite element analysis and experimental analysis, the three-dimensional system is discretization to D test points, and the time is discretization to T number of sampling points. The displacement response in time domain can be expressed by modal coordinate approximation [13], both the speed signal and the acceleration signal are also suitable to Eq. (20):

where ${\mathbf{X}}^{three}\left(t\right)\in {\mathbb{R}}^{D\times T}$, ${\mathbf{Y}}^{three}\left(t\right)\in {\mathbb{R}}^{D\times T},{\mathbf{Z}}^{three}\left(t\right)\in {\mathbb{R}}^{D\times T}$ are physical responses in three directions, and $\mathbf{U}\in {\mathbb{R}}^{D\times D}$, $\mathbf{V}\in {\mathbb{R}}^{D\times D}$, $\mathbf{O}\in {\mathbb{R}}^{D\times D}$ are mode shapes matrixes of $X$, $Y$ and $Z$ directions. Matrix $\mathbf{Q}\left(t\right)\in {\mathbb{R}}^{D\times T}$ is the modal response which contains the information of modal natural frequency. Theoretically infers $\mathbf{Q}\left(t\right)$ of three directions are the same.

In three-dimensional modal parameters analysis, modal shape identification is the most important and difficult. Modal shape matrix of one-dimensional structure is a one-dimensional vector, the proportion of each element is crucial. However, the size of element has no practical significance. For modal shape matrix of three-dimensional structure, the three one-dimensional modal shape vectors that in $X$, $Y$, $Z$ direction are arranged in sequence and the proportion of each direction and element must remain fixed. If the proportions are changed, the obtained modal shapes are no longer three-dimensional modal shapes. So, it’s a big change to identify modal shape for three-dimensional structures.

3.3. The Assembled three-dimensional modal shapes

According to the previous method [13], the three directions’ modal coordinate responses are same, so the algorithm decomposes modal response signals in one direction which has the maximum vibration. Then through substitution of $\mathbf{Q}\left(t\right)\in {\mathbb{R}}^{D\times T}$ into the other two-dimension displacement response signals using Moore-Penrose matrix inverse in Eq. (21), modal shapes of these two dimensional can be identified. And three dimensional modal shapes $\left[\mathbf{U}\mathbf{V}\mathbf{O}\right]$ are assembled:

This paper proposes to assemble three directions’ modal response firstly, assume the assembled matrix is the whole response signal of three-dimensional structure. Then setting and solving the PCA model of the whole response signal, the modal coordinate responses and the three-dimensional modal shapes are calculated directly. The new displacement response in time domain can be freshly expressed by modal coordinate approximation:

where $\mathbf{A}\left(t\right)\in {\mathbb{R}}^{3D\times T}$ is the assembled modal signal matrix, $\left[\mathbf{U}\mathbf{V}\mathbf{O}\right]\in {\mathbb{R}}^{3D\times D}$ is the three-dimensional modal shapes and $\mathbf{S}\left(t\right)\in {\mathbb{R}}^{D\times T}$ is the modal coordinate response.

The method of least-square inverse matrix substitution firstly gets modal coordinate response of the largest response direction by matrix decomposition, then assume that is the global modal coordinate response to calculate modal shapes of the other two direction. However, the modal coordinate response of the largest response direction is not equal to the global modal coordinate response, there must be error in the final three-dimensional modal shapes. As for direct matrix assemble, it calculates the global modal coordinate response and three-dimensional modal shapes matrix directly by matrix decomposition, thus the new method avoids matrix inversion. In the calculation process, matrix inversion error and Ill-posed problems [20] are inevitable in the calculation, so the new method has higher identification accuracy than the matrix substitution method. Comparing Eq. (21) and Eq. (22), Table 2 is the performance comparison of two methods.

3.4. OMA for three-dimensional structure by SIPCE and Direct matrix assemble

In order to identify operational modal parameters of three-dimensional structure, a strategy is adopted as below and the process of OMA for three-dimensional structure is shown in Fig. 4 below.

Fig. 4Process of OMA for three-dimensional structure by SIPCE

The new method for three-dimensional structure replaces traditional PCA with SIPCE to solve the model. The new method sharply reduces space-time complexity, avoids singular values and ill-posed problem. In addition, the new method calculates modal shapes directly to avoid matrix inversion substitution, thus reducing the identify error. To summarize, SIPCE algorithm and the new modal assemble method can be well used in OMA for three-dimensional structures, and it has higher identification accuracy and lower time-space complexity than previous method.

4.1. Simulation data generation

The simulation data mainly pays attention to the cylindrical shell with two simple supported ends, and the excitation is Gauss white noise. The parameters setting of cylindrical shell is: the thickness is 0.005 m, length is 0.37 m, radius is 0.1825 m, elasticity modulus is 205 GPa, Poisson’s ratio of material is 0.3, density of material is 7850 kg/m³, and mode damping ratios are 0.03, 0.05 and 0.1 respectively.

The cylindrical shell is a continuum structure, and it must be discretized in order to calculate the modal and vibration response of the structure by finite element method. The more the number of discretized units, the more accurate the calculated modal and vibration response. Along the axial cylindrical shell distributes evenly 38 circles, each circle distributes evenly 115 observation spots which are shown in Fig. 5. Thus, the total of observation spots $d=$ 4370. The uniform white noises excitation is applied at each response point to generate the simulation data. In order to characterize the cylindrical shell structure of the complex three-dimensional mode shape with high accuracy, so many observation spots are chosen in the simulation. Operational modal parameter identification only requires that the number of response points is more than the number of the modal of the structure. So, for a real structure under normal working conditions, SOBI based three-dimensional operational modal parameter identification does not need so many sensors. The sensor points are more than or equal to the number of the modals in three directions of the structure. Of course, the more of the sensor points, the more accurate of the identification mode shapes. At the same time, the sampling frequency is placed at 5120 Hz, and the sampling time is set to 1 s, $T=$ 5120. The boundary conditions of the cylindrical shell are simply supported at both ends. At last, response signals of three directions are calculated by LMS Virtual.Lab [21] using finite element analysis (FEA). An observation spot is selected randomly, such as the 1118th spot whose response signals of three directions shown in Fig. 6 as below. Response data are divided into two groups: data without measurement noises and data with 5 % Gauss measurement noises.

Fig. 5Three-dimensional picture of the cylindrical shell

Fig. 6The response signals of three directions

a) Response of $X$ direction about the 1118th observation spot

b) Response of $Y$ direction about the 1118th observation spot

c) Response of $Z$ direction about the 1118th observation spot

4.2. Evaluation criteria and parameter setting

Theoretical values and FEA values are set as standard, which are calculated under the condition of undamped real modal shape and undamped real modal nature frequency, and PCA algorithm is based on EVD method.

In order to quantitatively compare the modal vibration mode, modal assurance criterion (MAC) is investigated [22], which is an important standard when considering the effectiveness of modal shape identification. It is computed by:

where ${\overrightarrow{\phi }}_i$ is the identified shape, ${\overrightarrow{\varphi }}_i$ represents the real shape, and (${\overrightarrow{\phi }}_i^T{\overrightarrow{\varphi }}_i$) represents the inner product of the two vectors. MAC values oscillate between 0 (no coincidence) and 1 (complete coincidence). And the more MAC value is close to 1, the more accurate estimated shape is.

1) Set precision threshold $\alpha =$ 0.000001;

2) Set maximum number of iteration $J_{\mathrm{m}\mathrm{a}\mathrm{x}}=$ 100;

3) Set truncation error $\eta =$ 0.01;

4)Computer configuration:

Operating system: 64-bit Windows 10,

Manufacturer: HP,

System Type: HP Z240 SFF Workstation,

Processor: Intel(R) Core(TM) i5-6500 CPU @ 3.20 GHz (4 CPUs),

Memory: 8192 MB RAM,

Programing language and environment: MATLAB 2014a.

4.3. Simulation verification results

The modal vibration modes and natural frequencies by the FEA method are used as the real value, show in Fig. 7.

Fig. 7Real mode shapes calculated by FEA

a) The 1th real modal shape

b) The 2th real modal shape

c) The 3th real modal shape

d) The 4th real modal shape

e) The 5th real modal shape

f) The 7th real modal shape

When the damping ratio of the cylindrical shell is 0.03, the mode frequency identified by the purposed algorithm is shown in Fig. 8 below. Table 3 is the identified frequency by the proposed algorithm comparison of real frequency by FEA. Fig. 9 is the identified modal shape, Table 4 is the MAC values of mode shapes.

When the damping ratio of the cylindrical shell is 0.03, and the 5 % measure noise is added into the response signals, the mode frequency identified by the purposed algorithm is shown in Fig. 10 below. Table 5 is the identified frequency by the purposed algorithm comparison of real frequency by FEA. Fig. 11 is the identified modal shapes. Table 6 is the MAC values of mode shapes.

Fig. 8The identified mode frequencies when the damping ratio is 0.03

a) FFT of the 1th order

b) FFT of the 2th order

c) FFT of the 3th order

d) FFT of the 4th order

e) FFT of the 5th order

f) FFT of the 6th order

Fig. 9The identified mode shapes when the damping ratio is 0.03

a) Modal shape of the 1th order

b) Modal shape of the 2th order

c) Modal shape of the 3th order

d) Modal shape of the 4th order

e) Modal shape of the 5th order

f) Modal shape of the 6th order

Fig. 10The identified mode frequencies when the damping ratio is 0.03 with 5 % measurement noise

a) FFT of the 1th order

b) FFT of the 2th order

c) FFT of the 3th order

d) FFT of the 4th order

e) FFT of the 5th order

f) FFT of the 6th order

Table 7 is the comparison of time cost and memory cost between the two methods.

Fig. 11The identified mode shapes when the damping ratio is 0.03 with 5 % measurement noise

a) Modal shape of the 1th order

b) Modal shape of the 2th order

c) Modal shape of the 3th order

d) Modal shape of the 4th order

e) Modal shape of the 5th order

f) Modal shape of the 6th order

4.4. Results analysis

1) Fig. 5-8 and Table 3-6 show the new method (SIPCE with direct matrix assemble) identifies the three-dimensional structures well, and it’s robust to interference noise, the 5 % Gauss measure noise has little affected to the identified accuracy.

2) For three-dimensional modal shapes, Table 3-6 show the new method (SIPCE+ Direct matrix assemble) has higher identified accuracy in four modes (3-7 modes), it’s more accuracy than the previous method (PCA+ Least-square matrix substitution) in general. And Table 6 emphasizes the method of direct matrix assemble is the crux of the high accuracy.

3) Table 7 shows this SIPCE sharply reduces the time and memory consumption. When the response signal matrix is assembled, the size of assembled signal matrix is much larger than before. So, this paper uses SIPCE to solve the PCA model of directly matrix assembly. The SIPCE’s advantage of low space-time complexity is obvious.

4) The 6th order modal parameters are missing because of its small modal contribution in structure’s dynamic response measurements.

5) In the simulation verification of cylindrical shell with two simple supported ends, there are two modal shapes with phase difference corresponding to one order natural frequency. Three-dimensional structure is more complex than one-dimensional structure, so there is always large error for identified three-dimensional modal shapes.

The authors declare that there is no conflict of interest regarding the publication of this article and regarding the funding that they have received.