Abstract
This article presents a time-varying modal parameter identification method based on the novel information criterion (NIC) algorithm and a post-process method for time-varying modal parameter estimation. In the practical application of the time-varying 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 post-process procedure based on density-based 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 Euler-Bernoulli beam with a time-varying mass. Then the proposed approach is experimentally demonstrated by composite sandwich structure in a time-varying high temperature environment. The identified results illustrate that the proposed approach can obtain real modal frequencies in low signal-to-noise ratio (SNR) scenarios.
1. Introduction
In the past few decades, the time-varying modal parameter identification method has attracted substantial interest in the modal analysis and testing community. The time-varying 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 time-varying systems can be treated as the linear time-varying (LTV) system. Numerous methods have been developed under the LTV framework, such as time-frequency 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 time-frequency 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 subspace-based 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 subspace-based 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 finite-data-window 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 (IV-PAST) technique. Xu et al. [13] transformed the first-order state-space 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 post-process 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 density-based 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 time-varying 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 state-space 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 state-space equation of the discrete time-varying system can be expressed aswhere is the process noise and is the measurement noise. At time step , the Hankel matrices formed by the I-O 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 state-vector 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. Time-Varying 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 Moore-Penrose inversion; is the eigenvalue matrix constructed by time-varying conjugate complex eigenvalues , in which and ; and represents the time-varying 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 (Eps-neighborhood 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 Eps-neighborhood of point .
Definition 2 (core point). If the number of sample points within the Eps-neighborhood of a point is equal or larger than MinPts, than point is called a core point.
Definition 3 (directly density-reachable). A point is directly density-reachable from the point if is within the Eps-neighborhood of , and is a core point.
Definition 4 (density-reachable). A point is density-reachable from point if there is a point sequence in data set , and , is directly density-reachable from .
Definition 5 (density-connected). If there is a point that is density-reachable from both points and , then and are called density-connected.
Definition 6 (cluster). Given parameters Eps and MinPts, a cluster is the maximum nonempty set of density-connected points based on the density-reachable 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 density-reachable 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 density-reachable from .
The goal of DBSCAN algorithm is to find the maximum set of density-reachable 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 density-reachable 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 Eps-neighborhood MinPts: minimum number of points in the Eps-neighborhood 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 density-reachable 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 density-connected points for every directly density-reachable point within the Eps-neighborhood of all core points.REPEAT - UNTIL the Eps-neighborhood 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 Euler-Bernoulli beam with a time-varying 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 time-varying 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/m3, 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 time-varying 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 signal-to-noise 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 1st-order 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 time-varying structure in high temperature thermal environment. The object of this study is a sandwich structure, which is composed of fiber-reinforced 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 time-varying modal parameter identification method based on the NIC algorithm is introduced. A sifting and validation procedure for time-varying 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.