Modal Parameter Identification of Structures Using Reconstructed Displacements and Stochastic Subspace Identification

Abstract

A method of modal parameter identification of structures using reconstructed displacements was proposed in the present research. The proposed method was developed based on the stochastic subspace identification (SSI) approach and used reconstructed displacements of measured accelerations as inputs. These reconstructed displacements suppressed the high-frequency component of measured acceleration data. Therefore, in comparison to the acceleration-based modal analysis, the operational modal analysis obtained more reliable and stable identification parameters from displacements regardless of the model order. However, due to the difficulty of displacement measurement, different types of noise interferences occurred when an acceleration sensor was used, causing a trend term drift error in the integral displacement. A moving average low-frequency attenuation frequency-domain integral was used to reconstruct displacements, and the moving time window was used in combination with the SSI method to identify the structural modal parameters. First, measured accelerations were used to estimate displacements. Due to the interference of noise and the influence of initial conditions, the integral displacement inevitably had a drift term. The moving average method was then used in combination with a filter to effectively eliminate the random fluctuation interference in measurement data and reduce the influence of random errors. Real displacement results of a structure were obtained through multiple smoothing, filtering, and integration. Finally, using reconstructed displacements as inputs, the improved SSI method was employed to identify the modal parameters of the structure.

1. Introduction

Modal analysis is an experimental method used for structural parameter identification (natural frequency, damping ratio, mode shape). This method is used in vibration response calculations, the root cause analysis of vibration problems, damage detection, and virtual-acoustic simulations. In addition, it is also used for adding flexibility to multibody analyses and speeding up structural durability. Therefore, modal analysis is very effective for the evaluation of structural changes to any types of responses [1,2].

The conventional modal analysis or the experimental modal analysis (EMA) requires a known excitation force on a system to obtain structural vibration response data, processes input and output signals to obtain the frequency response function, and finally, calculate modal parameters of the system by curve fitting [3]. This algorithm can better identify modal parameters of a system depending on its adaptability to the test environment; thus, it cannot be used under operating conditions. However, it is important to develop a method for the real-time and in situ acquisition of the state information of structures and the evaluation of their safety and residual life. In contrast, in the operational modal analysis (OMA), modal parameters of a structure are obtained only by measuring its output responses [4,5]. Orlowitz et al. [6] found no significant difference between the modal parameters of a Plexiglas plate obtained by OMA and EMA under similar boundary conditions.

Stochastic subspace identification (SSI), which was developed by Van Overschee and De Moor in 1991, is a powerful time-domain output-only modal parameter identification method [7,8,9,10,11]. SSI-based methods are less susceptible to prediction errors and computational efficiency [12]. Hence, continuous efforts have been made to improve the efficiency of SSI. Priori et al. [13] studied the identification of user-defined parameters in a pure output, data-driven random subspace through the asymmetric partition of the Hankel matrix. Wen et al. [14] proposed a sliding window-fuzzy C-means clustering algorithm in combination with deterministic-stochastic subspace identification (SC-CDSI), and reported that the proposed SC-CDSI identification algorithm could achieve the intelligent identification of online tracking of structural frequencies. Zhou et al. [15] improved the accuracy of modal parameter identification to determine the system matrix order of a damping ratio matrix. Gres et al. [16] estimated periodic information from measured data using a non-steady-state Kalman filter, and then removed the periodic information from the original output signal by an orthogonal projection. Tran et al. [17] proposed the combined use of different modal parameter identification methods, such as frequency domain decomposition, observer Kalman filter identification algorithm, and combined deterministic SSI.

Dynamic displacement-based modal parameter identification methods are more useful than acceleration-based modal parameter identification methods, especially for low-frequency structures [18,19,20]. However, acceleration-based modal parameter identification methods are widely used because of difficulties in displacement measurement and the appearance of the trend term interference in acceleration integration [21].

Therefore, it is still a key problem to reduce the interference of the trend term drift error and low-frequency noise errors in acceleration integration. Lee et al. [22] presented a displacement reconstruction method for low-frequency structures and introduced overlapping time windows to improve the reconstruction displacement accuracy and suppress low-frequency noises. However, when the target frequency was relatively low, it easily oscillated in the frequency domain. Hong et al. [23] used two finite impulse response (FIR) filters to reconstruct dynamic displacements. The accuracy of the proposed filter was investigated in detail using the transfer functions of the discrete filter for different displacement sampling rates and target accuracy. The results show that it was difficult to control the influence of low-frequency noise on the integration accuracy. Park et al. [24] used FEM-FIR filter de-noising to estimate the temporal and spatial derivatives of displacements in an elastic solid and the integral displacement for system identification. Brandt et al. [25] proposed a low-frequency cut-off algorithm to reconstruct displacements from measured accelerations and obtained the time-domain waveform of displacements according to the relationship between acceleration and displacement spectra.

In the present work, an improved strategy-based SSI method was proposed to identify the structural modal parameters using reconstructed displacements. The proposed strategy consisted of three steps. In the first step, displacements were estimated from measured accelerations. Due to the interference of noise and the influence of initial conditions, the integral displacement inevitably had a drift term. In the second step, the moving average method was employed to effectively eliminate random fluctuation interferences in measured data and reduce the influence of random errors. Real displacement results of a structure were obtained through multiple smoothing, filtering, and integration. Finally, the modal parameters of the structure were identified by the proposed SSI method using reconstructed displacements as input. The developed displacement-based SSI method is more stable and accurate for low-frequency-dominated structures.

2. Displacement Reconstruction Using Measured Acceleration Data

2.1. Mathematic Model for Acceleration Measurement

Interference terms, such as random interferences and direct current (DC) components, generally exist in tested acceleration signals. Random interferences vary with time, and DC components are the drifts of an accelerometer relative to the baseline under the influence of environmental factors [26]. Therefore, the measured acceleration can be expressed as:

$$a(t)=a_s(t)+a_n(t)+C$$

(1)

where $a_s(t)$ is the structural information of the measured acceleration, $a_n(t)$ is the random interference part of the measured acceleration, and C is the DC component of the measured acceleration.

The structural displacement signal can be obtained by the quadratic integration of the measured acceleration under the influence of the initial velocity and displacement:

$$x(t)=\int v(t)+x(0)={\displaystyle \iint a_s(t)dt}dt+{\displaystyle \iint a_n(t)dt}dt+\frac{1}{2}C{t}^2+\lambda t+v(0)t+\eta +x(0)$$

(2)

where $v(0)t$ is the first integral of the initial velocity, $x(0)$ is the initial displacement, $\frac{1}{2}C{t}^2$ is the integral term of DC components, and λ and η are constant terms of the first and second integrals, respectively.

The obtained displacement can be divided into three parts:

$$\tilde{x}(t)={\displaystyle \iint a_s(t)dtdt, \widehat{x}(t)={\displaystyle \iint a_n(t)dtdt}}, \alpha (t)=\frac{1}{2}C{t}^2+\lambda t+v(0)t+\eta +x(0)$$

(3)

where the first part represents the real displacement of the structure, and the second part represents the influence of environmental noise and measurement, which is small in the practical tests and can be removed by multiple smoothing and filtering. The third part represents the combined influence of DC components, initial conditions, and integral constant terms (it can be removed by subtracting the average value and fitting the polynomial by the least square method). Therefore, the main goal of the following analysis is the removal of the second and third parts of Equation (3).

2.2. Elimination of Trend Items and De-Noising

The least square method of polynomial fitting is most commonly used to eliminate the interference of trend terms [19]. In this work, a moving average method of least squares was used in combination with a filter to eliminate the interference of random fluctuations in measured data; thus, the influence of random errors was reduced, and real displacement results of the structure were obtained.

First, based on the moving average method of least squares, a polynomial function of degree m was obtained [15]:

$$Y(t)=a_0+a_{1}t+a_2{t}^2+\cdots +a_m{t}^m$$

(4)

The original data with the minimum mean square error were obtained by the data smoothing approach:

$$A(a_0,a_1,\cdots ,a_m)={\displaystyle \sum _{i=-u}^u{[{\displaystyle \sum _{j=0}^m{a}_j{t}_i^j-Y_i}]}^2}$$

(5)

where u is the sliding order and m is the smoothing number. In order to determine undetermined coefficients and minimize the error sum of squares between the function Y(t) and discrete data, the partial derivative of $a_k(k=0,1,\dots ,m)$ in Equation (5) was considered to be zero. Hence, the linear equations of element m + 1 were obtained as:

$${\displaystyle \sum _{i=-u}^u{Y}_i{t}_i^k={\displaystyle \sum _{j=0}^m{a}_j{\displaystyle \sum _{i=-u}^u{t}_i^{k+j}}}}$$

(6)

In Figure 1a, a standard sinusoidal superposition signal with 0.8 octaves of random noises and the drift term is displayed, and Figure 1b presents data processed by the moving average method. It is noticeable that random noises and trend terms could be simultaneously removed. Therefore, the moving average method was applied to remove the interference of noise in the following analysis.

2.3. Digital Filtering and Frequency Domain Integration

The moving average method can effectively eliminate random fluctuations induced by high-frequency interferences through multiple corrections. In this work, a high-pass filter with low-frequency attenuation was applied to remove trend terms induced by low-frequency noises. The actual sampling signal can be expressed by the algorithm of Fourier transformation in the frequency domain:

$$X(k)=\frac{1}{N}{\displaystyle \sum _{r=0}^{N-1}x(r)e^{-j2\pi k\frac{r}{N}}}$$

(7)

where N is the number of sampled data, $x(r)$ is the sampling sequence signal of N points, and $x(k)$ is the discrete spectrum.

The frequency spectrum of the fast Fourier transform (FFT) of the acceleration sequence can be expressed as [21]:

$$a(k)=\frac{2}{N}{\displaystyle \sum _{r=0}^{N-1}a(r)e^{-j2\pi k\frac{r}{N}}=a_k+b_{k}j}$$

(8)

where $a(r)$ is the acceleration sequence with N number of sampled data and $a(k)$ is the complex sequence of $a(r)$ in the frequency domain after Fourier transformation. Further, $a(k)$ can be represented as:

$$a(k)=A_k\cos (2\pi f_{k}t+{\phi }_k)$$

(9)

where $A_k=\sqrt{a_k^2+b_k^2},{\phi }_k=\arctan \frac{b_k}{a_k}$, and $f_k$ is the frequency corresponding to $a(k)$.

The structural displacement in the time domain can be cumulatively expressed as:

$$s(t)={\displaystyle \sum _{k=1}^{N}s(k)=}{\displaystyle \sum _{k=1}^N\frac{A_k}{{\omega }_k^2}\cos (2\pi f_{k}t+{\phi }_k-\pi )}={\displaystyle \sum _{k=1}^N\chi \cos (2\pi f_{k}t+{\phi }_k-\pi )}$$

(10)

In order to eliminate low-frequency trend interferences, the low-frequency attenuation method, which could retain components near the target frequency to the maximum extent, was applied. Therefore, the displacement in the frequency domain can be calculated as [21]:

$$s(t)=-F^{-1}\left[\frac{{\omega }^2}{{\omega }^4+\delta }F(a(t))\right]$$

(11)

where δ is the regularization coefficient to control the trend term error of the measured acceleration:

$$\delta =\frac{1-\alpha }{\alpha }(2\pi f_t{)}^4$$

(12)

where $\alpha$ is the target accuracy coefficient and $f_t$ is the target frequency.

The low-frequency attenuation algorithm was applied by setting a high-pass filter with a target frequency. The amplitude–frequency response characteristics of a filter with a target frequency of 2 Hz are presented in Figure 2. It is noticeable that the amplitude of the transfer function approaches 1 in the high-frequency band and is equal to 0 in the low-frequency band. Therefore, this high-pass filter could effectively retain information near the target frequency and attenuate low-frequency information simultaneously. In this work, the analysis target precision coefficient was set to 0.95 to ensure better integration precision in the small target frequency. Meanwhile, the choice of target frequency is difficult for the low-frequency attenuation algorithm. The target frequency is affected by the level of low-frequency noises, and an appropriate value should be selected to eliminate the trend term error and drift error. In this paper, the value of the target frequency is selected to be 1 Hz, where good integration accuracies can be obtained in engineering applications.

3. Modal Parameter Identification

Modal parameter identification is an essential step in structural dynamic characteristic identification. In this work, reconstructed displacements are used as inputs for parameter identification. The SSI approach was applied to identify modal parameters by state space matrices based on the singular vector decomposition (SVD) of the Hankel matrix.

3.1. State Space Model of Structural Vibration

The first-order state space expression for motion of a discrete system with stochastic noise and errors can be represented as:

$$\{\begin{array}{l}{x}_{k+1}=A{x}_k+w_k\\ y_k=C{x}_k+v_k\end{array}$$

(13)

where $x_k=x(kt)\in {ℝ}^{2n}$ is the discrete state, $A=\exp (A_c\Delta t)\in {ℝ}^{2n\times 2n}$ is the state matrix under discrete time conditions, $C\in {ℝ}^{p\times 2n}$ is the output matrix, $w_k$ is the noise induced by environmental interference and modeling inaccuracy, and $v_k$ is the measured error.

Moreover,

$$x(t)=\left(\begin{array}{c}q(t)\\ \dot{q}(t)\end{array}\right),A_c=\left(\begin{array}{cc}0& I_n\\ -M^{-1}K& -M^{-1}C\end{array}\right),B_c=\left(\begin{array}{c}0\\ M^{-1}{B}_2\end{array}\right)$$

(14)

where $A_c\in {ℝ}^{2n\times 2n}$ is the state matrix, $B_c\in {ℝ}^{2n}$ is the input matrix, $x(t)\in {ℝ}^{2n}$ is the state vector, $M,C,K\in {ℝ}^{n\times n}$ are the mass, damping, and stiffness matrices respectively, $q(t)$ is the displacement vector, and $B_2$ is the description input location matrix.

The two different types of noise signals in Equation (13) are zero mean white noises and have a covariance matrix of:

$$E\left[\left(\begin{array}{c}{w}_r\\ v_r\end{array}\right)\left(\begin{array}{cc}{w}_s^T& v_s^T\end{array}\right)\right]=\left[\begin{array}{cc}Q& S\\ S^T& R\end{array}\right]{\delta }_{rs}$$

(15)

where E is the mathematical expectation, Q, R, S are the covariance matrices for noise sequence $w_k,v_k$, and ${\delta }_{rs}$ is the Kronecker delta.

3.2. Identification of the Hankel Matrix

In this section, the solutions of the state matrix A and the output matrix C are calculated. As the Hankel matrix of a vibrating structure includes modal information, the output state space model was constructed, and the Hankel matrix was obtained from displacement data [27]:

$$H_{0/2i-1}=\frac{1}{\sqrt{j}}\left(\begin{array}{cccc}{y}_0& y_1& \dots & y_{j-1}\\ y_1& y_2& \dots & y_j\\ \dots & \dots & \dots & \dots \\ y_{i-1}& y_i& \dots & y_{i+j-2}\\ y_i& y_{i+1}& \dots & y_{i+j-1}\\ y_{i+1}& y_{i+2}& \dots & y_{i+j}\\ \dots & \dots & \dots & \dots \\ y_{2i-1}& y_{2i}& \dots & y_{2i+j+1}\end{array}\right)=\frac{\chi }{\sqrt{j}}\left(\begin{array}{cccc}{a}_0& a_1& \dots & a_{j-1}\\ a_1& a_2& \dots & a_j\\ \dots & \dots & \dots & \dots \\ a_{i-1}& a_i& \dots & a_{i+j-2}\\ a_i& a_{i+1}& \dots & a_{i+j-1}\\ a_{i+1}& a_{i+2}& \dots & a_{i+j}\\ \dots & \dots & \dots & \dots \\ a_{2i-1}& a_{2i}& \dots & a_{2i+j+1}\end{array}\right)=\left(\frac{Y_{0/i-1}}{Y_{i/2i-1}}\right)=\left(\frac{Y_p}{Y_f}\right)$$

(16)

This Hankel matrix has 2i rows and j columns. In Equation (16), $y_k$ is the displacement vector, $a_k$ is the acceleration vector, $\chi$ is the acceleration integral displacement coefficient, and $Y_p\in {ℝ}^{i\times j}$, $Y_f\in {ℝ}^{i\times j}$ represent the output matrices’ ‘past’ and ‘future’, respectively.

Now, operating the two matrices in Equation (16), the block Toeplitz matrix composed of the output covariance matrix can be expressed as:

$$T_{1/i}=Y_f{Y}_p^T=\left[\begin{array}{cccc}C{A}^{i-1}& C{A}^{i-2}G& \dots & CG\\ C{A}^{i}G& C{A}^{i-1}G& \dots & CAG\\ \dots & \dots & \dots & \dots \\ C{A}^{2i-2}G& C{A}^{2i-3}G& \dots & C{A}^{i-1}G\end{array}\right]=\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ \dots \\ C{A}^{i-1}\end{array}\right]\left[\begin{array}{ccccc}{A}^{i-1}& A^{i-2}G& \dots & AG& G\end{array}\right]={\mathrm{O}}_i{\Delta }_i$$

(17)

where $O_i$ and ${\Delta }_i$ are the extended observability matrix and the controllability matrix, respectively.

Further, the singular vector decomposition of Equation (17) leads to:

$$T_{1/i}=US{V}^T=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}{V}_1^T\\ V_2^T\end{array}\right)=U_1^T{S}_1{V}_1^T$$

(18)

where U and V are orthogonal matrices and S is a diagonal matrix composed of singular values.

Hence, the observability matrix and the state sequence can be expressed as:

$$O_i=U_1{S}_1^{1/2}$$

(19)

$${\Delta }_i=S_1^{1/2}{V}_1^T$$

(20)

Therefore, the modal parameters of the system can be analyzed by calculating the state matrix A and the output matrix C based on the matrices $O_i$ and ${\Delta }_i$.

3.3. Extraction and Selection of Modes

Eigenvalue decomposition was then performed on the state matrix A of the system [28]:

$$A=\psi \mathrm{\Lambda }{\psi }^{-1}$$

(21)

where $\mathrm{\Lambda }=\mathrm{diag}({\lambda }_i)$, ${\lambda }_i$ is the eigenvalue of the discrete-time system, and Ψ is the eigenvector matrix of the system.

According to the relationship between the eigenvalues of the discrete-time system and the continuous system, mode frequencies, damping ratios, and strain modal shapes can be obtained as:

$$f=\frac{\sqrt{real({\lambda }_i^b)+imag({\lambda }_i^b)}}{2\pi }$$

(22)

$$\xi =\frac{-real({\lambda }_i^b)}{\sqrt{real({\lambda }_i^b)+imag({\lambda }_i^b)}}$$

(23)

$$\varphi =C\psi$$

(24)

where ${\lambda }_i^b=\frac{\ln {\lambda }_i}{\Delta t}$, ${\lambda }_i^b$ is the eigenvalue of the continuous system, $\Delta t$ is the time interval, and f, $\xi$, and $\varphi$ are the mode frequency, the damping ratios, and the modal shape, respectively.

3.4. Improved Stochastic Subspace Identification Algorithm

As the signal length affects the recognition accuracy and operation efficiency of SSI, it is important to control the signal length. In this work, an improved method of modal parameter identification is presented by adding a sliding window in the calculation of modal parameters. Generally, the window size determines the accuracy of recognition. If the window is too small, modals do not appear. On the contrary, if the window is large, the calculation process becomes complex and more error modals appear in the stability diagram. Therefore, both the window size and the step size are set through sliding windows to realize modal parameter identification with high efficiency (Figure 3).

Generally, the characteristics of the tested structure, the sampling frequency, and the noise level for the chosen size of the calculated window should be comprehensively considered. Here, specific analysis steps are applied to contain these factors based on the order determination of the stability diagram:

(1)

Assume that the initial window time is $t_1$, the sampling frequency is sf, and the initial window length is $L_1=sf\times t_1$.

(2)

Let the lag time, $\Delta t$, between previous and later windows and the data length correspond to the ith window, as $L_i=[t_1+(i-1)\times \Delta t]\times sf$.

(3)

Identify the corresponding stability diagrams of each window by the SSI algorithm. If the system order is N, the maximum number of stable points will be the same corresponding to each mode. The related number of stable points for the first three modes are $f_1$, $f_2$, and $f_3$, and the percentages, $G{F}_i=f_i/N$, of these modes can be obtained.

(4)

Change the value of $G{F}_i$ with the signal length. When the percentage exceeds 90%, the corresponding window time is supposed to be the appropriate one, otherwise the window time will be increased.

4. Numerical Simulation

In order to verify the accuracy of the acceleration integral displacement and displacement-based modal parameter identification, an ideal acceleration sequence composed of three harmonic signals was constructed and the difference between theoretical and integral displacement was analyzed. Further, the singular values of the Hankel matrix constructed from acceleration and displacement data were compared. Finally, a 3-DOF system was used to numerically compare the difference between acceleration-based and reconstructed displacement-based modal parameter identification.

4.1. Synthesized Ideal Discrete Acceleration

Consider a composite-simulated sinusoidal signal that includes the frequencies of 10, 15, and 20 Hz with the corresponding acceleration amplitudes of 10, 20, and 30 g, respectively. The acceleration signal can be expressed as:

$$a(t)=10\sin (20\pi t)+20\sin (30\pi t)+30\sin (40\pi t)$$

(25)

Now, taking the sampling frequency of the above signal sequence as 1000 Hz, the theoretical displacement signal can be obtained as:

$$x(t)=-\frac{10}{{(20\pi )}^2}\sin (20\pi t)-\frac{20}{{(30\pi )}^2}\sin (30\pi t)-\frac{30}{{(40\pi )}^2}\sin (40\pi t)$$

(26)

Figure 4 presents the effects of the improved algorithm on the integral displacement. Figure 4a depicts the difference between integral and standard displacements. The error was large at 0.2 s due to the uncertainty of the initial value and the truncation error, and it fitted well after 0.2 s. Displacements at 1.175–1.195 s are locally amplified in Figure 4b, and it is noticeable that these displacement signals had the same frequency and waveform, with slightly different amplitudes. The frequency spectrum shown in Figure 4c was obtained by performing FFT transformation on theoretical and reconstructed displacement signals. It is evident that theoretical and reconstructed displacement signals had the same frequency spectrum components, which included three main frequency components of 10, 15, and 20 Hz. However, the effect of interference on signal amplitude existed under 10 Hz.

Next, the integral error evaluation index was applied to calculate the error between integral and standard displacements (Figure 4a) [29]:

$$Err=\frac{|d(t_i)-y(t_i)|}{range(d)}$$

(27)

where $d(t_i)$ and $y(t_i)$ are the standard and integral displacements at time $t_i$ respectively, and range(d) is the range of the standard displacement in the whole period and represents the difference between the maximum and minimum standard displacements.

The error diagram of reconstructed and locally amplified displacement signals in the interval of 1–1.5 s is exhibited in Figure 4d. The error between reconstructed and locally amplified displacement signals was large at 0.2 s and then dropped fast to 0.02, proving the validity of the improved algorithm. The value of the error is related to the selection of target frequency and the precision coefficient, and the maximum value of error is shown in

Table 1. The maximum value of error after 0.5 s by changing the target frequency and target precision coefficient.

Maximum Error		Frequency
		0.2 Hz	0.5 Hz	1 Hz	4 Hz	8 Hz
Precision Coefficient	0.92	0.059	0.025	0.013	0.016	0.015
	0.95	0.054	0.023	0.012	0.015	0.016
	0.98	0.047	0.021	0.011	0.017	0.020

by changing these two parameters. It was found that the better precision was obtained until 1 Hz.

Generally, singular values of a matrix possess some important hidden information. Hankel matrices consisting of acceleration and displacement data were analyzed to verify the high recognition ability of the improved algorithm. First, 3000 data points were selected, and 1500 × 1500 Hankel matrices were constructed. The singular values and ranks of different Hankel matrices were obtained and normalized [19]. The 30th-order normalized singular value curves are displayed in Figure 5, where the blue and red curves present the results obtained based on displacement and acceleration data, respectively. The sudden drops in values of the curves correspond to the modal order of the system. It is noticeable that the curve obtained based on displacement data decreased faster than the one obtained based on acceleration data, indicating an easy identification of modal parameters. Moreover, as the displacement curve was closer to the horizontal line, its de-noising effect was better than that of the acceleration curve. Therefore, the Hankel matrix constructed by displacement data was more effective to identify low-frequency modes.

4.2. 3-DOF Mass-Damping-Spring System

In order to verify the accuracy of the improved SSI method, a spring oscillator system with three-degrees-of-freedom (DOF) was considered [30] (Figure 6). The parameters of the 3-DOF system were $m_1$ = 7 kg, $m_2$ = 10 kg, $m_3$ = 10 kg, $k_1$ = 1000 N/m, $k_2$ = 2000 N/m, and $k_3$ = 100 N/m. The mass matrix, M, and stiffness matrix, K, of the system were calculated, and the damping matrix was determined based on the Rayleigh damping assumption of $C=aM+bK$.

The difference between the conventional SSI method and the improved SSI method with reconstructed displacements was numerically simulated by the 3-DOF system. The coefficients of the damping ratio were a = 0.005 and b = 0.002. The natural frequencies of the first three orders of the system were obtained as ${\omega }_1$ = 0.46 Hz, ${\omega }_2$ = 1.21 Hz, and ${\omega }_3$ = 3.84 Hz, with the corresponding damping ratios of ${\xi }_1$ = 0.59%, ${\xi }_2$ = 0.33%, and ${\xi }_3$ = 0.45%, respectively. Gaussian white noises with a power of 60 dBW were then applied as the excitation, which was generated by the WGN function in Matlab. Figure 7 presents the time histories of the incentive force and its power spectral density. The acceleration and displacement data of the 3-DOF system were obtained by the Newmark-β method and used as input signals.

Figure 8 and Figure 9 present the time domain and power spectral density (PSD) of calculated accelerations and reconstructed displacements for the time interval of $\Delta t$ = 0.01 s and the point number of N = 5000. It is noticeable that reconstructed displacements had more low-order modes, whereas accelerations were more affected by high-order modes.

The modal parameters of the system were identified by calculated accelerations and reconstructed displacements. The frequencies, modals, and damping ratios of the system were identified by its stability diagrams. The x- and y-axes of the stability diagrams shown in Figure 10 represent the frequencies and mode orders of the system, respectively. When the difference between the characteristics of high and low orders was within a certain tolerance range, a stable mode appeared. The tolerances for the frequency, the damping ratio, and the modal shape were set to 1%, 5%, and 5% respectively, and the number of Hankel matrices was c × (n − c + 1), where c = 500 and n is the data length, and the order of the state space model ranged from 2 to 50. Moreover, the appropriate window size was 10 s and the corresponding sliding step was 5 s to satisfy 50% duplication of time-domain signals between adjacent windows.

Figure 10 displays the stability diagrams of modal frequencies based on calculated accelerations and reconstructed displacements, where the blue ‘○’ marks represent the poles of frequency stability, the red ‘△’ marks represent the poles of frequency stability and damping ratio stability, and the black ‘+’ marks represent the poles of frequency stability, damping ratio stability, and vibration mode stability.

It is noticeable that the first three modes were accurately identified by reconstructed displacements when the modal order was no more than five. On the contrary, the first-order model was identified by calculated accelerations when the modal order was close to 40, causing a loss of modes. Therefore, displacement-based modal parameter identification had better stability for low-frequency modals.

Figure 11 presents the variations of the frequency and damping ratio of the system identified by displacement signals in different periods (the first and second orders corresponded to 0–10 and 5–15 s, respectively), where the dashed lines represent the real frequency and damping ratio, and the circles, the boxes, and the triangles represent the frequency and damping of modal parameter identification using reconstructed displacements at different periods. It is clear that the accuracy of modal frequency identification was high, and the first three modes were accurately identified.

However, due to the high sensitivity of the damping ratio to noises, the accuracy of modal frequency identification was higher than that of modal damping ratio identification regardless of acceleration or displacement; hence, the damping ratio was larger than the estimated error due to the influence of inherent random errors. Although some fluctuations were noticed in the damping ratio of the third mode obtained by the proposed method, the results of the first two modes had good consistency, indicating the consistency and robustness of displacement-based modal parameter identification.

5. Experimental Verification

In this section, a dynamic experiment of a cantilever plate was performed to verify the accuracy of the proposed algorithm. In this experiment, the acceleration of the plate under a transverse force was measured, and the performance of modal parameter identification was judged by the developed strategy. Finally, the responses and modals obtained by reconstructed displacements were compared with tested data.

5.1. Experimental Setup

A carbon fiber-reinforced composite cantilever plate of 540 × 180 × 1.5 mm in size was used, and its physical parameters are presented in

Table 2. Properties of the carbon fiber-reinforced composite plate.

Carbon Fiber Composite
Young’s modulus (GPa)	Ex	Ey	Ez
	47.45	60.3	3.9
Shear modulus (GPa)	Gxy	Gyz	Gxz
	72.9	1.5	62.35
Poisson’s ratio	νxy	νyz	νxz
	0.3	0.4	0.3
Density (Kg/m3)	1800

.

The experimental platform and process are presented in Figure 12a,b, respectively. The input signal of the system was generated by an Agilent 33220A signal generator (1) and amplified by a YE5872 power amplifier (2) as the input signal of an electric shaker (3). The carbon fiber-reinforced composite cantilever plate (4) was fixed on the shaker by a specially designed fixture (5). The signal acquisition system consisted of a PCB M352C65 ceramic shear ICP accelerometer (6), an acceleration signal acquisition system (7), and a computer (8). The accelerometer was fixed on the cantilever plate to measure the acceleration of the plate. The modal frequency and damping ratio of the cantilever plate were identified by two sensors—one was set at the free end (acc1) and the other was placed in the middle of the cantilever plate (acc2).

5.2. Reconstructed Displacements Calculated from Accelerations

Environmental stochastic effects were also considered in the experiment, and an electric vibration table was used to generate noise excitations with an amplitude of 10 V. Figure 13 presents the time-domain signal and time-frequency distribution of acc1 calculated by short-time Fourier (STFT) transformation, with a sampling frequency of 1024 Hz. The spectrum diagram in Figure 13b was calculated by a Hamming window with a segment length of $2^{10}$ and a step size of 1. It is noticeable that excited frequencies were the same in the entire time domain with white noise excitations. Acceleration signals in the interval of 20–30 s were locally amplified to analyze the characteristics of system responses (Figure 14, where Figure 14a,b present the time-domain signal and the spectrum diagram, respectively). The proposed frequency-domain integration algorithm was then applied to obtain the dynamic displacement responses of the structure in the interval of 20–30 s. The target frequency for reconstruction was selected as 1 Hz. Figure 15a compares the time histories of measured accelerations and reconstructed displacements. It is clear from Figure 15b that higher-order modal components for reconstructed displacement were dramatically reduced as compared to those for measured accelerations. Therefore, reconstructed displacements could eliminate peak responses in the high-frequency range of measured accelerations and also reduce the non-stationarity of tested data.

Moreover, the modal parameters of the cantilever plate were obtained through a sinusoidal frequency sweep experiment. The frequency range was set as 1–100 Hz, and the excitation amplitude was 500 mV. Figure 16 displays the spectrum of measured accelerations obtained from the sinusoidal frequency sweep experiment. The first four frequencies were detected as 4.304, 20.839, 30.518, and 55.839 Hz respectively, and the corresponding damping ratios were calculated by the half-power points of resonant peaks as 2.10%, 3.68%, 2.30%, and 3.88%, respectively.

5.3. Modal Parameter Identification

Modal parameter identification with measured accelerations and reconstruction displacements was performed with the same Hankel matrix in the same period of 20–30 s, and the corresponding results are presented in Figure 17. The range of the state space model order was 1–60, the line number of the Hankel matrix was 600, and the frequency, damping, and modal shape thresholds were 1%, 5%, and 5%, respectively. Moreover, the appropriate sliding window size was 2 s, and the corresponding sliding step was 1 s. Figure 17a displays the stability diagram obtained from modal parameter identification based on measured accelerations in the period of 0–2 s. It is noticeable that the first three modal frequencies were accurately identified from measured accelerations; however, a strong dependency of modal frequencies on the model order was detected. For example, the frequency and damping ratio of the first mode could be determined when the model order reached 30. The stability diagram obtained from modal parameter identification based on reconstructed displacements in the period of 0–2 s is presented in Figure 17b. Therefore, reconstructed displacement-based modal parameter identification could avoid the loss of low-frequency modes of the plate.

In order to avoid the error caused by random factors, the modal analysis was performed 30 times with different acceleration signals and reconstruction displacements. The mean values of calculated modal frequencies and damping ratios are presented in

Table 3. Frequency and damping ratios of the cantilever plate.

		Calculated by Sinusoidal Sweeps	Calculated by Accelerations	Calculated by Displacements	By the Method of He [1]
First-order	Frequency (Hz)	4.304	4.013	4.407	4.535
	Damping ratio	2.10%	1.42%	2.31%	2.90%
Second-order	Frequency (Hz)	20.839	20.241	20.189	20.351
	Damping ratio	3.68%	2.71%	3.49%	4.20%
Third-order	Frequency (Hz)	30.518	30.256	30.188	30.258
	Damping ratio	2.30%	1.55%	1.50%	3.50%
Fourth-order	Frequency (Hz)	55.839	51.821	55.449	53.556
	Damping ratio	3.88%	4.80%	3.50%	5.10%

. It is evident that the displacement-based modal analysis had higher accuracy as compared to the conventional acceleration-based modal parameter identification method. Meanwhile, to establish the effectiveness of the proposed method, the energy-oriented categorization of modal components, SSI [1], was used to identify the modal parameters of the structure, and the modal parameters are presented in Table 3.

6. Conclusions

An improved modal parameter identification method using reconstructed displacements was proposed in this work. The main findings are presented below:

• In comparison to the acceleration-based modal parameter identification method, the proposed method based on reconstructed displacements could suppress high-frequency components and reduce the non-stationary characteristics of measured data more effectively. Therefore, this method could provide more reliable and accurate parameter identification results for small model order selection.
• In comparison to other integration methods, the proposed method had no error accumulation and long period drift of reconstructed displacements, resulting in higher computational efficiency.
• Experimental results of the cantilever plate revealed that the improved SSI method could improve the recognition efficiency and also identify more stable results for low-frequency structures.