*Article* **Response-Only Parametric Estimation of Structural Systems Using a Modified Stochastic Subspace Identification Technique**

**Chang-Sheng Lin \* and Yi-Xiu Wu**

Department of Vehicle Engineering, National Pingtung University of Science and Technology, Pingtung 912301, Taiwan; m10938002@mail.npust.edu.tw

**\*** Correspondence: changsheng@mail.npust.edu.tw

**Abstract:** The present paper is a study of output-only modal estimation based on the stochastic subspace identification technique (SSI) to avoid the restrictions of well-controlled laboratory conditions when performing experimental modal analysis and aims to develop the appropriate algorithms for ambient modal estimation. The conventional SSI technique, including two types of covariance-driven and data-driven algorithms, is employed for parametric identification of a system subjected to stationary white excitation. By introducing the procedure of solving the system matrix in SSI-COV in conjunction with SSI-DATA, the SSI technique can be efficiently performed without using the original large-dimension data matrix, through the singular value decomposition of the improved projection matrix. In addition, the computational efficiency of the SSI technique is also improved by extracting two predictive-state matrixes with recursive relationship from the same original predictive-state matrix, and then omitting the step of reevaluating the predictive-state matrix at the next-time moment. Numerical simulations and experimental verification illustrate and confirm that the present method can accurately implement modal estimation from stationary response data only.

**Keywords:** operational modal analysis; ambient modal analysis; stochastic subspace identification; singular value decomposition; stationary white noise

#### **1. Introduction**

The dynamic characteristics of a structural system, such as natural frequencies, damping ratios, and mode shapes, can be investigated through numerical and experimental analysis. The response of a structural system is measured with a known excitation in modal testing, which is usually performed under well-controlled laboratory conditions. However, performing experimental modal analysis in real operating conditions may be possible, even for large and complex mechanical systems with real boundary conditions [1,2]. The modal parameters obtained theoretically under the free boundary condition can be calculated by mathematical modeling to obtain the characteristics under arbitrary boundary constraints [3]. However, experimental results obtained under specified boundary conditions cannot be converted to other dynamic characteristics under the constraints of other boundaries. Therefore, it is difficult to perform modal testing under practical boundary conditions. The hammer excitation testing method is generally used to measure the frequency response function of the structural system, and then parametric estimation is performed to understand the dynamic characteristics of the structural system [4].

The system-identification methods described above are generally used to systematically determine or improve a mathematical model for a physical system and are implemented by measuring both observed structural excitation and corresponding response data. However, an obvious difference exists between the operating conditions of realistic structures in practical work and a controlled-environment laboratory in modal testing [5]. Dynamic characteristics cannot fully represent the system mode under real-world operating conditions; thus, it is necessary to study how to perform modal identification of systems in authentic operating environments [6].

**Citation:** Lin, C.-S.; Wu, Y.-X. Response-Only Parametric Estimation of Structural Systems Using a Modified Stochastic Subspace Identification Technique. *Appl. Sci.* **2021**, *11*, 11751. https://doi.org/ 10.3390/app112411751

Academic Editors: Nikos D. Lagaros and Cecilia Surace

Received: 25 September 2021 Accepted: 6 December 2021 Published: 10 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Operational modal analysis [7], which is also called "ambient modal analysis", or "output-only modal analysis" [8], is extensively used in modal estimation of large structures under environmental and operational loads [9], such as vehicle suspension systems [10], offshore wind power facilities [11], and stadium structures [12]. Many identification methods have been extensively employed for modal extraction based on ambient response. In 1993, the so-called natural excitation technique (NExT) was proposed and initially used for modal estimation of structures in wind engineering, by assuming that ambient excitation is stationary white noise [13]. It was employed to replace free or impulse response in conventional modal estimation methods in the time domain. Subsequently, if ambient excitation can be expressed as a product model of stationary white noise and an envelope function describing the same variation of time history as excitation amplitude, the corresponding response of a structural system can be converted approximately into free response through the correlation technique [14] or random decrement technique [15]. Modal estimation can then be carried out, using the parametric estimation technique in the time domain. In addition, by introducing the correlation matrix between ambient response data to the procedure of ERA/DC, the ERA/DC can be effectively applied to modal identification of structures subjected to stationary white excitation, even to the practical recorded excitation of an earthquake [16].

In recent years, Stochastic Subspace Identification (SSI), applied with NExT, has been widely employed to modal estimation of structures under ambient vibration [17]. The SSI method is a time-domain modal-estimation method under the assumption of stationary white noise for ambient excitation and can be directly applied to modal estimation from ambient response records only [18]. There is no need for excitation measurement; thus, it is suitable for the analysis of ambient vibration. In addition, among the algorithms for structural health monitoring (SHM) to perform modal identification of structural systems, SSI is a reliable time-domain technique using extended observability matrices [19]. Numerous studies have specifically concentrated on realistic applications of SSI in recent years. The SSI-COV method uses the calculation of correlation function through the output data and then constructs a correlation matrix. The observability matrix can be obtained by using the singular value decomposition (SVD) of the correlation function matrix, and then the modal parameters can be estimated. In 1993, SSI-DATA was proposed, based on the concept of Kalman filter and space-vector projection [18]. Through the projected output matrix obtained by projecting the output vector of the future into the output vector space of the past, we substitute the projected output matrix into the original correlation function matrix. The modal parameters can be estimated from the observability matrix obtained by SVD of the projected matrix [20].

SSI-DATA is relatively complete in the derivation process under the signal length limitation of general response data, but there are some cases where the calculation efficiency is poor. In this study, we introduce correlation function calculations in the SSI-COV system matrix method into the SSI-DATA algorithm. Through the SVD of the improved projection matrix, low computational efficiency due to the large matrix dimension can be avoided. By extracting two predictive-state matrixes with recursive relationships from the same original predictive-state matrix, the efficiency of computation can be improved, and the step of reevaluating the predictive-state matrix at the next-time moment can then be omitted.

#### **2. Stochastic Subspace Identification Method**

The analysis of the stochastic subspace identification (SSI) method is based on the framework of the state-space model. To treat the measurement data with the SSI method, this method is derived from the continuous-time domain to discrete time. Because the SSI method can be used to process the output-only system, which is different from the deterministic state space, we consider the input to be a stationary random process that can be expressed as a random, discrete-time, state-space equation.

Since the identification process of the SSI method can be implemented from the output measurement data only, the ambient excitation is assumed to be white noise input, without considering external force input. Therefore, the external force and the noise can be combined as white noise. To apply the measurement data to the SSI method, we can construct a Hankel matrix [*H*] composed of the measurement data, from which the relationship between the different measurement channels and different sampling times are as follows:

$$[H]\_{2li \times j} = \begin{bmatrix} \underline{y}\_0 & \underline{y}\_1 & \cdots & \underline{y}\_{j-1} \\ \underline{y}\_1 & \underline{y}\_2 & \cdots & \underline{y}\_j \\ & \vdots & \vdots & \ddots & \vdots \\ & & \underline{y}\_{i-1} & \underline{y}\_i & \cdots & \underline{y}\_{i+j-2} \\ & & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ & & \underline{y}\_i & \underline{y}\_{i+1} & \cdots & \underline{y}\_{i+j-1} \\ & & \underline{y}\_{i+1} & \underline{y}\_{i+2} & \cdots & \underline{y}\_{i+j+1} \\ & & \vdots & \vdots & \ddots & \vdots \\ & & \underline{y}\_{2i-1} & \underline{y}\_{2i} & \cdots & \underline{y}\_{2i+j-2} \end{bmatrix}\_{2li \times j} = \begin{bmatrix} [Y\_j]\_{(li) \times j} \\ \cdots \\ [Y\_j]\_{(li) \times j} \end{bmatrix} = [a]\_{2li \times 2n} [\dot{X}]\_{2n \times j} \tag{1}$$

where the upper half of this matrix is called "the past" and denoted [*Yp*], and the lower half of the matrix is called "the future" and is denoted [*Yf* ] [16].

A conditional mean for Gaussian processes can be completely described by its covariance. Since the shifted data matrices are also defined as covariance, the projection can be calculated directly. Note that the state matrix estimated by Kalman filter, and the state-space model can be constructed by the measured output vector used to estimate the predictive-state matrix [*X*ˆ*i*]. The projection matrix [Ω] can be expressed as a product of the observability matrix [*α*] and the predictive-state matrix [*X*ˆ*i*] of the Kalman filter in the following [18]:

$$\begin{aligned} \left[\Omega\right] \quad &= E\left( \left[\boldsymbol{Y}\_f\right] \left[\boldsymbol{Y}\_p\right]^T \right) \left( E\left( \left[\boldsymbol{Y}\_p\right] \left[\boldsymbol{Y}\_p\right]^T \right) \right)^{\stackrel{\stackrel{\text{c}}{\rightleftharpoons}}{}} \left[\boldsymbol{Y}\_p\right] \\ &= \left[\boldsymbol{a}\right]\_{li \times 2n} \left[\left.\dot{\mathbb{X}}\_i\right]\_{2n \times j} \end{aligned} \tag{2}$$

where

$$[\mathfrak{a}]\_{li \times 2n} = \begin{bmatrix} [\mathbb{C}] \\ [\mathbb{C}][A] \\ \dots \\ [\mathbb{C}][A]^{i-1} \end{bmatrix}\_{li \times 2n} \tag{3}$$

$$[\hat{X}\_i]\_{2n \times j} = \left[ \begin{array}{cc} \underbrace{\mathfrak{x}}\_{i} & \underbrace{\mathfrak{x}}\_{i+1} & \cdots & \underbrace{\mathfrak{x}}\_{i+j-1} \end{array} \right]\_{2n \times j} \tag{4}$$

[*C*] is the output/observation matrix; ⊕ is Moore–Penrose pseudoinverse; [*A*] is the system matrix; *E*[·] is the expectation operator. In the first line of Equation (2), the first four matrices in the product introduce the covariance between channels at different time delays, and the last matrix in this product defines the conditions. By using the SVD analysis and choosing the effective singular-value number, [Ω] can be expressed in minimum order realization as:

$$\begin{aligned} \left[\Omega\right] &= \quad \left[\![\![\Lambda]\!]\!\right] \left[\![\![V]\!]\!\right]^T\\ &= \left(\begin{array}{cc} \left[\![\![I\_1\}\!]\!\right]\!\right) \left(\begin{array}{cc} \left[\![\Lambda\_1]\!\right] & 0\\ 0 & \left[\![\Lambda\_2]\!\right]\!\approx 0 \end{array}\right) \left(\begin{array}{cc} \left[\![V\_1\right]\!\right]^T\\ & \left[\![V\_2]\!\right]^T \end{array}\right) \\ &\approx \left[\![I\_1\!]\!\right]\_{li\times2n} \left[\![\Lambda\_1]\!\right]\_{2n\times2n} \left[V\_1\right]^T \end{aligned} \tag{5}$$

where [*U*] and [*V*] are both unitary matrixes, and [Δ] is a matrix containing singular values. The dimension of [Δ1] can, in general, be employed to estimate the system order or number of poles. However, in practical work, the partial diagonal terms of the singular-value

matrix [Δ] may be nonzero, produced by noise from the procedure of data acquisition and numerical truncation.

Through the elimination of partial matrix [Δ2], consisting of the smaller singular values, a minimum realization is obtained that results in a minimum order system representing the structural system. In Equation (5), we can, therefore, choose the number of effective singular values to obtain the minimum order realization through the SVD analysis of [Ω]. From Equations (2) and (5), with appropriate partitioning of [*α*] and [*X*ˆ*i*], the following equations can be written:

$$\begin{aligned} [\boldsymbol{\kappa}] &= [\boldsymbol{\mathcal{U}}\_1][\boldsymbol{\Delta}\_1]^\mathfrak{s} \\ [\boldsymbol{\mathcal{X}}\_i] &= [\boldsymbol{\Delta}\_1]^\mathfrak{t} [\boldsymbol{V}\_1]^\mathfrak{T} \end{aligned} \tag{6}$$

where *s* + *t* = 1. Indeed, one possible choice is [*α*]=[*U*1][Δ1] 1/2 and [*X*ˆ*i*]=[Δ1] 1/2[*V*1] *T*, which appears to make both [*α*] and [*X*ˆ*i*] balanced.

However, poor computational efficiency may occur, caused by relatively large dimensions of [Ω]. In this paper, we construct the data matrix [Ω][Ω] *<sup>T</sup>* composed of [Ω], and perform the SVD analysis of [Ω][Ω] *<sup>T</sup>* to determine the order of a structural system to be identified. It can be shown that the eigenvalue of [Ω][Ω] *<sup>T</sup>* is the square roots of the eigenvalues of ([Ω][Ω] *<sup>T</sup>*)([Ω][Ω] *<sup>T</sup>*), and that the corresponding eigenvectors of [Ω][Ω] *T* are the same as those of ([Ω][Ω] *<sup>T</sup>*)([Ω][Ω] *<sup>T</sup>*). The dimension *li* <sup>×</sup> *<sup>j</sup>* of [Ω] can be reduced to the dimension *li* × *li* of ([Ω][Ω] *<sup>T</sup>*)([Ω][Ω] *<sup>T</sup>*), where *j* >> *li*. Based on the above, the efficiency of modal estimation can be improved, and system order can be determined through the SVD analysis of ([Ω][Ω] *<sup>T</sup>*)([Ω][Ω] *T*).

In addition, to further improve the efficiency of the SSI method, we consider extracting the predictive-state matrixes [*X*ˆ *extract*1] and [*X*ˆ *extract*2] with a recursive relationship directly from the original predictive-state matrix [*X*ˆ*i*], as described next. From Equation (2), the predictive-state matrix [*X*ˆ*i*] can be obtained from the observation matrix [*α*] in the following:

$$\begin{aligned} [\mathcal{X}\_i] &= [a]^{\oplus}[\Omega] \\ &= \begin{bmatrix} \hat{\mathfrak{X}}\_i & \hat{\mathfrak{X}}\_{i+1} & \cdots & \hat{\mathfrak{X}}\_{i+j-3} & \hat{\mathfrak{X}}\_{i+j-2} & \hat{\mathfrak{X}}\_{i+j-1} \end{bmatrix}\_{2n \times j} \end{aligned} \tag{7}$$

From the measured stationary responses at *n* stations on a structure under test, we define a system matrix [*A*], such that

$$[A]\left[\hat{X}\_{\text{extract1}}\right] = \left[\hat{X}\_{\text{extract2}}\right] \tag{8}$$

where [*X*ˆ *extract*1] is a predictive-state matrix of measured response from [*X*ˆ*i*], and [*X*ˆ *extract*2] is a predictive-state matrix of time-delayed response from [*X*ˆ*i*] as follows

$$\begin{aligned} \left[ \left[ \mathcal{X}\_{\textit{extract1}} \right] = \left[ \underline{\mathfrak{k}}\_{i} \,\, \underline{\mathfrak{k}}\_{i+1} \,\, \cdots \,\, \underline{\mathfrak{k}}\_{i+j-3} \,\, \underline{\mathfrak{k}}\_{i+j-2} \right]\_{2n \times j-1} \\ \left[ \left[ \mathcal{X}\_{\textit{extract2}} \right] = \left[ \underline{\mathfrak{k}}\_{i+1} \,\, \cdots \,\, \underline{\mathfrak{k}}\_{i+j-3} \,\, \underline{\mathfrak{k}}\_{i+j-2} \,\, \underline{\mathfrak{k}}\_{i+j-1} \right]\_{2n \times j-1} \end{aligned} \tag{9}$$

Therefore, following almost the same procedure as used in Equation (7), the system matrix [*A*] can be obtained through the least-squares method. By extracting the predictive-state matrixes [*X*ˆ *extract*1] and [*X*ˆ *extract*2] with a recursive relationship directly from the original predictive-state matrix [*X*ˆ*i*], we can then avoid the step of reevaluating the predictive-state matrix at the next-time moment in the conventional SSI method, which can further improve the computational efficiency of the SSI method.

We can further solve the eigenproblem of the system matrix [*A*] to obtain the dynamic characteristics of the system, and the characteristic equation can be written as:

$$[A][\Psi] = [\Psi][\Lambda] \tag{10}$$

where [Ψ] consists of eigenvectors, i.e., mode shapes, and [Λ] contains eigenvalues *λi*. The relationship between discrete-time matrix [*A*] and continuous-time matrix [*AS*] can be expressed as:

$$\mathbb{E}\left[A\right] = \mathfrak{e}^{\left[A\_{\mathcal{S}}\right]\Delta t} \in \mathbb{R}^{2n \times 2n} \tag{11}$$

Denote the eigenvalues of [*A*] and [*AS*] as *λ<sup>i</sup>* and *λsi*, respectively. The relationship between *λ<sup>i</sup>* and *λsi* can then be expressed as

$$
\lambda\_{\bar{i}} = \mathfrak{e}^{\lambda\_{\bar{i}i} \Delta t} \tag{12}
$$

Through the eigenvalue analysis associated with the continuous-time system matrix, [*AS*], the eigenvalues *λsi* can be obtained as:

$$
\lambda\_{\rm si} = -\omega\_i \mathfrak{z}\_i \pm j\omega\_i \sqrt{1 - \mathfrak{z}\_i^{\mathfrak{z}}} \tag{13}
$$

Set the eigenvalues *λsi* of continuous-time system matrix [*AS*]

$$
\lambda\_{si} = a\_i + jb\_i \tag{14}
$$

The natural frequencies *ω<sup>i</sup>* and damping ratios *ξ<sup>i</sup>* of the structural system can be obtained as:

$$\begin{array}{l} \omega\_{i} = \sqrt{a\_{i}^{2} + b\_{i}^{2}}\\ \xi\_{i} = -\frac{a\_{i}}{\sqrt{a\_{i}^{2} + b\_{i}^{2}}} \end{array} \tag{15}$$

Consequently, the parametric estimation of structures can be implemented through the eigenvalue analysis associated with the system matrix, [*AS*], once the system matrix [*AS*] is obtained through the least-squares estimate from measured response data.

#### **3. Numerical Simulations and Experimental Verification**

#### *3.1. Six DOF Chain Model of a Cantilever Beam*

To illustrate and confirm the validity of the proposed method in this paper, we first consider a numerical example of a chain model with six degrees of freedom (6-DOF) to simulate a cantilever beam, as shown schematically in Figure 1. The masses *m*1, *m*2, *m*3, *m*4, *m*5, and *m*<sup>6</sup> for the 6-DOF chain model are equal to 2, 2, 2, 2, 3, and 4 kg. The stiffnesses *k*1, *k*2, *k*3, *k*4, *k*5, and *k*<sup>6</sup> are equal to 1, 1, 1, 1, 2, and 3 N/m, respectively. The damping matrix [*C*] of the system is in the form of [*C*] = 0.1[*M*] + 0.001[*K*] <sup>N</sup>·<sup>s</sup> <sup>m</sup> , which indicates that this structure contains proportional damping because of damping matrix [*C*], expressed as the linear combination of the mass matrix [*M*] and stiffness matrix [*K*]. The excitation force is simulated as stationary white noise, which is approximately generated as a zeromean band-pass noise [21], whose frequency range is from 0 to 50 Hz, and the standard deviation, i.e., power spectrum density, is 0.04 <sup>N</sup>2·s/rad. The sampling interval is chosen as Δ*t* = 0.01 s, and the sampling period, as shown in Figure 2, is *T* = *Nt*·Δ*t* = 1310.72 s, where *Nt* was chosen as 217. The cut-off frequency, *ω<sup>c</sup>* is 314.15 rad/s, and the resolution in frequency domain <sup>Δ</sup>*<sup>ω</sup>* is 4.79 <sup>×</sup> <sup>10</sup>−<sup>3</sup> rad/s. The sampling interval is chosen as 0.01 s, and the sampling period is 150 s. Through Newmark's method, the displacement responses of the system are obtained and then employed for modal estimation through the modified SSI method.

**Figure 1.** Schematic plot of a 6-DOF chain model of a cantilever beam.

**Figure 2.** A sample function of stationary white noise in the time domain.

By introducing the procedure of solving the system matrix in SSI-COV, in conjunction with SSI-DATA, we can increase the computational efficiency without using the original large-dimension data matrix, through the singular value decomposition of the improved projection matrix. To make the number of modes evaluated through the SSI algorithm equal to or greater than the number of modes to be identified, the dimensions of the Hankel matrix must be not less than the system order to be identified (*li* ≥ 2*n*). Through the channel expansion technique, we set the number of expansion channels to 20.

In practical engineering analysis, a continuum structure has an infinite number of degrees of freedom and modes; thus, the dimensions of the Hankel matrix depend on the number of modes to be identified. The projection matrix, [Ω] is obtained by Hankel matrix calculation. The SVD analysis of the projection matrix [Ω] can then be performed, and the number of singular values is employed to determine the order of system to be identified. In Figures 3 and 4, the distribution of the singular values of [Ω][Ω] *<sup>T</sup>* shows a relatively more obvious drop than those of [Ω] around the singular value number 12 and can be further employed to estimate the system order and the number of modes to be identified. In addition, the SVD analysis of [Ω][Ω] *<sup>T</sup>* can be employed to reduce a greater number of calculations, rather than performing the SVD of [Ω] or system matrix [*A*]. A comparison of the computational efficiency of SVD between [Ω] and [Ω][Ω] *<sup>T</sup>* for different sampling

points is shown in Table 1, in which we clearly see the computation time of SVD of [Ω] is much longer than that of [Ω][Ω] *<sup>T</sup>*, especially for [Ω] with high dimensions. In addition, through the stabilization diagram of different polynomial orders constructed by identified natural frequency, corresponding to eigenvalues from the frequency response function matrix with different polynomial orders, estimation results can then be sorted as either structural or fictitious modes. In addition, the number of structural modes can clearly be seen to be six, as shown in Figure 5.

**Table 1.** Comparison of computation time of SVD of [Ω] and [Ω][Ω] *<sup>T</sup>* for different sampling points as well as matrix multiplication processing time for [Ω][Ω] *T*.


Central processing unit (CPU) is Intel® Core™ i7-9700K CPU @ 3.60 GHz. Random-access memory (RAM) is 64.0 GB.

**Figure 3.** Distribution of the singular values of projection matrix [Ω] from the stationary responses of a 6-DOF chain model of a cantilever beam.

**Figure 4.** Distribution of the singular values of the data matrix, [Ω][Ω] *<sup>T</sup>*, constructed by projection matrix, [Ω], from the stationary responses of a 6-DOF chain model of a cantilever beam.

**Figure 5.** Typical plot of the stabilization diagram of stationary response of first DOF of a 6-DOF chain model of a cantilever beam.

In the stabilization diagram, a clear "location" of modal frequency obviously exists, even though no obvious peaks appear among close modes in the frequency response function because of modal interference. The system matrix, [*A*], can be obtained from Equations (9)–(11), and the modal parameters can then be estimated through the eigenvalue analysis of system matrix, [*A*], as summarized in Tables 2 and 3, where we clearly see the well-estimated structural modal parameters through the Modified Stochastic Subspace Identification (MSSI) method of both [Ω] and [Ω][Ω] *<sup>T</sup>*. Note that the MSSI method is a combination of the conventional SSI method and SVD analysis. A comparison between the exact and identified mode shapes is shown in Figure 6, and the corresponding Modal Assurance Criterion (MAC) [22] values evaluated are shown in Figure 7, where good agreement is observed.


**Table 2.** Identification results of 6-DOF chain model through Modified Stochastic Subspace Identification (MSSI) method from projection matrix [Ω].

**Table 3.** Identification results of 6-DOF chain model through Modified Stochastic Subspace Identification (MSSI) method from the data matrix [Ω][Ω] *<sup>T</sup>* constructed by projection matrix.


**Figure 6.** *Cont*.

**Figure 6.** Comparison between the identified mode shapes and the exact mode shapes of a 6-DOF chain model of a cantilever beam subjected to stationary white noise: (**a**) 1st bending mode shape; (**b**) 2nd bending mode shape; (**c**) 3rd bending mode shape; (**d**) 4th bending mode shape; (**e**) 5th bending mode shape; (**f**) 6th bending mode shape.

**Figure 7.** Typical plot of the Modal Assurance Criterion (MAC) for the identified mode shapes and exact mode shapes of a 6-DOF chain model of a cantilever beam.

#### *3.2. Six DOF Railway Vehicle Model with Modal Interference*

In the previous examples, only a proportionally damped structure is considered. Actually, the hypothesis of proportional damping, although frequently considered in the literature, is difficult to be found in practice in real mechanical systems [23,24]. To study the feasibility and validity of the proposed method for relatively complex structural systems, we consider a linear 6-DOF railway vehicle model with viscous damping [25], as shown schematically in Figure 8. The mass matrix *M*, stiffness matrix *K*, and the damping matrix *C* of the system are given as follows [26]:

$$M = \begin{bmatrix} 1200 & 0 & 0 & 0 & 0 & 0 \\ 0 & 850 & 0 & 0 & 0 & 0 \\ 0 & 0 & 4125 & 0 & 0 & 0 \\ 0 & 0 & 0 & 125000 & 0 & 0 \\ 0 & 0 & 0 & 0 & 850 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1220 \end{bmatrix} \text{kg}$$

$$K = \begin{bmatrix} k\_1 + k\_2 & -k\_2 & 0 & 0 & 0 & 0 \\ -k\_2 & k\_2 + k\_3 & -k\_3 & -k\_3L & 0 & 0 \\ 0 & -k\_3 & k\_3 + k\_4 & k\_3L - k\_4L & -k\_4 & 0 \\ 0 & -k\_3L & k\_3L - k\_4L & k\_3L^2 + k\_4L^2 & k\_4L & 0 \\ 0 & 0 & -k\_4 & k\_4L & k\_4 + k\_5 & -k\_5 \\ 0 & 0 & 0 & 0 & -k\_5 & k\_5 + k\_6 \end{bmatrix} \text{N/m}$$

$$C = \begin{bmatrix} c\_1 + c\_2 & -c\_2 & 0 & 0 & 0 & 0 \\ -c\_2 & c\_2 + c\_3 & -c\_3 & -c\_3L & 0 & 0 \\ 0 & -c\_3 & c\_3 + c\_4 & c\_3L - c\_4L & -c\_4 & 0 \\ 0 & -c\_3L & c\_3L - c\_4L & c\_3L^2 + c\_4L^2 & c\_4L & 0 \\ 0 & 0 & -k\_4 & c\_4L & c\_4 + c\_5 & -c\_5 \\ 0 & 0 & 0 & 0 & -c\_5 & c\_5 + c\_6 \end{bmatrix} \text{N\cdot\'\'\'m}$$

in which *<sup>K</sup>*<sup>1</sup> = <sup>3</sup> × 107 N/m, *<sup>K</sup>*<sup>2</sup> = <sup>10</sup><sup>6</sup> N/m, *<sup>K</sup>*<sup>3</sup> = <sup>6</sup> × 106 N/m, *<sup>K</sup>*<sup>4</sup> = <sup>6</sup> × 106 N/m, *<sup>K</sup>*<sup>5</sup> = 106 N/m, *<sup>K</sup>*<sup>6</sup> = <sup>3</sup> × 107 N/m, *<sup>L</sup>* = 8.53 m, *<sup>c</sup>*<sup>1</sup> = 0, *<sup>c</sup>*<sup>2</sup> = <sup>6</sup> × 103 <sup>N</sup>·s/m, *<sup>c</sup>*<sup>3</sup> = 1.8 × <sup>10</sup><sup>4</sup> <sup>N</sup>·s/m, *<sup>c</sup>*<sup>4</sup> = 1.8 × 104 <sup>N</sup>·s/m, *<sup>c</sup>*<sup>5</sup> = <sup>6</sup> × 103 <sup>N</sup>·s/m, and *<sup>c</sup>*<sup>6</sup> = 0. This railway structural system has the features of relatively high modal damping levels for third and fourth modes, and a pair of closely spaced fifth and sixth modes. Note that this railway structure has non-proportional damping, because damping matrix, *C*, cannot be expressed as the linear combination of the mass matrix, *M*, and stiffness matrix, *K*. The stationary white noise in the previous numerical example is still used as the excitation force acting the 6th mass of the railway vehicle system. The sampling interval is set as 0.01 s, and the sampling period is 110 s. Then, using Newmark's method, the displacement responses obtained are employed for modal estimation using the modified SSI method proposed in this paper.

**Figure 8.** Schematic plot of the railway vehicle model.

It should be noted that the number of modes to be identified serves to determine the dimensions of the Hankel matrix. Using the channel expansion technique, we set the number of expansion channels at 20. The dimensions of the Hankel matrix are not less than the order of the system to be identified, thus satisfying the condition that the number of modes evaluated in the SSI algorithm is not less than the number of modes to be identified. The SVD analysis of the projection matrix, which is obtained from the Hankel matrix calculation, can then be implemented, and the number of singular values is employed to determine the order of the system to be identified. Around the 12th singular value, an obvious drop showed in the distribution of singular values associated with the data matrix [Ω][Ω] *<sup>T</sup>* constructed by projection matrix [Ω] from stationary response data, as shown in Figure 9. It can be estimated that the order of the system, i.e., the number of modes to be identified, is 12.

**Figure 9.** Distribution of the singular values associated with a data matrix, [Ω][Ω] *<sup>T</sup>*, constructed by projection matrix, [Ω], from stationary responses of a railway vehicle model..

Compared with the efficiency when performing the SVD of the system matrix [*A*], computation can be reduced by implementing SVD analysis of data matrix [Ω][Ω] *<sup>T</sup>* constructed from the projection matrix [Ω]. As shown in Figure 10, from the stabilization diagram corresponding to different modal orders, it can be observed that the number of structural modes to be identified is six, and serious modal interference can be observed between the last two close modes. The system matrix, [*A*], can be found using Equations (9)–(11), and the modal parameters can then be determined through Equations (12)–(15). Table 4 presents the well-implemented modal estimation through eigenvalue analysis of the system matrix, [*A*]. A comparison between the exact and identified mode shapes is shown in Figure 11, and the corresponding MAC values evaluated are shown in Figure 12, where good agreement is observed.

**Table 4.** Identification results of railway vehicle through Modified Stochastic Subspace Identification (MSSI) method.


**Figure 10.** Typical plot of the stabilization diagram of stationary responses of 1st DOF of a 6-DOF chain model of a railway vehicle.

**Figure 11.** Comparison between the identified mode shapes and the exact mode shapes of 6-DOF chain model of a railway vehicle subjected to stationary white noise: (**a**) 1st mode shape; (**b**) 2nd mode shape; (**c**) 3rd mode shape; (**d**) 4th mode shape; (**e**) 5th mode shape; (**f**) 6th mode shape.

**Figure 12.** Typical plot of the Modal Assurance Criterion (MAC) for the identified mode shapes and exact mode shapes of 6-DOF chain model of a railway vehicle.

Due to the experimental restrictions of economic cost and structural geometry in the practical measurement of structural response, the number of sensors sited on the structures to record response data is not usually sufficient for the overall degrees of freedom of a structure to be identified. This may cause problems of incomplete measurement of the degrees of freedom of the identified mode-shape vector. To address this issue, we also implement modal estimation from the incomplete modal-information data obtained due to insufficient measurement channels; thus, only response data of the first, third, and fifth DOF of the railway system subjected to ambient excitation in the previous numerical example are employed to implement modal estimation. The modal estimation results are shown in Table 5. The modal parameters estimated by the proposed method and the exact results are in good agreement, because the errors in both natural frequencies and damping ratios are less than 1%. Thus, we confirm the effectiveness of the proposed method under the likely practical conditions of insufficient measurement information.


**Table 5.** Identification results of railway vehicle through Modified Stochastic Subspace Identification (MSSI) method from incomplete modal-information measurement data.

#### *3.3. Experimental Validation of a Cantilever Beam*

To further validate the effectiveness of the method proposed in this paper, an actual beam structure of free boundary is used for the experiment, as shown in Figure 13. Brüel & Kjær RT Pro Photon 7.41 data acquisition system, PCB 208C02 force sensor (with a sensitivity of 112,410 mV/kN, a measurement range of 0.4448 kN, a low frequency response (−5%) of 0.001 Hz, and an upper frequency limit of 36,000 Hz), and Polytec OFV-5000 Modular Vibrometer (having a frequency range of DC-24 MHz with velocities up to ±10 m/s and displacements from the picometer to meter range) are used for measuring response signal. The simulated nonstationary excitation through Teledyne LeCroy T3AFG40 signal generator is imported into the Modal Shop K20070E01 vibration exciter, and the vibration shaker excites the cantilever beam structure. A typical white noise with an approximately consistent-power spectral density was synchronized with the shaker voltage time history and recorded, as shown in Figure 14, as provided from the waveform source panel in the RT Pro Photon 7.41 data acquisition system. Currently, an OMA-based roving different directions of sensor head of laser doppler vibrometer is performed to measure the actual modal properties of the beam structure by the data acquisition device, including four channels, only twenty measurement positions on the actual aluminum alloy beam were marked, as shown in Figure 15, and the shaker excitation impacts acted as the third location of the beam.

**Figure 13.** Schematic diagram of the experimental setup.

**Figure 14.** Power spectrum associated with stationary white noise.

**Figure 15.** Twenty measurement positions marked on the actual aluminum alloy beam.

The length, width, and height of the beam structure used for this paper are 512.8, 25.5, and 3.0 mm, respectively, the mass is 104.97 g, and the material is 5052-0 aluminum alloy. In addition, the "exact" modal parameters are obtained from the data of auto power spectrum of stationary response and cross power spectrum of stationary excitation and response by using ME'Scope software. Note that "exact" stationary excitation is a built-in stationary white noise generated from Brüel & Kjær RT Pro Photon 7.41 data acquisition system is imported into the Modal Shop K20070E01 electrodynamic shaker to excite the cantilever beam structure. The natural frequencies of the first four modes are about 10.85, 70.32, 176.31, and 290.31 Hz, as listed in Table 6. Finally, the modal parameters obtained are used to compare the identification results of MSSI, as shown in Table 6 and Figure 16. It is observed that the frequency errors are less than 15%, Figure 16 shows the identified mode shapes which are approximately coincident with "exact" mode shapes, and the MAC values are larger than 0.77. This means that the proposed method is effective on modal identification in practical application.

**Table 6.** Identification results of a practical cantilever beam through Modified Stochastic Subspace Identification (MSSI) method.


**Figure 16.** Comparison between the identified mode shapes and the "exact" mode shapes of a practical cantilever beam subjected to stationary white noise: (**a**) 1st mode shape; (**b**) 2nd mode shape; (**c**) 3rd mode shape; (**d**) 4th mode shape.

#### **4. Conclusions**

The topic of this paper was a study of ambient modal analysis based on the stochastic subspace identification technique (SSI). The paper aimed to develop the appropriate algorithms for output-only modal analysis to overcome difficulties when performing experimental modal analysis (EMA). As a modification of SSI, we introduced the procedure of solving the system matrix in SSI-COV in conjunction with SSI-DATA, allowing modal estimation to be well implemented. A system matrix can, therefore, be obtained directly from the observability matrix without evaluating the predictive-state matrix, and this will improve the efficiency of computation.

In addition, we extracted predictive-state matrixes with recursive relationships directly from the same original predictive-state matrix, and then omitted the step of reevaluating the predictive-state matrix at the next-time moment to improve the computational efficiency of the SSI method. In addition, through the SVD analysis of a data matrix [Ω][Ω] *<sup>T</sup>*, evaluated by the projection matrix, [Ω], the modal estimation can be effectively performed, and the corresponding computational efficiency can be improved.

By solving the system matrix through the observability matrix and constructing a new predictive-state matrix composed from the original measured data matrix, the procedure of modal estimation can be simplified, and the modal parameters can be effectively identified, even for a structural system having closely spaced modes and relatively high damping. Furthermore, the proposed modified SSI algorithm is applicable to the parametric estimation of structures with incomplete modal information obtained from insufficient measurement channels. In addition, the computational efficiency of the SSI method can be improved due to the non-uniqueness of the observability matrix. However, the need for white noise

excitation is still a main limitation to be resolved in the proposed method from ambient response, and many mechanical systems are expected to be excited by significantly different frequency content, in particular, by specific harmonics [27,28]. The actual limitations and the applicability of the proposed method to real mechanical systems could be considered for discussion in future work. Through numerical simulations and experimental verification, we illustrated and validated the effectiveness of the proposed method for modal estimation of structural systems from stationary ambient response data only.

**Author Contributions:** C.-S.L. conceived of the proposed idea, designed the experiments and verified the data along with Y.-X.W. The draft preparation was originally written by the first and corresponding author C.-S.L. The article was mainly revised by C.-S.L. and reviewed and edited in part by Y.-X.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by Ministry of Science and Technology of Taiwan under the grant MOST 110-2221-E-020-008-.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research was supported in part by Ministry of Science and Technology of Taiwan under the grant MOST 110-2221-E-020-008-. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions in revising the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

