Research Article  Open Access
TimeVarying Modal Parameters Identification by Subspace Tracking Algorithm and Its Validation Method
Abstract
This article presents a timevarying modal parameter identification method based on the novel information criterion (NIC) algorithm and a postprocess method for timevarying modal parameter estimation. In the practical application of the timevarying modal parameter identification algorithm, the identified results contain both real modal parameters and aberrant ones caused by the measurement noise. In order to improve the quality of the identified results as well as sifting and validating the real modal parameters, a postprocess procedure based on densitybased spatial clustering of applications with noise (DBSCAN) algorithm is introduced. The efficiency of the proposed approach is first verified through a numerical simulation of a cantilever EulerBernoulli beam with a timevarying mass. Then the proposed approach is experimentally demonstrated by composite sandwich structure in a timevarying high temperature environment. The identified results illustrate that the proposed approach can obtain real modal frequencies in low signaltonoise ratio (SNR) scenarios.
1. Introduction
In the past few decades, the timevarying modal parameter identification method has attracted substantial interest in the modal analysis and testing community. The timevarying behavior is caused by the system parameters changing, such as mass properties of a flight vehicle and stiffness properties of a thermal protection system in the high temperature environment. In practice, many timevarying systems can be treated as the linear timevarying (LTV) system. Numerous methods have been developed under the LTV framework, such as timefrequency analysis method [1, 2]. Short time Fourier transform (STFT) has been used as an effective signal process tool in the modal analysis community, and the influence of Fourier analysis parameters has been studied in [3]. Compared to the Fourier analysis, the wavelet transform (WT) has an adjustable window size and location, and its effectiveness has been demonstrated in [4, 5].
Besides the timefrequency analysis method, the subspace identification method has also been widely used. Liu extended the concept of the modal parameters into LTV systems [6]. Subsequently, Liu and Deng [7] developed a subspacebased algorithm for pseudo natural frequencies (PNF) and experimentally verified it by an axially moving cantilever beam. The modal parameters, which are assumed to be changed slowly, can be extracted from a series of matrices constructed by the response data. Different from the above algorithm, Juang and Pappa [8] proposed a subspacebased identification algorithm that requires only one experiment. In this algorithm, the subspace is constructed by the singular value decomposition (SVD). For the online tracking, the SVD costs too much computational time. In order to overcome such shortage, Pang et al. [9] utilized the projection approximation subspace tracking (PAST) [10] technique instead of the SVD. Due to the problem of data saturation, PAST algorithm may lose its tracking ability. To solve this problem, Yu et al. [11] presented a finitedatawindow least squares technique. Goethals et al. [12] presented a recursive method to identify the vibroacoustic behavior of an airplane utilizing the instrumental variable projection approximation subspace tracking (IVPAST) technique. Xu et al. [13] transformed the firstorder statespace equations into linear equations by using the orthogonality of the wavelet scaling functions.
As stated in [11], the measurement noise may affect the accuracy of the identified results; several methods have been addressed to refine the accuracy. Qin et al. [15] addressed an improved stabilization diagram for stochastic subspace identification. The continuous wavelet transform (CWT) has a property of multiresolution, on which numerous methods are based [16, 17]. Ni et al. [18] improved the identification precision by applying the wavelet denoising technique, in which the effectiveness of this method is demonstrated by modal parameter identification of a spacecraft. Other methods focus on the postprocess procedure. Liu and Deng [7] developed a PNFs selection method based on the assumption that the PNFs vary slowly in a short time. The PNFs are selected by comparing with given values and the values are updated by the selected PNFs. Recently, Zhou et al. [2, 19] proposed an identification method based on the fuzzy clustering algorithm to improve the precision of the identified results. The goal of this method is to determine the membership of different orders, and the computational modes are separated from the physical modes. However, some parameters need to be artificially specified when using the aforementioned two methods. The method proposed by Liu and Deng requires the estimation of physical modal parameters, while the method proposed by Zhou et al. needs to specify the number of clusters.
The purpose of this paper is to improve the quality of the final identification results through the NIC based subspace tracking algorithm and the densitybased spatial clustering of applications with noise (DBSCAN) [20] algorithm. The sifting and validation approach proposed in this article can avoid specifying the number of clusters or the estimation of the physical modal parameters artificially.
The remainder of this paper is organized as follows. In Section 2, the modal parameter estimation method through NIC based subspace tracking algorithm is briefly introduced. In Section 3, DBSCAN algorithm and its application on modal parameter validation are introduced. In Section 4, a numerical simulation of a cantilever beam with timevarying mass is performed, and the identification results obtained by using NIC algorithm and other three different types of subspace tracking algorithms are discussed. In Section 5, an experimental example is performed to verify the effectiveness of the proposed method. Finally some conclusions are provided in Section 6.
2. Modal Parameter Extraction Based on Subspace Tracking
2.1. NIC Recursive Algorithm
The motion equation of a degree of freedom LTV system can be written as where , , and are mass, damping, and stiffness matrices of the system, respectively, is the displacement vector, is the input shape matrix, is the input vector, in which the notation denotes dimensional real vector space, denotes the real matrix space, and is the number of input force excitation signals.
The corresponding output equation of the same degree of freedom LTV system can be expressed aswhere is the output vector, in which is the number of measured signals, and , , and are the output matrices for the measurement of acceleration, velocity, and displacement, respectively.
The system state vector can be defined as And then (1) can be transformed into statespace motion equation, and the output equation can be written as follows:where is the system matrix, is the input matrix, is the output matrix, and is the feedthrough matrix, which are given as
Consider adding a noise term into (4) and noticing that the feedthrough matrix is the zero matrix in this study; then the statespace equation of the discrete timevarying system can be expressed aswhere is the process noise and is the measurement noise. At time step , the Hankel matrices formed by the IO data can be written aswhere is the block row number of the Hankel matrix, which must insure that the rank of the Hankel matrix is greater than the system order . Then we havewhereis the generalized observability matrix, is the statevector matrix, and is defined as . An estimation of is the principle subspace matrix of . If can be obtained by a recursive form, the principle subspace tracking algorithm can be implicated to estimate .
The Hankel matrices at time step can be written aswhereand then we havewhereand . Data vector is fed into the subspace tracking algorithm to obtain the estimation of the generalized observability matrix .
Consider the following NIC framework for principal subspace analysis (PSA):where is the covariance matrix of random data vector . The gradient of with respect to the signal subspace isThen the gradient ascent rule for updating iswhere is the learning step. Because the exact covariance matrix is often unknown, here we use its estimation :Assume that can be approximated by for all ; this approximation is also utilized in the PAST algorithm [9]. Substituting (17) into (16), we haveDefine and . Utilizing the matrix inversion lemma to , we havewhere is called the gain vector. Using the property that , we can writeNow, we have the recursive form of , and the tracking algorithm can be summarized as follows:
Initializations ( is a small positive number; is an identity matrix). ( null matrix). is initialized by the SVD of described by (8).
Updating Procedure
A specific derivation of NIC algorithm can be seen in [21].
2.2. TimeVarying Modal Parameter Extraction
After is determined, the discrete system matrix can be formed as follows:where is the first block rows of matrix ; is the last block of ; the superscript “+” represents MoorePenrose inversion; is the eigenvalue matrix constructed by timevarying conjugate complex eigenvalues , in which and ; and represents the timevarying eigenvector matrix. Finally, the pseudo natural frequency and damping ratio can be obtained by the following equations:where represents the th order pseudo natural frequency and is the corresponding damping ratio and the superscript “” represents the real part of the complex eigenvalue .
3. Validation Method Based on DBSCAN Algorithm
Clustering is a method that groups the data into different datasets (clusters); the elements in the same cluster have the familiar properties and are different from the elements in other clusters. Clustering method is different from classification method since the properties of the clusters are not given. Many clustering algorithms have been developed, for example, means, fuzzy means (FCM), and DBSCAN algorithms. Unlike means or FCM algorithm, DBSCAN algorithm is based on the density of the dataset and has the notion of noise. DBSCAN algorithm can determine arbitrarily shaped clusters from the dataset with noise. Application of the DBSCAN algorithm on modal parameter validation is based on one assumption that the number of identified physical modes is larger than that of the computational modes: that is, in the identified PNFs figure, the density of the physical modes area is larger than that of the computational modes area.
DBSCAN algorithm requires only two parameters: Eps and MinPts; here is the related definitions.
Definition 1 (Epsneighborhood of a point ). For any point in the sample dataset , the spherical region with as the center point and Eps as the radius is called the Epsneighborhood of point .
Definition 2 (core point). If the number of sample points within the Epsneighborhood of a point is equal or larger than MinPts, than point is called a core point.
Definition 3 (directly densityreachable). A point is directly densityreachable from the point if is within the Epsneighborhood of , and is a core point.
Definition 4 (densityreachable). A point is densityreachable from point if there is a point sequence in data set , and , is directly densityreachable from .
Definition 5 (densityconnected). If there is a point that is densityreachable from both points and , then and are called densityconnected.
Definition 6 (cluster). Given parameters Eps and MinPts, a cluster is the maximum nonempty set of densityconnected points based on the densityreachable property. Points that do not belong to any cluster are called noise.
According to the definitions given before, we can derive lemmas as follows.
Lemma 7. If point is a core point, and is a point set that is densityreachable from , then is a cluster.
Lemma 8. Let be a cluster and point be any point in it; then equals the set of points which are densityreachable from .
The goal of DBSCAN algorithm is to find the maximum set of densityreachable points. It can be seen from the two lemmas above that given parameters Eps and MinPts we can find a cluster from a dataset through two steps: first, locate a core point as a seed from dataset ; second, find all points that are densityreachable from the seed; then these points and the seed form a complete cluster. The detailed process of DBSCAN algorithm can be described as follows.
Description of DBSCAN Algorithm Input parameters and data: Eps: radius of Epsneighborhood MinPts: minimum number of points in the Epsneighborhood of a core point : dataset to be clustered Output data: clustered datasets Procedure:Find an arbitrary point in dataset D; check if it satisfies Definition 2.IF is a core point, THEN retrieve all points that are directly densityreachable from , forming a cluster.ELSE the selected point does not satisfy Definition 2, THEN mark as an interfering point, and exit this loop to find the next point.REPEAT – UNTIL all the points in are checked.Retrieve all densityconnected points for every directly densityreachable point within the Epsneighborhood of all core points.REPEAT  UNTIL the Epsneighborhood of all core points is processed. The remaining points that do not belong to any cluster are marked as noise.
According to [20], the average run time complexity of DBSCAN algorithm is . While the PNFs are too long, it is recommended to implement the DBSCAN algorithm on several subsections and then assemble the results of these subsections into a whole result. The major steps of the proposed method are shown in Figure 1.
4. Simulation Example and Performance Evaluation
In order to validate the proposed method, a cantilever EulerBernoulli beam with a timevarying mass is considered, which has been previously used in [22] by Zhao and Yu. As it is shown in Figure 2, the left half of the beam has a constant density and the right half of that has a timevarying density , which is described aswhere indicates the full simulation duration, indicates the time that the density begins to change, and is a constant value that relates to the speed of density variation. Here we assume that the relative velocity of the expelled mass with respect to the body coordinate system is zero. In this example, Young’s modulus, density , and Poisson’s ratio are equal to Pa, 0.3, and 7800 kg/m^{3}, respectively; , , and in (24) are 0.9, 3, and 10, respectively. The excitation force acting on the free end of the beam is a white Gaussian noise load with root mean square (RMS) equal to 1 N. The time step of the simulation is chosen to be 0.0005 s.
The beam model contains 31 nodes and 30 elements in total. Time finite element method for linear timevarying structures proposed by Zhao and Yu [14] is adopted to calculate the vibration responses of the beam. Acceleration responses of nodes 1, 7, 13, 19, and 25 are used for modal frequencies estimation, which are shown in Figure 3.
(a)
(b)
(c)
(d)
(e)
Figure 4 shows the spectrogram of the acceleration responses; it can be seen from the figure that there exist 4 active modes in the range of 0–1000 Hz. In order to demonstrate the antinoise ability of the proposed method, the response signals are assumed to be corrupted with zero mean Gaussian white noise, and the signaltonoise ratio (SNR) is defined aswhere and represent the standard deviation of the response signal and the added noise signal, respectively. The added noise can be treated as measurement noise. The parameters for all the tracking algorithms used in this example are provided: forgetting factor and Hankel matrix block row parameter which is elaborated in [23]. The data imported for DBSCAN is divided into 10 blocks, and each block is processed with the same parameters: MinPts = 100 and Eps = 0.002. Subspace tracking algorithms adopted in this example are PAST [10], approximated power iteration (API) [24], natural power iteration (NPI) [25], and the NIC algorithm introduced in Section 2.
(a)
(b)
(c)
(d)
(e)
Figure 5(a) shows the comparison of the reference results calculated by the theoretical modal analysis method and the results estimated by NIC algorithm with SNR = 10. Figure 6(a) shows the same comparison with SNR = 0.5. Black lines in these figures are the reference values and the cross marker points are the identified results sorted by magnitude. Since the response signals are contaminated with noise, the identified values between each pair of lines are inaccurate. Figures 5(b) and 6(b) demonstrate the improved results by the validation method introduced in Section 3. As shown in Figures 5(b) and 6(b), the aberrant points are removed, and the PNFs are automatically sifted.
(a)
(b)
(a)
(b)
In order to achieve quantized results, the mean absolute percentage error (MAPE) is defined aswhere denotes the true value and denotes the identified value. represents the total sample number.
Table 1 shows the comparison of the calculated corresponding true values with the values identified by four types of subspace tracking algorithms. The “original” represents the results sorted by magnitude, and “validated” indicates the results validated and sifted by the proposed method. For high SNR scenarios, the MAPEs of the NIC results are similar to the MAPEs of the PAST results, and the original results of the four algorithms agree well with the corresponding true values; therefore, the proposed method has little improvement. For low SNR scenarios, the NIC algorithm performs better than other three algorithms. Due to the antinoise ability limitation of the subspace tracking algorithms, the identified values fluctuate significantly, and the corresponding MAPEs increase while the SNR decreases. By implementing the proposed validation method, the MAPEs decrease observably; for example, the MAPE of the 1storder NPI result with SNR = 0.5 decreases from 23.80% to 6.32%.

In this section, we carry out the results identified by four subspace tracking algorithms with different level of noise. The results indicate that the NIC based subspace tracking algorithm and the validation method introduced in this paper can improve the accuracy of the final results, especially for low SNR scenarios.
5. Experimental Application
As mentioned in the preliminary section, the NIC based subspace tracking algorithm and the proposed validation method has shown improvement on identification precision by simulation results. In this section, we consider a timevarying structure in high temperature thermal environment. The object of this study is a sandwich structure, which is composed of fiberreinforced mullite matrix composite panels and elastic ceramic foams layer. As shown in Figure 7, the random excitation is applied on point 17. The excitation is emitted by a signal generator and can be treated as a white Gaussian noise with its RMS equal to 1 N. Five accelerometers are placed on points 1, 5, 21, 25 and 45. The schematic of experimental setups is shown in Figure 8. The signals are acquired in the heating process, and the temperature curves are shown in Figure 9. The sample frequency of the measured signals is 2000 Hz.
The temperature of the heated side increased from about 22°C to 800°C in about 120 seconds. The parameters used in NIC algorithm are as follows: forgetting factor is 0.94 and Hankel matrix block row parameter . Figure 10(a) shows the results identified by NIC algorithm, the PNFs are sorted by magnitude, and the color indicates the membership of different modes. Three major modes can be seen in Figure 10(a) located around 240, 350, and 480 Hz, and some values are not located around these lines; these points are computational modes. In Figure 10(b), most of the computational points are removed, the PNFs are grouped into three clusters, and each cluster represents a single mode. The variation of the PNFs is easier to observe in Figure 10(b). In the first half of the time history, the PNFs decrease in different degrees, while in the last half of time history, the PNFs increase obviously.
(a)
(b)
6. Conclusion
In this paper, a timevarying modal parameter identification method based on the NIC algorithm is introduced. A sifting and validation procedure for timevarying modal parameters is proposed. The results of numerical simulation demonstrate that the validation procedure introduced in this paper can improve the precision of the final result. Four types of subspace tracking algorithms are compared in this paper, and the proposed validation procedure is suitable for all these tracking methods. As shown in Table 1, when the SNR is too low (SNR < 10), the results of the NIC algorithm can reach the lowest MAPEs with the same calculation parameters used in other tracking algorithms. After the sifting and validation procedure, the inaccurate points are removed, the PNFs are automatically sorted, and the MAPEs of the final results are improved. The experimental example shows that the proposed method can identify PNFs of structures in high temperature and is effective in sifting and validating PNFs.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This present work was supported by the National Natural Science Foundation of China under Grant no. 11372084.
References
 S. Conforto and T. D. Alessio, “Spectral analysis for nonstationary signals from mechanical measurements: a parametric approach,” Mechanical Systems and Signal Processing, vol. 13, no. 3, pp. 395–411, 1999. View at: Publisher Site  Google Scholar
 S.D. Zhou, L. Liu, W. Yang, and Z.S. Ma, “Operational modal identification of timevarying structures via a vector multistage recursive approach in hybrid time and frequency domain,” Shock and Vibration, vol. 2015, Article ID 397364, 13 pages, 2015. View at: Publisher Site  Google Scholar
 M. Lamb and V. Rouillard, “Assessing the influence of fourier analysis parameters on shorttime modal parameter extraction,” Journal of Vibration and Acoustics, vol. 134, no. 3, Article ID 031002, 2012. View at: Publisher Site  Google Scholar
 R. Ghanem and F. Romeo, “A waveletbased approach for the identification of linear timevarying dynamical systems,” Journal of Sound and Vibration, vol. 234, no. 4, pp. 555–576, 2000. View at: Publisher Site  Google Scholar
 S. Nagarajaiah and B. Basu, “Output only modal identification and structural damage detection using time frequency & wavelet techniques,” Earthquake Engineering and Engineering Vibration, vol. 8, no. 4, pp. 583–605, 2009. View at: Publisher Site  Google Scholar
 K. Liu, “Extension of modal analysis to linear timevarying systems,” Journal of Sound and Vibration, vol. 226, no. 1, pp. 149–167, 1999. View at: Publisher Site  Google Scholar
 K. Liu and L. Deng, “Identification of pseudonatural frequencies of an axially moving cantilever beam using a subspacebased algorithm,” Mechanical Systems and Signal Processing, vol. 20, no. 1, pp. 94–113, 2006. View at: Publisher Site  Google Scholar
 J.N. Juang and R. S. Pappa, “An eigensystem realization algorithm for modal parameter identification and model reduction,” Journal of Guidance, Control, and Dynamics, vol. 8, no. 5, pp. 620–627, 1985. View at: Publisher Site  Google Scholar
 S. Pang, K. Yu, and J. Zou, “Improved subspace method with application in linear tievarying structural modal parameter identification,” Chinese Journal of Applied Mechanics, vol. 22, no. 2, pp. 184–189, 2005. View at: Google Scholar
 B. Yang, “Projection approximation subspace tracking,” IEEE Transactions on Signal Processing, vol. 43, no. 1, pp. 95–107, 1995. View at: Publisher Site  Google Scholar
 K. Yu, K. Yang, and Y. Bai, “Experimental investigation on the timevarying modal parameters of a trapezoidal plate in temperaturevarying environments by subspace trackingbased method,” Journal of Vibration and Control, vol. 21, no. 16, pp. 3305–3319, 2014. View at: Publisher Site  Google Scholar
 I. Goethals, L. Mevel, A. Benveniste, and B. D. Moor, “Recursive outputonly subspace identification for inflight flutter monitoring,” in Proceedings of the In Proceedings of the 22nd International Modal Analysis Conference, vol. 7, Dearborn, Mich, USA, 2004. View at: Google Scholar
 X. Xu, Z. Y. Shi, and Q. You, “Identification of linear timevarying systems using a waveletbased statespace method,” Mechanical Systems and Signal Processing, vol. 26, no. 1, pp. 91–103, 2012. View at: Publisher Site  Google Scholar
 R. Zhao and K. Yu, “Hamilton's law of variable mass system and time finite element formulations for timevarying structures based on the law,” International Journal for Numerical Methods in Engineering, vol. 99, no. 10, pp. 711–736, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 S. Qin, J. Kang, and Q. Wang, “Operational Modal Analysis Based on Subspace Algorithm with an Improved Stabilization Diagram Method,” Shock and Vibration, vol. 2016, Article ID 7598965, 10 pages, 2016. View at: Publisher Site  Google Scholar
 T. Kijewski and A. Kareem, “Wavelet transforms for system identification in civil engineering,” ComputerAided Civil and Infrastructure Engineering, vol. 18, no. 5, pp. 339–355, 2003. View at: Publisher Site  Google Scholar
 B. Yan and A. Miyamoto, “A comparative study of modal parameter identification based on wavelet and HilbertHuang transforms,” ComputerAided Civil and Infrastructure Engineering, vol. 21, no. 1, pp. 9–23, 2006. View at: Publisher Site  Google Scholar
 Z. Ni, R. Mu, G. Xun, and Z. Wu, “Timevarying modal parameters identification of a spacecraft with rotating flexible appendage by recursive algorithm,” Acta Astronautica, vol. 118, pp. 49–61, 2016. View at: Publisher Site  Google Scholar
 S.D. Zhou, W. Heylen, P. Sas, and L. Liu, “Parametric modal identification of timevarying structures and the validation approach of modal parameters,” Mechanical Systems and Signal Processing, vol. 47, no. 12, pp. 94–119, 2014. View at: Publisher Site  Google Scholar
 M. Ester, H.P. Kriegel, J. Sander, and X. Xu, “A densitybased algorithm for discovering clusters in large spatial databases with noise,” in Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining (KDD '96), pp. 226–231, 1996. View at: Google Scholar
 Y. Miao and Y. Hua, “Fast subspace tracking and neural network learning by a novel information criterion,” IEEE Transactions on Signal Processing, vol. 46, no. 7, pp. 1967–1979, 1998. View at: Publisher Site  Google Scholar
 R. Zhao and K. Yu, “An efficient transient analysis method for linear timevarying structures based on multilevel substructuring method,” Computers & Structures, vol. 146, pp. 76–90, 2015. View at: Publisher Site  Google Scholar
 F. Tasker, A. Bosse, and S. Fisher, “Realtime modal parameter estimation using subspace methods: Theory,” Mechanical Systems and Signal Processing, vol. 12, no. 6, pp. 797–808, 1998. View at: Publisher Site  Google Scholar
 R. Badeau, B. David, and G. Richard, “Fast approximated power iteration subspace tracking,” IEEE Transactions on Signal Processing, vol. 53, no. 8, part 1, pp. 2931–2941, 2005. View at: Publisher Site  Google Scholar  MathSciNet
 Y. Hua, Y. Xiang, T. Chen, K. AbedMeraim, and Y. Miao, “New look at the power method for fast subspace tracking,” Digital Signal Processing, vol. 9, no. 4, pp. 297–314, 1999. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Haotian Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.