Operational modal analysis of three-dimensional structures by second-order blind identification and least square generalized inverse

2.1. The model of blind source separation

Supposing that observed signals $\mathbf{X}\left(t\right)=\left[{\overrightarrow{x}}_1\right(t),{\overrightarrow{x}}_2(t),\cdots ,{\overrightarrow{x}}_n(t){]}^T\in {\mathbb{R}}^{n\times T}$ mixed by several unobserved signals $\mathbf{S}\left(t\right)=\left[{\overrightarrow{s}}_1\right(t),{\overrightarrow{s}}_2(t),\cdots ,{\overrightarrow{s}}_m(t){]}^T\in {\mathbb{R}}^{m\times T}$, the linear blind source separation model under the condition of non-noise is expressed as:

where, $\mathbf{A}\in {\mathbb{R}}^{n\times m}$ is a mixing matrix. In the present of measurement noises, the source identification model can be expressed as:

In order to recover source signals $\mathbf{S}\left(t\right)$ from observed mixing signals $\mathbf{X}\left(t\right)$, it is necessary to find a separated matrix $\mathbf{W}\in {\mathbb{R}}^{m\times n}$ to meet the following equation:

where, $\mathbf{Y}\left(t\right)$ is an estimation matrix of $\mathbf{S}\left(t\right)$, Fig. 1 shows the principle of blind source separation.

Fig. 1The model of blind source separation

2.2. Using SOBI to solve blind source separation model

The main idea of the blind source separation is that recovering the unobserved source signals from multiple observed mixing signals without any prior information. SOBI is a blind source separation technique, which is different from independent component analysis [16]. By using joint diagonalization of a set of covariance matrices, SOBI can estimate the separated matrix effectively [10].

Second order blind identification mainly uses the joint diagonalization of a set of covariance matrices to estimate source signals accurately. There are two hypotheses: The mixing matrix is a column full rank matrix and the sources are not related to each other.

According to the hypotheses, the whiten matrix $\mathbf{V}\in {\mathbb{R}}^{n\times n}$ is calculated by eigenvalue decomposition as below:

where $\mathbf{D}\mathbf{i}\mathbf{a}\mathbf{g}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(d_1,\cdots ,d_n)\in {\mathbb{R}}^{n\times n}$ is a diagonal matrix composed by eigenvalues $d_1,\cdots ,d_n$, $\mathbf{E}\mathbf{i}\mathbf{g}\mathbf{e}\mathbf{n}\in {\mathbb{R}}^{n\times n}$ is a eigenvector matrix of the covariance matrix ${\mathbf{C}}_{\mathbf{X}\mathbf{X}}=E\left\{\mathbf{X}\right(t\left){\mathbf{X}}^T\right(t\left)\right\}\in {\mathbb{R}}^{n\times n}$. Thus, the observed matrix $\mathbf{X}\left(t\right)$ is preprocessed to $\mathbf{Z}\left(t\right)\in {\mathbb{R}}^{n\times T}$:

When the matrix $\mathbf{V}$ is a whiten matrix, the covariance matrix of $\mathbf{Z}\left(t\right)$ is a unit matrix, and the second-order correlation between each component has been removed.

The set of covariance matrices is estimated as follows:

where $\tau \in \left\{{\tau }_j\right|j=1,2,\cdots ,k\}$ is a fixed time delay.

For all the $\mathbf{R}\left({\tau }_j\right)\in {\mathbb{R}}^{n\times n}$, using joint diagonalization algorithm, the orthogonal matrix $\mathbf{U}\in {\mathbb{R}}^{n\times m}$ is estimated:

where ${\mathrm{\Lambda }}_j\in {\mathbb{R}}^{m\times m}$ is a diagonal matrix.

Therefore, the source signals can be estimated as $\mathbf{Y}\left(t\right)$:

And the separated matrix $\mathbf{W}$ is:

An estimation matrix of mixing matrix $\mathbf{A}$ is:

2.3. Uncertainty factors in using SOBI to solve blind source separation model

It is easy to find some inevitable ambiguity or uncertainty factors from Eq. (1), which are summarized as followings:

a) The variance of the separated signal is uncertain.

Suppose $\mathbf{D}$ is a diagonal matrix, the linear instantaneous mixture modal can be expressed as:

So, the amplitude of separated signal is inconsistent with the amplitude of source signals.

b) The order of the separated signal is uncertain.

Suppose $\mathbf{P}$ is a permutation matrix, the linear instantaneous mixture modal can be expressed as:

where $\mathbf{P}\mathbf{S}\left(t\right)$ is a new source signals matrix after reordering and $\mathbf{A}{\mathbf{P}}^{-1}$ is a new mixing matrix.

The reason is that it is impossible to determinate specific values of the mixing matrix $\mathbf{A}$ and source signals $\mathbf{S}\left(t\right)$ simultaneously without any prior knowledge.

c) The number of separated components is uncertain

It is hard to identify and separate independent components, if the contribution of independent components is not sufficient. So, without any prior knowledge, to determine the number of separated components or sources by the SOBI algorithm is difficult.

2.4. Modal coordinates decomposition of one-dimensional structural vibration response signal

OMA extracts modal parameters (including mode shapes, natural frequencies and damping ratios) only from response signals when the structures at working condition. For $n$-degree of freedom time-invariant vibration systems, the kinetic equation is:

In Eq. (13), $\mathbf{M}\in {\mathbb{R}}^{n\times n}$, $\mathbf{C}\in {\mathbb{R}}^{n\times n}$, $\mathbf{K}\in {\mathbb{R}}^{n\times n}$ represent the mass matrix, the damping matrix and the stiffness matrix, respectively. $\mathbf{F}\left(t\right)\in {\mathbb{R}}^{n\times T}$ is an external excitation, as well as $\mathbf{X}\left(t\right)\in {\mathbb{R}}^{n\times T}$ is a displacement matrix, $\dot{\mathbf{X}}\left(t\right)$ and $\ddot{\mathbf{X}}\left(t\right)$ are the first derivative and the second derivative of $\mathbf{X}\left(t\right)$ respectively.

For proportional damped vibration systems under ambient excitation, it’s random vibration response that can be decomposed in modal coordinates:

In Eq. (14), $\mathbf{\Phi }\in R^{n\times m}$ represents modal shape matrix making a series of mode shapes ${\overrightarrow{\varphi }}_i$ and $\mathbf{Q}\left(\mathbf{t}\right)\in {\mathbb{R}}^{m\times T}$ is the vector matrix composed of each order modal response ${\overrightarrow{q}}_i\left(t\right)$. Through the modal response, the mode frequency ${\omega }_i$ and modal damping ratio ${\xi }_i$ can be calculated.

By the way, normalized modal shapes and the modal response matrix meet:

where $E\{\cdot \}$ represents expectation, and ${\mathrm{\Lambda }}_{n\times n}$ is a diagonal matrix whose size is $n\times n$.

2.5. OMA for one-dimensional structures by SOBI and its uncertainty factors

According to the mode theory, it is well known that each of the modal coordinates vector is independent. Combining with the SOBI algorithm, the response of structures decomposes into mode coordinate vectors and a mode shape matrix. Modal coordinates vectors can be regarded as source signals, and responses of systems can be recognized as observed signals, as well as the modal shape matrix is the mixing matrix. Therefore, SOBI can identify the modal parameters, and the process of identifying as showed in Fig. 2.

Fig. 2The process of identifying modal parameters by SOBI

On account of three assumptions of the SOBI method, the identified modal parameters have the following characteristics.

Every order of the modal shape has the different amplitude. In case of the SOBI method, the energy of separation matrix is not unique, and separated components also lose amplitude information. Unlike the principal component analysis method [17], the SOBI method cannot get the contribution ratio information of each modal. The modal shape is a relative quantity rather than an absolute value. So, in order to compare the modal shape with the real modal shape, the separation matrix and the modal shape identified by the SOBI method should be normalized.

The order of modal parameter identified by SOBI is uncertain. Modals identified by the SOBI method are not in accordance with the order from small to large. In fact, the first separated output modal is the one whose independence is the strongest rather than the first order modal parameter. Therefore, in order to compare natural frequencies with real natural frequencies, modal parameters need to be reordered by modal frequencies.

It is hard to identify and separate independent components, and to determine the numbers of modal parameters if the contribution of independent components is not sufficient. So, in the SOBI algorithm based output-only modal analysis, if the modal parameter is a small contribution to independence, the SOBI method may barely identify them and will cause modal parameters missing.

3.1. Problems and differences in OMA for one-dimensional to three-dimensional structures

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.

As the three-dimensional structures are more complex, it is inevitable that the errors of the identified modal parameters will be the greater than the identified modal parameters in one-dimensional structures.

3.2. Modal coordinates decomposition of three-dimensional structural vibration response signals

For complex three-dimensional continuous systems under ambient excitation, its random response in time domain can be decomposed in modal coordinates:

where ${\overrightarrow{\psi }}_i$, ${\overrightarrow{o}}_i$, ${\overrightarrow{b}}_i$ are the $i$th ($i=1,2,\dots ,\infty$) mode shape vector of ${\mathbf{X}}_{three}$,${\mathbf{Y}}_{three}$and ${\mathbf{Z}}_{three}$ direction. ${\overrightarrow{h}}_i\left(t\right)$ is the $i$th modal coordinate response which contains the information of the modal natural frequency.

In finite element analysis and experimental analysis, the three-dimensional systems are discretized to $d$ test points, and the time is discretized to $T$ number of sampling points, and $t=1,2,\dots ,T$. Therefore, the vibration response in the domain of time can be approximated as:

where $l$ is the number of modal truncation calculated by finite element analysis. ${\mathbf{H}}_{l\times T}\left(t\right)$ of three directions are the same, so the modal coordinate response just should be calculated only once.

3.3. Right pseudo inverse and its general reversion the minimal norm solution of least square

Modal coordinates responses ${\mathbf{H}}_{l\times T}\left(t\right)$ are the same in three directions. $T$ is the discretized number, and $l$ is the the number of modal truncation. So $T>l$ and $rank\left({\mathbf{H}}_{l\times T}\right(t\left)\right)=l$, the matrix ${\mathbf{H}}_{l\times T}\left(t\right){\mathbf{H}}_{T\times l}^T\left(t\right)\in {\mathbb{R}}^{l\times l}$is reversible. The right pseudo inverse matrix defined as [19]:

where the right pseudo inverse matrix meets ${\mathbf{H}}_{l\times T}\left(t\right){\mathbf{H}}_{T\times l}^{+}\left(t\right)={\mathbf{I}}_{l\times l}$, the right pseudo inverse matrix is uniquely determined and right pseudo inverse is associated with its general reversion the minimal norm solution of least square.

After decomposing the response in one direction based on the SOBI modal analysis algorithm, this study requires to take use of the mode shape of the first direction, and calculates the modal parameters of the other two directions by the general reversion of least square algorithm.

The least squares generalized inverse method is the unbiased optimal estimate of the sum of squares of errors and the minimum sense [19, 20]. Supposing the vibration response $X$ direction is the largest, the mode shape matrix ${\mathbf{\Psi }}_{d\times l}$ and the modal coordinate response ${\mathbf{H}}_{l\times T}\left(t\right)$ of $X$ direction is known, according to Eq. (18), the vibration response of $Y$ direction:

Multiplied by ${\mathbf{H}}_{T\times l}^T\left(t\right)$ on both sides of Eq. (20):

So, the mode shape matrix ${\mathbf{O}}_{d\times l}$ of $Y$ direction can be expressed as:

In the same way, the mode shape matrix of $Z$ direction also can be got:

3.4. OMA for three-dimensional structures by SOBI and general reversion of least square

The select of a direction decomposed by the SOBI algorithm: in the first place, one direction of three-dimensional structures is decomposed by SOBI, however if selecting the different directions to decompose, it will affect the final effect of modal parameters. So, which directions are decomposed needs to be explained theoretically and demonstrated by experiments.

The modal response matrix returns into other two directions: since the modal response is the same in three directions, decomposing the response of one direction based on SOBI modal analysis algorithm which requires to return into the other two directions. The modal response is a matrix. How to return into the other directions? In practice, the response of other two directions is multiplied by the inverse of the modal response matrix.

In order to identify operational modal parameters of three-dimensional structures, a strategy is adopted as below and the process of OMA for three-dimensional structures is shown in Fig. 3 below (Supposing that vibration response of $X$ direction is largest).

1) The purposed method takes advantages of SOBI algorithm to decompose a direction of the response signals of displacement which vibration response is the largest. Supposing the vibration response $X$ direction is the largest, we can get the mode shape matrix ${\mathbf{\Psi }}_{d\times l}$ and modal coordinate response ${\mathbf{H}}_{l\times T}\left(t\right)$ of $X$ direction.

2) The purposed method calculates the mode shape matrix of other two directions by using the general reversion of least square algorithm. The mode shape matrix ${\mathbf{O}}_{d\times l}$ of $Y$ direction can be calculated with Eq. (22) and the mode shape matrix ${\mathbf{B}}_{d\times l}$ of $Z$ direction can be calculated with Eq. (23).

3) The modal parameters of three-dimensional structures are identified by the matrix assembly method. In other words, $[{\overrightarrow{\psi }}_i,{\overrightarrow{o}}_i,{\overrightarrow{b}}_i]$ is the $i$th mode shape of three-dimensional structures, and ${\overrightarrow{h}}_i\left(t\right)$ is the $i$th modal coordinates response.

Fig. 3Process of OMA for three-dimensional structure

4.1. Simulation data generation of three-dimensional structures

The simulation data mainly pays attention to the cylindrical shell. The parameters settings of cylindrical shell are listed as follows: the thickness is 0.005 m, the length is 0.37 m, the radius is 0.1825 m, the elasticity modulus is 205 GPa, the poison’s ratio of materials is 0.3, the density of materials is 7850 kg/m³, mode damping ratios are 0.03 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. 4. 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 [18] 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. 5.

Response data are divided into two groups: data without measurement noises and data with 5 % Gauss measurement noises.

Fig. 43-D picture of the cylindrical shell

Fig. 5The 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

In order to evaluate the effect of identification about the new method of three-dimensional structures, the mode shapes and natural frequencies are calculated by the finite element analysis (FEA) method as the real modal parameters to compare with the identified modal parameters. Modal assurance criterion (MAC) is an important criterion to reflect the effectiveness of the modal identification by the new method. The modal assurance criterion is defined as:

In Eq. (24), ${\overrightarrow{\varphi }}_j$ is a theoretical value by the FEA method of the $i$th order mode shape and ${\overrightarrow{\stackrel{-}{\varphi }}}_j$ represents the identified value by the new method of $i$th order mode shape. It is worth noting that the value of MAC ranges from 0 to 1, and the value closer to 1 states the higher accuracy of the identified mode shape.

4.3. Simulation verification results

From Fig. 5, the vibration response of the $X$ direction and the $Y$ direction is larger than the $Z$ direction, so according to Fig. 3 about the process of OMA for three-dimensional structures. At first the $X$ direction is decomposed by SOBI method. Then with the aid of the general reversion of least square algorithm, the other two directions are calculated. At last, identify the modal parameters of three-dimensional structure are identified by matrix assembly method. At the same time, another method is introduced to assemble three directions directly and then identify the modal parameters of three-dimensional structures directly is compared with the proposed method.

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

Fig. 6Real 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. 7. Table 1 and Table 2 are the identified frequency by the purposed algorithm comparison of the real frequency by FEA and the MAC of mode shapes respectively.

According to the identified modal shapes by the SOBI method, Table 3 shows the MAC values between each identified modal shape. Each modal shape of its self MAC value is 1, and one modal shape of MAC with other modal shapes are small and almost close to 0.

When the damping ratio of the cylindrical shell is 0.03 and with 5 % measurement noises, the mode frequency is identified by the purposed algorithm that is shown in Fig. 8 below. Table 4 and Table 5 are the identified frequency by the purposed algorithm comparison of real frequency by the FEA method and the MAC of mode shapes respectively.

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

a) FFT of the 1th separated component

b) FFT of the 2th separated component

c) FFT of the 3th separated component

d) FFT of the 4th separated component

e) FFT of the 5th separated component

f) FFT of the 6th separated component

Fig. 8The identified mode frequencies when the damping ratio is 0.03 and with 5 % measurement noises

a) FFT of the 1th separated component

b) FFT of the 2th separated component

c) FFT of the 3th separated component

d) FFT of the 4th separated component

e) FFT of the 5th separated component

f) FFT of the 6th separated component

When the damping ratio of the cylindrical shell is 0.03, Table 6 shows the mode shapes identified by finite element analysis, the purposed algorithm without noises and the purposed algorithm with 5 % measurement noises.

When the damping ratio of the cylindrical shell is 0.1, the mode frequency identified by the purposed algorithm is shown in Fig. 9. Table 7 and Table 8 is the identified frequency by the purposed algorithm comparison of real frequency by FEA and the MAC of mode shapes respectively. And “–” represents that the mode parameters are not identified.

Table 6Mode shapes when the damping ratio is 0.03

Method	Mode shapes
	1th order	2th order	3th order	4th order	5th order	7th order
FEA
This method
This method (5 % noise)

Fig. 9The identified mode frequencies when the damping ratio is 0.1

a) FFT of the 1th separated component

b) FFT of the 2th separated component

c) FFT of the 3th separated component

d) FFT of the 4th separated component

e) FFT of the 5th separated component

f) FFT of the 6th separated component

When the damping ratio of the cylindrical shell is 0.1 and with 5 % measurement noises, the mode frequency identified by the purposed algorithm is shown in Fig. 10. Table 9 and Table 10 are the identified frequency by the purposed algorithm comparison of real frequency by FEA and the MAC of mode shapes respectively.

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

a) FFT of the 1th separated component

b) FFT of the 2th separated component

c) FFT of the 3th separated component

d) FFT of the 4th separated component

e) FFT of the 5th separated component

f) FFT of the 6th separated component

When the damping ratio of the cylindrical shell is 0.1, Table 11 shows the mode shapes identified by finite element method, the purposed algorithm without noises and the purposed algorithm with 5 % measurement noises.

Table 11mode shapes when the damping ratio is 0.1

Method	Mode shapes
	1th order	2th order	3th order	4th order	5th order	7th order
FEA
This method					–
This method (5 % noise)				–	–	–

4.4. Results analysis

1) From the result identified by the purposed method, Fig. 7 and Fig. 9 show the identified modal frequencies with different damping ratios, by comparing with the frequencies calculated by the FEA method in Table 1 and Table 7. The modal frequencies identification precision changes when damping ratios increases. However, the relative errors are under 5 %. Table 6 and Table 11 show the mode shapes under different methods.

2) From Table 1, Table 4, Table 7 and Table 9, the order of identified modal frequencies and real modal frequencies are not one-to-one correspondence. The reason is that the order of the mode coordinate vector is uncertainty according to Eq. (14).

3) From all the results, the 6th order modal parameters are not identified because of its small modal contribution in dynamic response measurement signals.

4) With the increasing of damping ratio, the identified mode shapes are becoming deformed. It is more difficult to identify the modal parameters because the larger of damping ratio, the smaller of response of the structure.

5) Contrasting Fig. 7 with Fig. 8, and from Table 4 and Table 9. In conclusion, this method is robustness to Gauss measurement noises disturbances. Because second order statistics rather than higher-order statistics method is used in SOBI method.