Time-Varying Modal Parameters Identification by Subspace Tracking Algorithm and Its Validation Method

This article presents a time-varying modal parameter identification method based on the novel information criterion (NIC) algorithm and a post-processmethod for time-varyingmodal 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.


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 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 finite-datawindow 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.

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 M(), C 0 (), and K() ∈  × are mass, damping, and stiffness matrices of the system, respectively, q() ∈  ×1 is the displacement vector, B 0 ∈  ×  is the input shape matrix, u() ∈    ×1 is the input vector, in which the notation  ×1 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 as where y() ∈    ×1 is the output vector, in which   is the number of measured signals, and C  , C V , and C  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 A() ∈  2×2 is the system matrix, B() ∈  2×  is the input matrix, C() ∈    ×2 is the output matrix, and D() ∈    ×  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 as where w  () is the process noise and k  () is the measurement noise.At time step , the Hankel matrices formed by the I-O data can be written as where  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 have where is the generalized observability matrix, X() is the statevector matrix, and U ⊥ () is defined as can be obtained by a recursive form, the principle subspace tracking algorithm can be implicated to estimate Γ().
The Hankel matrices at time step  + 1 can be written as where and then we have where and ( + 1) = u  ( + 1)(U()U  ()) −1 u( + 1).Data vector p() is fed into the subspace tracking algorithm to obtain the estimation of the generalized observability matrix Γ().
Updating Procedure , A specific derivation of NIC algorithm can be seen in [21].

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 Epsneighborhood 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 Epsneighborhood of , and  is a core point.
Definition 4 (density-reachable).A point  is density-reachable from point  if there is a point sequence  1 , . . .,   in data set , and  1 = ,   = ,  +1 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 Epsneighborhood of a core point : dataset to be clustered Output data: clustered datasets Procedure: (1) Find an arbitrary point  in dataset D; check if it satisfies Definition 2. (2) IF  is a core point, THEN retrieve all points that are directly density-reachable from , forming a cluster.
(3) ELSE the selected point does not satisfy Definition 2, THEN mark  as an interfering point, and exit this loop to find the next point.(4) REPEAT (1)-( 4) UNTIL all the points in  are checked.
(5) Retrieve all density-connected points for every directly density-reachable point within the Epsneighborhood of all core points.(6) REPEAT (5)-( 6) 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 ( log ).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.

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  0 and the right half of that has a time-varying density (), which is described as where   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  0 , and Poisson's ratio are equal to 2.1 × 10 11 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 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.
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 as where srs and sn 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  = 0.96 and Hankel matrix block row parameter  = 100 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.
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 In order to achieve quantized results, the mean absolute percentage error (MAPE) is defined as where   denotes the true value and f 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.

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

Time (s)
1st-order estimated 2nd-order estimated 3rd-order estimated 1st-order calculated 2nd-order calculated 3rd-order calculated Frequency (Hz)  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 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   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.

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.

Figure 1 :
Figure 1: Flowchart of the proposed method.

Figure 4 :
Figure 4: STFT analysis of the acceleration responses.
(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.

Figure 5 :
Figure 5: Reference results and estimated results by NIC algorithm with SNR = 10: (a) without the validation method; (b) with the validation method.

Figure 6 :
Figure 6: Reference results and estimated results by NIC algorithm with SNR = 0.5: (a) without the validation method; (b) with the validation method.

Figure 7 :
Figure 7: Schematic diagram of the test structure (mm).

Table 1 :
Error comparison of identified pseudo-natural frequencies.