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

Authors: Haotian Zhou; Kaiping Yu; Yushu Chen; Rui Zhao; Yunhe Bai
Journal: Shock and Vibration
 (2018)

---

## 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.


---

## Body

## 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 Section2, 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 ap-degree of freedom LTV system can be written as(1)Mtq¨t+C0tq˙t+Ktqt=B0ut, where M(t), C0(t), and K(t)∈Rp×p are mass, damping, and stiffness matrices of the system, respectively, q(t)∈Rp×1 is the displacement vector, B0∈Rp×ni is the input shape matrix, u(t)∈Rni×1 is the input vector, in which the notation Rn×1 denotes n-dimensional real vector space, Rn×n denotes the n×n real matrix space, and ni is the number of input force excitation signals.The corresponding output equation of the samep-degree of freedom LTV system can be expressed as(2)yt=Caq¨t+Cvq˙t+Cdqt,where y(t)∈Rno×1 is the output vector, in which no is the number of measured signals, and Ca, Cv, and Cd are the output matrices for the measurement of acceleration, velocity, and displacement, respectively.The system state vector can be defined as(3)xt=qtq˙t. And then (1) can be transformed into state-space motion equation, and the output equation can be written as follows:(4)x˙t=Atxt+Btut,yt=Ctxt+Dtut,where A(t)∈R2p×2p is the system matrix, B(t)∈R2p×ni is the input matrix, C(t)∈Rno×2p is the output matrix, and D(t)∈Rno×ni is the feedthrough matrix, which are given as(5)At=0I-M-1tKt-M-1tC0t,Bt=0M-1tB0,Ct=Cd-CaM-1tKtCv-CaM-1tC0t,Dt=CaM-1tB0.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(6)xk+1=Akxk+Bkuk+wnk,yk=Ckxk+vnk,where wn(k) is the process noise and vn(k) is the measurement noise. At time step k, the Hankel matrices formed by the I-O data can be written as(7)Uk=u1u2⋯uku2u3⋯uk+1⋮⋮⋱⋮uMuM+1⋯uk+M-1,Yk=y1y2⋯yky2y3⋯yk+1⋮⋮⋱⋮yMyM+1⋯yk+M-1,where M 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 n. Then we have(8)YkU⊥k=Yk-YkUTkUkUTk-1Uk=ΓkXkU⊥k,where(9)Γk=CkCk+1Ak⋮Ck+M-1Ak+M-2⋯Akis the generalized observability matrix, X(k) is the state-vector matrix, and U⊥(k) is defined as U⊥(k)=I-UT(k)U(k)UT(k)-1U(k). An estimation of Γ(k) is the principle subspace matrix of Y(k)U⊥(k). If Y(k)U⊥(k) can be obtained by a recursive form, the principle subspace tracking algorithm can be implicated to estimate Γ(k).The Hankel matrices at time stepk+1 can be written as(10)Uk+1=Uku-k+1,Yk+1=Yky-k+1,where(11)u-k+1=uTk+1uTk+2⋯uTk+MT,y-k+1=yTk+1yTk+2⋯yTk+MT,and then we have(12)Yk+1U⊥k+1Yk+1=Yk+1U⊥k+1Yk+1U⊥k+1T=YkU⊥kYk+pk+1pTk+1,where(13)pk+1=11+αk+1YkUTkUkUTk-1u-k+1-y-k+1,and α(k+1)=u-T(k+1)U(k)UT(k)-1u-(k+1). Data vector p(k) is fed into the subspace tracking algorithm to obtain the estimation of the generalized observability matrix Γ(k).Consider the following NIC framework for principal subspace analysis (PSA):(14)JW=12trlog⁡WTRppW-trWTW,where Rpp=EppT is the covariance matrix of random data vector p. The gradient of JW with respect to the signal subspace W is(15)∇JW=RppWWTRppW-1-W.Then the gradient ascent rule for updating W(k) is(16)Wk=1-ηWk-1+ηR^ppkWk-1WTk-1RppkWk-1,where η is the learning step. Because the exact covariance matrix Rpp is often unknown, here we use its estimation R^pp(k):(17)R^ppk=1k∑i=1kβk-ipipTi=βk-1R^ppk-1k+pkpTkk.Assume that WT(k-1)p(i) can be approximated by q(i)≜WT(i)p(i) for all 1≤i≤k; this approximation is also utilized in the PAST algorithm [9]. Substituting (17) into (16), we have(18)Wk=1-ηWk-1+η∑i=1kβk-ipiqTi∑i=1kβk-iqiqTi-1.Define P(k)=∑i=1kβk-iq(i)qT(i)-1 and Γ(k)=∑i=1kβk-ip(i)qT(i)P(k). Utilizing the matrix inversion lemma to P(k), we have(19)Pk=β∑i=1k-1βk-i-1qiqTi+qkqTk-1=βP-1k-1+qkqTk-1=Pk-1-gkqTkPk-1β,where g(k)=P(k-1)q(k)/β+qT(k)P(k-1)q(k) is called the gain vector. Using the property that g(k)qT(k)P(k-1)=g(k)qT(k)P(k-1)T=P(k-1)q(k)gT(k), we can write(20)Γk=β∑i=1k-1βk-i-1piqTi+pkqTkPk=∑i=1k-1βk-i-1piqTiPk-1-∑i=1k-1βk-i-1piqTiPk-1qkgTk+pkqTkPk=Γk-1+pk-Γk-1qkgTk.Now, we have the recursive form of Γ(k), and the tracking algorithm can be summarized as follows:Initializations P ( 0 ) = δ I r (δ is a small positive number; Ir is an r×r identity matrix). Γ ( 0 ) = 0 (n×r null matrix). W ( k ) is initialized by the SVD of Y(k)U⊥(k) described by (8).Updating Procedure (21) q k = W T k - 1 p k , g k = P k - 1 q k β + q T k P k - 1 q k , P k = P k - 1 - g k q T k P k - 1 β , Γ k = Γ k - 1 + p k - Γ k - 1 q k g T k , W k = 1 - η W k - 1 + η Γ k .A specific derivation of NIC algorithm can be seen in [21].


### 2.2. Time-Varying Modal Parameter Extraction

AfterΓ(k+1) is determined, the discrete system matrix can be formed as follows:(22)Ak=Γ1+k+1Γ2k=ΦkΛkΦTk,where Γ1(k+1) is the first M-1 block rows of matrix Γ(k+1); Γ2(k) is the last M-1 block of Γ(k); the superscript “+” represents Moore-Penrose inversion; Λ(k)=diag(λ1(k),…,λn(k)) is the eigenvalue matrix constructed by time-varying conjugate complex eigenvalues λi(k), in which λi(k)=exp⁡(-ςi(k)Δt±jωi(k)Δt) and j=-1; and Φ(k) represents the time-varying eigenvector matrix. Finally, the pseudo natural frequency and damping ratio can be obtained by the following equations:(23)ωik=ln⁡λik2Δtπ,ςik=-ln⁡λiRkΔtωik,where ωi(k) represents the ith order pseudo natural frequency and ςi(k) is the corresponding damping ratio and the superscript “R” represents the real part of the complex eigenvalue λi(k).


## 2.1. NIC Recursive Algorithm

The motion equation of ap-degree of freedom LTV system can be written as(1)Mtq¨t+C0tq˙t+Ktqt=B0ut, where M(t), C0(t), and K(t)∈Rp×p are mass, damping, and stiffness matrices of the system, respectively, q(t)∈Rp×1 is the displacement vector, B0∈Rp×ni is the input shape matrix, u(t)∈Rni×1 is the input vector, in which the notation Rn×1 denotes n-dimensional real vector space, Rn×n denotes the n×n real matrix space, and ni is the number of input force excitation signals.The corresponding output equation of the samep-degree of freedom LTV system can be expressed as(2)yt=Caq¨t+Cvq˙t+Cdqt,where y(t)∈Rno×1 is the output vector, in which no is the number of measured signals, and Ca, Cv, and Cd are the output matrices for the measurement of acceleration, velocity, and displacement, respectively.The system state vector can be defined as(3)xt=qtq˙t. And then (1) can be transformed into state-space motion equation, and the output equation can be written as follows:(4)x˙t=Atxt+Btut,yt=Ctxt+Dtut,where A(t)∈R2p×2p is the system matrix, B(t)∈R2p×ni is the input matrix, C(t)∈Rno×2p is the output matrix, and D(t)∈Rno×ni is the feedthrough matrix, which are given as(5)At=0I-M-1tKt-M-1tC0t,Bt=0M-1tB0,Ct=Cd-CaM-1tKtCv-CaM-1tC0t,Dt=CaM-1tB0.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(6)xk+1=Akxk+Bkuk+wnk,yk=Ckxk+vnk,where wn(k) is the process noise and vn(k) is the measurement noise. At time step k, the Hankel matrices formed by the I-O data can be written as(7)Uk=u1u2⋯uku2u3⋯uk+1⋮⋮⋱⋮uMuM+1⋯uk+M-1,Yk=y1y2⋯yky2y3⋯yk+1⋮⋮⋱⋮yMyM+1⋯yk+M-1,where M 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 n. Then we have(8)YkU⊥k=Yk-YkUTkUkUTk-1Uk=ΓkXkU⊥k,where(9)Γk=CkCk+1Ak⋮Ck+M-1Ak+M-2⋯Akis the generalized observability matrix, X(k) is the state-vector matrix, and U⊥(k) is defined as U⊥(k)=I-UT(k)U(k)UT(k)-1U(k). An estimation of Γ(k) is the principle subspace matrix of Y(k)U⊥(k). If Y(k)U⊥(k) can be obtained by a recursive form, the principle subspace tracking algorithm can be implicated to estimate Γ(k).The Hankel matrices at time stepk+1 can be written as(10)Uk+1=Uku-k+1,Yk+1=Yky-k+1,where(11)u-k+1=uTk+1uTk+2⋯uTk+MT,y-k+1=yTk+1yTk+2⋯yTk+MT,and then we have(12)Yk+1U⊥k+1Yk+1=Yk+1U⊥k+1Yk+1U⊥k+1T=YkU⊥kYk+pk+1pTk+1,where(13)pk+1=11+αk+1YkUTkUkUTk-1u-k+1-y-k+1,and α(k+1)=u-T(k+1)U(k)UT(k)-1u-(k+1). Data vector p(k) is fed into the subspace tracking algorithm to obtain the estimation of the generalized observability matrix Γ(k).Consider the following NIC framework for principal subspace analysis (PSA):(14)JW=12trlog⁡WTRppW-trWTW,where Rpp=EppT is the covariance matrix of random data vector p. The gradient of JW with respect to the signal subspace W is(15)∇JW=RppWWTRppW-1-W.Then the gradient ascent rule for updating W(k) is(16)Wk=1-ηWk-1+ηR^ppkWk-1WTk-1RppkWk-1,where η is the learning step. Because the exact covariance matrix Rpp is often unknown, here we use its estimation R^pp(k):(17)R^ppk=1k∑i=1kβk-ipipTi=βk-1R^ppk-1k+pkpTkk.Assume that WT(k-1)p(i) can be approximated by q(i)≜WT(i)p(i) for all 1≤i≤k; this approximation is also utilized in the PAST algorithm [9]. Substituting (17) into (16), we have(18)Wk=1-ηWk-1+η∑i=1kβk-ipiqTi∑i=1kβk-iqiqTi-1.Define P(k)=∑i=1kβk-iq(i)qT(i)-1 and Γ(k)=∑i=1kβk-ip(i)qT(i)P(k). Utilizing the matrix inversion lemma to P(k), we have(19)Pk=β∑i=1k-1βk-i-1qiqTi+qkqTk-1=βP-1k-1+qkqTk-1=Pk-1-gkqTkPk-1β,where g(k)=P(k-1)q(k)/β+qT(k)P(k-1)q(k) is called the gain vector. Using the property that g(k)qT(k)P(k-1)=g(k)qT(k)P(k-1)T=P(k-1)q(k)gT(k), we can write(20)Γk=β∑i=1k-1βk-i-1piqTi+pkqTkPk=∑i=1k-1βk-i-1piqTiPk-1-∑i=1k-1βk-i-1piqTiPk-1qkgTk+pkqTkPk=Γk-1+pk-Γk-1qkgTk.Now, we have the recursive form of Γ(k), and the tracking algorithm can be summarized as follows:Initializations P ( 0 ) = δ I r (δ is a small positive number; Ir is an r×r identity matrix). Γ ( 0 ) = 0 (n×r null matrix). W ( k ) is initialized by the SVD of Y(k)U⊥(k) described by (8).Updating Procedure (21) q k = W T k - 1 p k , g k = P k - 1 q k β + q T k P k - 1 q k , P k = P k - 1 - g k q T k P k - 1 β , Γ k = Γ k - 1 + p k - Γ k - 1 q k g T k , W k = 1 - η W k - 1 + η Γ k .A specific derivation of NIC algorithm can be seen in [21].


## 2.2. Time-Varying Modal Parameter Extraction

AfterΓ(k+1) is determined, the discrete system matrix can be formed as follows:(22)Ak=Γ1+k+1Γ2k=ΦkΛkΦTk,where Γ1(k+1) is the first M-1 block rows of matrix Γ(k+1); Γ2(k) is the last M-1 block of Γ(k); the superscript “+” represents Moore-Penrose inversion; Λ(k)=diag(λ1(k),…,λn(k)) is the eigenvalue matrix constructed by time-varying conjugate complex eigenvalues λi(k), in which λi(k)=exp⁡(-ςi(k)Δt±jωi(k)Δt) and j=-1; and Φ(k) represents the time-varying eigenvector matrix. Finally, the pseudo natural frequency and damping ratio can be obtained by the following equations:(23)ωik=ln⁡λik2Δtπ,ςik=-ln⁡λiRkΔtωik,where ωi(k) represents the ith order pseudo natural frequency and ςi(k) is the corresponding damping ratio and the superscript “R” represents the real part of the complex eigenvalue λi(k).


## 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,k-means, fuzzy C-means (FCM), and DBSCAN algorithms. Unlike k-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 pointp).
For any pointp in the sample dataset D, the spherical region with p as the center point and Eps as the radius is called the Eps-neighborhood of point p.Definition 2 (core point).
If the number of sample points within the Eps-neighborhood of a pointp is equal or larger than MinPts, than point p is called a core point.Definition 3 (directly density-reachable).
A pointp is directly density-reachable from the point q if p is within the Eps-neighborhood of q, and q is a core point.Definition 4 (density-reachable).
A pointp is density-reachable from point q if there is a point sequence p1,…,pn in data set D, and p1=q,pn=p, pi+1 is directly density-reachable from pi.Definition 5 (density-connected).
If there is a pointo that is density-reachable from both points p and q, then p and q 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 pointp is a core point, and O is a point set that is density-reachable from p, then O is a cluster.Lemma 8.
LetC be a cluster and point p be any point in it; then C equals the set of points which are density-reachable from p.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 datasetD through two steps: first, locate a core point as a seed from dataset D; 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 D: dataset to be clustered Output data: clustered datasets Procedure:
( 1 )
Find an arbitrary pointp in datasetD; check if it satisfies Definition 2. ( 2 )
IFp is a core point, THEN retrieve all points that are directly density-reachable from p, forming a cluster. ( 3 )
ELSE the selected point does not satisfy Definition2, THEN mark p as an interfering point, and exit this loop to find the next point. ( 4 )
REPEAT(1)–(4) UNTIL all the points in D are checked. ( 5 )
Retrieve all density-connected points for every directly density-reachable point within the Eps-neighborhood 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 o(nlog⁡n). 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.Figure 1
Flowchart of the proposed method.


## 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 ρ0 and the right half of that has a time-varying density ρ(t), which is described as(24)ρt=ρ0,0≤t≤tc,ρ01-βt-tc,tc<t≤tf,where tf indicates the full simulation duration, tc 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×1011 Pa, 0.3, and 7800 kg/m3, respectively; β, tc, and te 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.Figure 2
Schematic of a cantilever beam with a time-varying mass [14].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 3
Acceleration responses of the beam.
(a) (b) (c) (d) (e)Figure4 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(25)SNR=srssn,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 M=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 4
STFT analysis of the acceleration responses.
(a) (b) (c) (d) (e)Figure5(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.Figure 5
Reference results and estimated results by NIC algorithm with SNR = 10: (a) without the validation method; (b) with the validation method.
(a) (b)Figure 6
Reference results and estimated results by NIC algorithm with SNR = 0.5: (a) without the validation method; (b) with the validation method.
(a) (b)In order to achieve quantized results, the mean absolute percentage error (MAPE) is defined as(26)MAPE=1N∑i=1Nfi-f~ifi×100%,where fi denotes the true value and f~i denotes the identified value. N represents the total sample number.Table1 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%.Table 1
Error comparison of identified pseudo-natural frequencies.
SNR Tracking algorithm MAPE (%) 1st order 2nd order 3rd order Original Validated Original Validated Original Validated 30 PAST 1.75 1.42 0.24 0.24 0.62 0.62 API 1.75 1.45 0.24 0.24 0.62 0.62 NPI 2.50 1.45 0.24 0.24 0.62 0.62 NIC 1.81 1.49 0.24 0.24 0.62 0.62  10 PAST 1.92 1.61 0.24 0.24 0.62 0.62 API 1.95 1.66 0.25 0.25 0.75 0.62 NPI 2.70 1.67 0.35 0.25 0.81 0.63 NIC 1.95 1.62 0.24 0.24 0.62 0.62  1 PAST 4.36 3.67 0.68 0.43 2.15 0.67 API 6.10 5.06 0.61 0.50 1.95 0.73 NPI 6.48 3.88 0.54 0.52 2.16 0.64 NIC 3.89 3.57 0.69 0.43 1.96 0.68  0.8 PAST 4.40 4.10 0.85 0.46 3.53 0.69 API 11.29 6.77 1.76 0.61 2.95 0.72 NPI 10.48 4.20 1.82 0.66 3.77 0.75 NIC 4.17 3.85 1.19 0.51 3.35 0.70  0.5 PAST 9.04 5.60 1.43 0.71 4.21 0.74 API 23.07 11.8 2.81 1.07 4.11 0.78 NPI 23.80 6.32 3.00 1.01 5.58 0.63 NIC 7.78 5.73 1.22 0.73 3.71 0.72In 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 Figure7, 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.Figure 7
Schematic diagram of the test structure (mm).Figure 8
Schematic of the test setup.Figure 9
Temperature curves of the test structure.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 parameterM=400. 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.Figure 10
Identified result of NIC algorithm: (a) without the validation procedure; (b) with the validation procedure.
(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 Table1, 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.



---

*Source: 4378793-2018-03-27.xml*