Singular Spectrum Analysis for Modal Estimation from Stationary Response Only

Conventional experimental modal analysis uses excitation and response information to estimate the frequency response function. However, many engineering structures face excitation signals that are difficult to measure, so output-only modal estimation is an important issue. In this paper, singular spectrum analysis is employed to construct a Hankel matrix of appropriate dimensions based on the measured response data, and the observability of the system state space model is used to treat the Hankel matrix as three components containing system characteristics, excitation and noise. Singular value decomposition is used to factorize the data matrix and use the characteristics of the left and right singular matrices to reduce the dimension of the data matrix to improve calculation efficiency. Furthermore, the singular spectrum is employed to estimate the minimum order to reconstruct the Hankel matrix; then, the excitation and noise components can be removed, and the system observability matrix can be obtained. By appropriately a factorizing system observability matrix, we obtain the system matrix to estimate the modal parameters. In addition, the fictitious modes produced by increasing the order of the matrix can be eliminated through the stabilization diagram.


Introduction
In the design of civil structures or mechanical systems, the actual working conditions such as earthquakes on buildings, the rotational speed of shafts and the operation of machine tools are often considered. If the Fourier spectrum of operational conditions has rich frequency content around the structure modes, it will cause the structure to resonate and reduce the life cycle of the structure, so the industry will include modal analysis as a reference in the design. The traditional experimental modal analysis (EMA) belongs to the frequency-domain analysis method. The dynamic information of systems extracted from frequency response function can be obtained by knowing both the excitation and response of a system. Meanwhile, in structural systems with small damping and without severe modal interference, the identification results are good in general. However, in many practical engineering structures, the structure often vibrates under unknown ambient excitation, which is often difficult to measure. Therefore, how to estimate the modal parameters of such structural systems subjected to unknown excitation for further structural control, damage detection or finite element model updating is a very important issue in structural engineering.
In 1965, Ho and Kalman [1] proposed the minimum realization algorithm, which uses Markov parameters composed of impulse response functions to obtain the state space representation of the minimum order. This method can accurately describe the system; it does not, however, effectively work on assessment with noise influence. In 1974, Zeiger and McEwen [2] proposed the concept of combining singular value decomposition (SVD) with the minimum realization method; in the 1980s, Juang and Pappa [3] developed beams and used SSA combined with a damage index to locate the possible damage in a structure. Li et al. [19] conducted a spectrum analysis on the measured responses of arch dams excited by earthquakes and used ensemble empirical mode decomposition combined with a wavelet threshold and a singular spectrum analysis to filter and highlight the dynamic characteristics of arch dams. Fitzgerald et al. [20] designed a piezoelectric energyharvesting device (EHD) in the form of a cantilever beam and used the SSA method to filter relevant interference signals, thereby accurately obtaining the signal of the measured object. Trendafilova [21] performed SSA on the impulse response of a single degree of freedom in the beam structure, decomposed the impulse response into linear and nonlinear oscillatory signals by using different orders of a singular spectrum and selected a proper order of singular spectrum to reconstruct the original signal. In general, SSA is equivalent to low-pass filtering the original signal, filtering out high-frequency noise and aperiodic signals in the signal and then smoothing the original signal so as to extract the information that mainly contributes to the dynamic characteristics of the system. This method means that SSA can preserve the characteristics of important dynamic parameters of the system in the signal and remove the irregular and random components in the time series.
This research mainly uses singular spectrum analysis (SSA) with state space for outputonly system modal parameter estimation. We construct the original response data into a Hankel matrix, for which the singular spectrum to estimate the main singular values of the matrix, and use the grouping and averaging in the SSA procedure, respectively, to remove the singular values with lower contributions and reconstruct a Hankel matrix. By introducing the observability property of the state space, a Hankel matrix with appropriate dimensions is applicable to reduce the influence of excitation and noise on the response by singular spectrum analysis, thereby effectively improving the accuracy of system identification. Finally, numerical simulations confirm the validity and reliability of singular spectrum analysis for modal estimation.

Research Methods
Singular Spectrum Analysis (SSA) is a signal decomposition method. This method is used to decompose the original signal into multiple components, determine the components to find the signal trend, eliminate the influence of noise, reconstruct the signal to make it smooth, and can determine the order of the system. The four steps singular spectrum analysis are embedding, singular value decomposition (SVD), grouping and averaging.

State Space Equation
A state space equation is applicable to describe a physical system in terms of inputs, outputs and states. This paper uses a linear time-invariant system and starts by considering the measurement data contaminated with noise in practical engineering. The discrete state space equation can be expressed as: Equation (1)  in the state space is observability, which represents the state matrix of the system inferring from the output data. The output equation based on observability can be expressed as: where [α] is called the observability matrix, which represents the matrix of the initial state of the system at each moment of the system. By factorizing the observability matrix into different time series, the system matrix [A] can then be obtained.

Embedding
Embedding is to create a trajectory matrix from the original signal. The trajectory matrix is used to arrange the measurement response signals in sections. The anti-diagonal elements correspond to the same time point, and the representative trajectory matrix is the Hankel matrix.
Consider a known time signal with measurement degrees of freedom p and signal length m, and the data matrix can be represented as Y p×m = y 1 y 2 · · · y m . We used this signal to perform segmental arrangement to create a Hankel matrix H, where the number of row vectors is K and the number of column vectors is L, as follows: where [α] and X k are the observability and state matrices of the system, respectively, and [U] is the Hankel matrix of the external force. If the randomness of noise and excitation is included in the response data, it is impossible to obtain accurate dynamic characteristics by using the response data only. Therefore, the influence of noise needs to be eliminated first. The noise may include measurement noise or model errors. Since the noise has no correlation with the dynamic characteristics of the system, the predecessors often assumed that the noise is stationary white, so that the noise term [ _ v] is uncorrelated to the dynamic characteristics of a system. To avoid the column vector not being extended enough to the original data length, the dimension of the Hankel matrix would be constructed under the restrictions of L ≤ m 2 and L ≤ K. Equation (4) indicates that the constructed Hankel matrix is composed of the information of structure itself, excitation and noise.

Singular Value Decomposition
Singular value decomposition (SVD) is a common decomposition of non-square matrix, which can decompose a matrix into multiple combinations, and is often applicable to signal processing and statistics fields. Through SVD, the factorization of H can be expressed as: where U and V are orthogonal matrices, and Σ is a singular value matrix consisting of zero off-diagonal entries and obvious non-zero singular values on the main diagonal entries. By arranging the singular values from large to small, i.e., s 1 > s 2 > · · · s pL , to obtain the singular spectrum of the signal, the relatively large singular value can be determined as the system component; otherwise, the smaller singular value will be viewed as noise.
To measure a more complete system signal, a longer acquisition time sample needs to be chosen, so the corresponding operation time due to singular value decomposition is increased. In general, the length K of the reference signal will be much larger than matrix column pL. To reduce the operation time of singular value decomposition, we construct the data matrix H pL×K H T pL×K in conjunction with the characteristics of the left and right singular matrices, and the dimension of the H pL×K H T pL×K is reduced from pL × K to pL × pL as follows: Subsequently, the right orthogonal matrix V T K×K corresponding to the matrix H pL×k is evaluated as follows: Through the matrix dimension reduction, we can define the Hankel matrix parameter K as m − L + 1, so that the parameters required to be adjusted are reduced in the later analysis, and then the system order is estimated by the stabilization diagram.

Grouping
By comparing the corresponding values in the singular spectrum, the order corresponding to the line segment with a larger slope in the singular spectrum is evaluated as the selected system order. The vectors of the corresponding order are employed to reconstruct the Hankel matrix, and then the singular value decomposition is performed as follows: where [S 1 ] i×i determines the order of the matrix to be retained, and [S 2 ] (pL−i)×(pL−K) determines the order of the matrix to be removed. The first few items with larger singular values are used to reconstruct a new Hankel matrix Ĥ , so as to retain the dynamic characteristics of the system and reduce the excitation and noise components: where Ĥ is the reconstructed Hank matrix, and i is the number of pre-reserved components.

Averaging
In the process of reconstructing the new Hank matrix Ĥ , the calculation of SVD will result in that the responses of the anti-diagonals will not be exactly equal. Since this matrix does not conform to the form of the Hank matrix, it is necessary to average the values of the anti-diagonals. We first reconstructed the signal and then reconstructed it into a Hank matrix. The correction equation is as follows: We reconstructed the above corrected time series g (n) into a Hankel matrix, as shown in the following: With the minimum order of the singular spectrum estimation matrix, the reconstructed Hank matrix [Ω] has eliminated the low-contribution components in the signal and then eliminated most of the excitation and noise information in the signal. By comparing Equation (11) with Equation (4) and removing the terms of excitation and noise, the following equation is obtained: [ By using the SVD, [Ω] can be factorized into [α] and X k as follows: [Ω] where [α] = UΣ 1 and X k = Σ 0 V T . According to Equation (4), [α] contains system matrix [A], which can be expressed as: To Since [α] 1 is not a square matrix, directly using matrix inversion cannot obtain the system matrix. Therefore, the Moore-Penrose pseudoinverse method is used to perform the generalized inverse of [α] 1 , and then the inverse matrix of [α] 1 is obtained. In practice, the singular value decomposition of [α] 1 can be evaluated as follows: Substituting Equation (17) into Equation (16), [P] can be expressed as: The system matrix [P] obtained above contains the dynamic characteristics of the system, including natural frequencies, damping ratios and mode shapes. To obtain the dynamic characteristics, the eigenvalue decomposition of the system matrix is as follows: where [Z] = diag z 1 z 2 · · · z i is the i eigenvalues of the matrix [P], and [Ψ] = ψ 1 ψ 2 · · · ψ i is the corresponding eigenvectors, that is, the mode shapes of the structural system. The eigenvalues of discrete system have the following relationship with the eigenvalues of the continuous vibration system: where ∆t is the resolution of the discrete signal, and λ R k and λ I k are the real and imaginary parts of the eigenvalue λ k , respectively. Comparing with the eigenvalues of the vibration system, the natural frequency ω k and damping ratio ξ k of the system can be obtained as follows:

Stabilization Diagram
In theory, a continuum structure has an infinite number of degrees of freedom and modes. In practice, we cannot accurately estimate how many modes are required to describe the observed dynamic behavior of the structural system and therefore cannot accurately estimate the system order. In recent years, with the improvement in computer performance, scholars have proposed stabilization diagrams in the field of system identification, which can effectively confirm the order of an unknown system. The method of the stabilization diagram is to use to the general structure as a second-order vibration system; if it is an underdamped system, its eigenvalues will have conjugate characteristics, which can be applicable to identify the true modes; otherwise, they are identified as fictitious modes. Meanwhile, when solving the eigenvalues of the system matrix, whether the eigenvalues evaluated under different system orders are stable or not can be applicable to identify if they belong to true modes.
By setting the error range of natural frequency |( f n − f n−1 )/ f n | < 0.01 and damping ratio |(ξ n − ξ n−1 )/ξ n | < 0.05, where f n and ξ n are the natural frequency and damping ratio of the nth mode, respectively, we use the above error range to determine whether the n th and n + 1 th modes are a conjugate pair and then use the eigenvector corresponding to each eigenvalue to further estimate the system order. The Modal Assurance Criteria (MAC) extensively used in the experimental modal analysis is employed check for agreement between the identified exact eigenvectors, i.e., mode shapes. The eigenvectors corresponding to the conjugate-pair eigenvalues are generally consistent. We use MAC to confirm whether the pole is stable, so as to assist in confirming the true structural modes. The definition of the MAC is: where the superscript H is the conjugate transpose operator, and φ n is the n-th eigenvector. Furthermore, consider a structural system with complex modes, and the real and imaginary parts of the identified mode shape are denoted as φ R and φ I , respectively. We evaluate the variance S xx , S yy and covariance S xy of the eigenvectors φ R and φ I and then construct S xx , S yy and S xy into a variance-covariance matrix S as follows: By solving the eigenvalue problem associated with S, the eigenvalues λ 1,2 of S can be obtained as follows: The Modal Phase Collinearity (MPC) is defined by the above eigenvalues λ 1,2 and can be employed to identify whether the mode shape corresponds to a real mode, which is expressed as follows [22]: It is known from equation (26) that the MPC values range from 0 for a mode with completely uncorrelated phase angles to 1 for a monophase result. MPC quantifies that spatial consistency of the identification results and the degree of monophase behavior by comparing the real size of the eigenvalues of S [22]. For classical normal modes, the corresponding mode shapes are real or monophase vectors, and S of the real and imaginary parts of mode shapes has only one nonzero eigenvalue; on the other hand, the λ 1 and λ 2 of S will be approximately equal if the phase angles of the identified mode shape are uncorrelated [22]. In this paper, MPC analysis is employed to identify whether the mode is a real mode. MPC is set to 0.7 as a preliminary extracting threshold to be employed to filter some fictitious modes [23], and then the real modes can be confirmed by observing the existence and changes of the poles.

Numerical Simulations and Verifications
When the structure is subjected to dynamic testing under external force excitation, the modal parameters can be estimated from the measured data of the excitation and response of the structural system. However, it is often difficult to actually perform dynamic testing of large structures, and the modal information of real structures is often difficult to measure completely. Therefore, to confirm the validity of the developed algorithm, numerical simulations are usually used to verify the feasibility of the proposed method.

6-DOF Chain Model
A 6-DOF linear chain model system, as shown in Figure 1, is employed to demonstrate the effectiveness of the proposed method [24], and the corresponding mass matrix [M], stiffness matrix [K] and damping matrix [C] of the system are expressed as follows: The damping matrix     C of the system is proportional damping, i.e., the linear combination of mass matrix     M and stiffness matrix     K . The stationary white noise is generated using the spectrum approximation method as a zero-mean bandpass noise, as shown in Figure 2, for which the power spectral density constant is ⋅ 2 0.02 N sec/ rad , and the phase angle is a random variable uniformly distributed in 0~2π . To obtain the displacement response in time domain, the sampling interval is chosen as 0.01 s, and the sampling period is 131,072 s. To simulate the situation that the system is subjected to base excitation, the system is assumed to be initially at rest, the stationary white signal excitation is applied to the sixth degree of freedom of the system, and the corresponding displacement responses of the system can be obtained using Newmark's method.

Comparison of computation efficiency of singular spectrum analysis
To effectively obtain modal estimation results from the measured complete system signal, especially for damping identification, a longer acquisition time sample is often selected, so the corresponding singular value decomposition operation time is increased. . The stationary white noise is generated using the spectrum approximation method as a zero-mean bandpass noise, as shown in Figure 2, for which the power spectral density constant is 0.02 N 2 · s/rad, and the phase angle is a random variable uniformly distributed in 0∼2π. To obtain the displacement response in time domain, the sampling interval is chosen as 0.01 s, and the sampling period is 131,072 s. To simulate the situation that the system is subjected to base excitation, the system is assumed to be initially at rest, the stationary white signal excitation is applied to the sixth degree of freedom of the system, and the corresponding displacement responses of the system can be obtained using Newmark's method. The damping matrix     C of the system is proportional damping, i.e., the linear combination of mass matrix     M and stiffness matrix     K . The stationary white noise is generated using the spectrum approximation method as a zero-mean bandpass noise, as shown in Figure 2, for which the power spectral density constant is ⋅ 2 0.02 N sec/ rad , and the phase angle is a random variable uniformly distributed in 0~2π . To obtain the displacement response in time domain, the sampling interval is chosen as 0.01 s, and the sampling period is 131,072 s. To simulate the situation that the system is subjected to base excitation, the system is assumed to be initially at rest, the stationary white signal excitation is applied to the sixth degree of freedom of the system, and the corresponding displacement responses of the system can be obtained using Newmark's method.

Comparison of computation efficiency of singular spectrum analysis
To effectively obtain modal estimation results from the measured complete system signal, especially for damping identification, a longer acquisition time sample is often selected, so the corresponding singular value decomposition operation time is increased.

Comparison of computation efficiency of singular spectrum analysis
To effectively obtain modal estimation results from the measured complete system signal, especially for damping identification, a longer acquisition time sample is often selected, so the corresponding singular value decomposition operation time is increased. The characteristics of the left and right singular matrices can be used to reduce the dimension of singular value solution, so as to improve the computational efficiency. Thus, the selection length of the signal to be analyzed is less affected by the computer performance. The simulation of the 6-DOF system is performed to compare the difference in computation time to compare the difference between the original and improved algorithm procedure of SSA. The computer operating environment is as follows: the central processing unit (CPU) is an Intel ® Core™ i7-8086K CPU, Random access memory (RAM) is 64.0 GB and the computing software is MATLAB R2020b. Moreover, the dimension L of Hankel matrix is from 3 to 60, and K is 50,000. Two stabilization diagrams of the operation results are shown in Figure 3, from which six modes can be clearly observed, and there is only subtle differences in the identification results as shown in Table 1. The computation time of the original and improved algorithm procedure of SSA are 7108.07 s and 112.20 s, respectively. The computation time before and after the improvement in SSA is increased by about 62 times, so the overall computation efficiency is significantly improved. The characteristics of the left and right singular matrices can be used to reduce the dimension of singular value solution, so as to improve the computational efficiency. Thus, the selection length of the signal to be analyzed is less affected by the computer performance. The simulation of the 6-DOF system is performed to compare the difference in computation time to compare the difference between the original and improved algorithm procedure of SSA. The computer operating environment is as follows: the central processing unit (CPU) is an Intel ® Core™ i7-8086K CPU, Random access memory (RAM) is 64.0 GB and the computing software is MATLAB R2020b. Moreover, the dimension L of Hankel matrix is from 3 to 60, and K is 50,000. Two stabilization diagrams of the operation results are shown in Figure 3, from which six modes can be clearly observed, and there is only subtle differences in the identification results as shown in Table 1

Convergence analysis with various K (reference signal length)
Convergence analysis is carried out for the parameters K of the Hankel matrix. Here, the parameter L is fixed as 60, and each increment in K is 500. The stabilization diagram of the analysis is shown in Figure 4, and the corresponding convergence trend curves with various K for damping identification is shown in Figure 5. According to the stabilization diagram, it can be found that the influence of the signal length on the frequency is not significantly different from that of the damping. From the convergence trend curves with various K for damping identification, it can be found that the damping will be greatly affected by the reference signal length. In general, the longer the signal length is, the better the identification results are.

Convergence analysis with various K (reference signal length)
Convergence analysis is carried out for the parameters K of the Hankel matrix. Here, the parameter L is fixed as 60, and each increment in K is 500. The stabilization diagram of the analysis is shown in Figure 4, and the corresponding convergence trend curves with various K for damping identification is shown in Figure 5. According to the stabilization diagram, it can be found that the influence of the signal length on the frequency is not significantly different from that of the damping. From the convergence trend curves with various K for damping identification, it can be found that the damping will be greatly affected by the reference signal length. In general, the longer the signal length is, the better the identification results are.

Convergence analysis with various K (reference signal length)
Convergence analysis is carried out for the parameters K of the Hankel matrix. Here, the parameter L is fixed as 60, and each increment in K is 500. The stabilization diagram of the analysis is shown in Figure 4, and the corresponding convergence trend curves with various K for damping identification is shown in Figure 5. According to the stabilization diagram, it can be found that the influence of the signal length on the frequency is not significantly different from that of the damping. From the convergence trend curves with various K for damping identification, it can be found that the damping will be greatly affected by the reference signal length. In general, the longer the signal length is, the better the identification results are.

Convergence analysis with various L
Convergence analysis is performed for the order L of the matrix, and its parameter K is set as − +1 m L , because it can be known from the signal length convergence analysis that the larger the K is, the more stable the damping identification result is. The stabilization diagrams with various L of the 6-DOF chain model are shown in Figure 6, the corresponding damping convergence diagrams are shown in Figure 7, and the mode shape identification results are shown in Figure 8. It is observed from the stabilization diagram that the identification result of natural frequency is relatively stable, but in the damping convergence diagram, it can be observed that the damping identification error is relatively large at low orders, and the identification results have a relatively converging trend when the order is increased. Subsequent analyses use the selected parameters  Table 2. The maximum errors in the natural frequencies and damping ratios are −0.98% and −10.45%, respectively, and the MAC values corresponding to all mode shapes are close to 1.00.

Convergence analysis with various L
Convergence analysis is performed for the order L of the matrix, and its parameter K is set as m − L + 1, because it can be known from the signal length convergence analysis that the larger the K is, the more stable the damping identification result is. The stabilization diagrams with various L of the 6-DOF chain model are shown in Figure 6, the corresponding damping convergence diagrams are shown in Figure 7, and the mode shape identification results are shown in Figure 8. It is observed from the stabilization diagram that the identification result of natural frequency is relatively stable, but in the damping convergence diagram, it can be observed that the damping identification error is relatively large at low orders, and the identification results have a relatively converging trend when the order is increased. Subsequent analyses use the selected parameters L = 60 and K = 131, 013. The identification results of the 6-DOF chain model are listed in Table 2. The maximum errors in the natural frequencies and damping ratios are −0.98% and −10.45%, respectively, and the MAC values corresponding to all mode shapes are close to 1.00.

2-D Truss Structure
To examine the effectiveness of the present method for more complex structural systems, a two-dimensional truss model with proportional damping [25] is considered, as shown in Figure 9. This system has a total number of eight nodes, of which four are fully restrained, and hence the total number of active DOFs is eight (one horizontal and one vertical per each node). The mass, damping and stiffness matrices for this system are given as follows: [M] = diag 100 100 100 100 100 100 100 100 Kg (30)

2-D Truss Structure
To examine the effectiveness of the present method for more complex structural systems, a two-dimensional truss model with proportional damping [25] is considered, as shown in Figure 9. This system has a total number of eight nodes, of which four are fully restrained, and hence the total number of active DOFs is eight (one horizontal and one vertical per each node). The mass, damping and stiffness matrices for this system are given as follows:  In this example, the previous stationary white-noise (with the same power spectral density constant still used, as shown in Figure 10, and a random distribution of phase angle) was used as input acting horizontally and vertically at all active DOFs of the truss model, and then the corresponding displacement responses obtained by Newmark's method were used for modal estimation. The initial condition of this truss structural In this example, the previous stationary white-noise (with the same power spectral density constant still used, as shown in Figure 10, and a random distribution of phase angle) was used as input acting horizontally and vertically at all active DOFs of the truss model, and then the corresponding displacement responses obtained by Newmark's method were used for modal estimation. The initial condition of this truss structural system is assumed to be at rest. The choice of the sampling interval and the sampling period are the same as the previous example. SSA is employed to analyze the simulated response data, and the corresponding stabilization diagram is shown in Figure 11. The results of modal identification with the selection of L = 60 and K = 131, 013 are summarized in Table 3, which shows that the modal identification in this case is satisfactory. The maximum errors in natural frequencies and damping ratios are −0.43% and −16.20%, respectively. The identification results of mode shapes are shown in Figure 12. Observing the MAC values, which signify the consistency between the identified and the theoretical mode shapes, all modes are found to be identified accurately. system is assumed to be at rest. The choice of the sampling interval and the sampling period are the same as the previous example. SSA is employed to analyze the simulated response data, and the corresponding stabilization diagram is shown in Figure 11. The results of modal identification with the selection of = 60 L and = 131013 K are summarized in Table 3, which shows that the modal identification in this case is satisfactory. The maximum errors in natural frequencies and damping ratios are −0.43% and −16.20%, respectively. The identification results of mode shapes are shown in Figure 12. Observing the MAC values, which signify the consistency between the identified and the theoretical mode shapes, all modes are found to be identified accurately.    Figure 11. A stabilization diagram of 2-D truss structure subjected to stationary white noise excitation. Figure 11. A stabilization diagram of 2-D truss structure subjected to stationary white noise excitation.

7-DOF Model of a Sedan
To demonstrate this method effectiveness on relatively complex 3D structure systems, a sedan model with two groups of similar modes is considered [24]. The system is a 7-DOF system, as shown in Figure 13, including the bounce, pitch and roll motion of the body, as well as the bouncing motion of the four suspension systems of a sedan, in which the related parameters are shown in Ref. [24]. The mass matrix [M], stiffness matrix [K] and damping matrix [C] of the system are expressed as follows: [M] = diag m l y l x m 1 m 2 m 3 m 4 (33)

7-DOF Model of a Sedan
To demonstrate this method effectiveness on relatively complex 3D structure systems, a sedan model with two groups of similar modes is considered [24]. The system is a 7-DOF system, as shown in Figure 13, including the bounce, pitch and roll motion of the body, as well as the bouncing motion of the four suspension systems of a sedan, in which the related parameters are shown in Ref. [24]. The mass matrix     M , stiffness matrix     K and damping matrix     C of the system are expressed as follows:   Considering that the vehicle will experience bumps on the road when driving, the simulated stationary white noise serves as the excitation input acts on to the DOFs of the four suspension components and bouncing of car body of a sedan. The condition of this sedan structural system is assumed to be initially at rest. The choice of the sampling in- Considering that the vehicle will experience bumps on the road when driving, the simulated stationary white noise serves as the excitation input acts on to the DOFs of the four suspension components and bouncing of car body of a sedan. The condition of this sedan structural system is assumed to be initially at rest. The choice of the sampling interval and the sampling period are the same as the previous example. Figure 14 is the stabilization diagram of the 7-DOF sedan model, roughly showing 5 frequencies of about 5, 7, 19, 70 and 81 rad/s. Two groups of closely spaced higher modes cannot be effectively estimated from the power spectrum only. By using the distribution of poles in the stabilization diagram, it can be observed that there are two obvious poles in 70 rad/s. The order of the Hankel matrix is needed to be increased to possibly observe two severely closely spaced modes around 81 rad/s. The results of modal identification with the selection of L = 60 and K = 131, 013 for the 7-DOF sedan model by SSA are shown in Table 4 and Figure 15. The maximum errors in natural frequencies and damping ratios are −6.89% and −26.38%, respectively, and the lowest MAC in mode shapes is 0.66. terval and the sampling period are the same as the previous example. Figure 14 is the stabilization diagram of the 7-DOF sedan model, roughly showing 5 frequencies of about 5, 7, 19, 70 and 81 rad/s. Two groups of closely spaced higher modes cannot be effectively estimated from the power spectrum only. By using the distribution of poles in the stabilization diagram, it can be observed that there are two obvious poles in 70 rad/s. The order of the Hankel matrix is needed to be increased to possibly observe two severely closely spaced modes around 81 rad/s. The results of modal identification with the selection of = 60 L and = 131013 K for the 7-DOF sedan model by SSA are shown in Table 4 and Figure 15. The maximum errors in natural frequencies and damping ratios are −6.89% and −26.38%, respectively, and the lowest MAC in mode shapes is 0.66. Figure 14. A stabilization diagram of 7-DOF model of a sedan subjected to stationary white noise excitation. Figure 14. A stabilization diagram of 7-DOF model of a sedan subjected to stationary white noise excitation. terval and the sampling period are the same as the previous example. Figure 14 is the stabilization diagram of the 7-DOF sedan model, roughly showing 5 frequencies of about 5, 7, 19, 70 and 81 rad/s. Two groups of closely spaced higher modes cannot be effectively estimated from the power spectrum only. By using the distribution of poles in the stabilization diagram, it can be observed that there are two obvious poles in 70 rad/s. The order of the Hankel matrix is needed to be increased to possibly observe two severely closely spaced modes around 81 rad/s. The results of modal identification with the selection of = 60 L and = 131013 K for the 7-DOF sedan model by SSA are shown in Table 4 and Figure 15. The maximum errors in natural frequencies and damping ratios are −6.89% and −26.38%, respectively, and the lowest MAC in mode shapes is 0.66.

7-DOF Highly Coupled System
In this case, a highly coupled system [26] is used to verify the effectiveness of the proposed method. A schematic representation of this model is shown in Figure 16. The mass matrix     M , stiffness matrix     K and damping matrix     C of the system are expressed as follows: Figure 15. Comparison between the identified mode shapes and the exact mode shapes of the 7-DOF model of a sedan subjected to stationary white-noise input.
In this example, we still use the previous stationary white-noise (with the same power spectral density constant as shown in Figure 17 and random distribution of phase angle) as input acting all DOFs of this model, and then the corresponding displacement responses obtained by Newmark's method are used for modal estimation. We assume that the initial condition of this truss structural system is at rest. The choice of the sampling interval is chosen as 0.001 s, and the sampling period is 131.072 s. SSA is employed to analyze the simulated response data, and the corresponding stabilization diagram is shown in Figure 18. Table 5 and Figure 19 summarize the results of the modal identification, with the selection of = 90 L and = 130983 K in this case, and present the well-identified modal parameters of the system. The maximum errors in natural frequencies and damping ratios are −0.98% and 28.90%, respectively, and the lowest MAC in mode shapes is 0.98. In this example, we still use the previous stationary white-noise (with the same power spectral density constant as shown in Figure 17 and random distribution of phase angle) as input acting all DOFs of this model, and then the corresponding displacement responses obtained by Newmark's method are used for modal estimation. We assume that the initial condition of this truss structural system is at rest. The choice of the sampling interval is chosen as 0.001 s, and the sampling period is 131.072 s. SSA is employed to analyze the simulated response data, and the corresponding stabilization diagram is shown in Figure 18. Table 5 and Figure 19 summarize the results of the modal identification, with the selection of L = 90 and K = 130983 in this case, and present the well-identified modal parameters of the system. The maximum errors in natural frequencies and damping ratios are −0.98% and 28.90%, respectively, and the lowest MAC in mode shapes is 0.98.      Figure 19. A model of 7-DOF highly coupled system.

Conclusions
The singular spectrum analysis (SSA) in the present paper can be effectively applicable to modal estimation from stationary response only. In this paper, the system response data is composed of Hankel matrix, and the appropriate factorization of data matrix under consideration is to process through embedding, singular value decomposition (SVD), grouping, and averaging in SSA. The dynamic information of systems is effectively extracted, the influence of external force and noise on the data matrix is reduced and the corresponding observability matrix of the system to be identified can be obtained. Then, the modal parameters are estimated through the appropriate factorization of the observability matrix. In addition, the characteristics of the left and right singular matrices in SVD can be employed to effectively reduce the dimension of the data matrix. Therefore, the computation efficiency of the response data matrix constructed from long-time data samples is significantly improved and without the restriction of computer performance to affect the length of the chosen response data. In additional, when SSA is applicable to identify modal parameters, increasing the system order can make SVD more effective in removing non-structural information, but it may yield some fictitious modes due to numerical computation; the stabilization diagram is employed to determine whether they are fictitious modes and to determine the system modes.