# Instantaneous Modal Parameter Identification of Linear Time-Varying Systems Based on Chirplet Adaptive Decomposition

Authors: Jie Zhang; Zhiyu Shi
Journal: Shock and Vibration
 (2019)

---

## Abstract

Instantaneous modal parameter identification of time-varying dynamic systems is a useful but challenging task, especially in the identification of damping ratio. This paper presents a method for modal parameter identification of linear time-varying systems by combining adaptive time-frequency decomposition and signal energy analysis. In this framework, the adaptive linear chirplet transform is applied in time-frequency analysis of acceleration response for its higher energy concentration, and the response of each mode can be adaptively decomposed via an adaptive Kalman filter. Then, the damping ratio of the time-varying systems is identified based on energy analysis of component response signal. The proposed method can not only improve the accuracy of instantaneous frequency extraction but also ensure the antinoise ability in identifying the damping ratio. The efficiency of the method is first verified through a numerical simulation of a three-degree-of-freedom time-varying structure. Then, the method is demonstrated by comparing with the traditional wavelet and time-domain peak method. The identified results illustrate that the proposed method can obtain more accurate modal parameters in low signal-to-noise ratio (SNR) scenarios.


---

## Body

## 1. Introduction

System identification techniques have received significant attention in civil, aerospace, and mechanical engineering in the past few decades [1]. Most of them assume that the analyzed systems are linear time-invariant (LTI), i.e., the output for such systems does not change with a delay in the input. However, in reality, this assumption is not valid for many engineering structures which are damaged or subjected to strong excitation. Thus conventional modal analysis methods and experimental techniques are not appropriate for these linear time-varying (LTV) systems [2]. In such cases, for the purpose of health monitoring or damage identification, it is necessary to develop reliable approaches for the identification of time-varying characteristics.In terms of LTV system identification, many methods have been developed for various time-varying cases. These methods generally belong to two categories: time-domain methods and time-frequency-domain methods [3].In the time domain, identification methods based on the time-varying model such as state-space model and time series model have become the main research objects since the time-domain description of LTV systems is very convenient. Liu [4] proposed a subspace-based identification technique and used eigenvalues of the estimated transition matrices to characterize the dynamic properties. The subspace is constructed by the singular value decomposition (SVD) which costs too much computational time. Authors in [5–7] further presented improved algorithm to overcome such drawback. Pang et al. [7] utilized the projection approximation subspace tracking (PAST) technique instead of the SVD. Zhou et al. [8] presented a postprocess method to improve the quality of identification results in low signal-to-noise ratio scenarios. Meanwhile, parameter identification is established based on the autoregressive-moving-average (ARMA) model which can match the structures of real-world processes. In order to process nonstationary signals and identify the parameters of time-varying structures, many approaches were also proposed based on the time-dependent autoregressive moving average (TARMA) model [9]. The TARMA model can be utilized as an efficient tool to reveal the dynamic characteristics of nonstationary signals due to its simplicity and effectiveness. However, a continuing challenge is how to establish a precise TARMA model. Typically, the recursive least-squares (RLS) approach and the Kalman filter [10] were applied to solve the problem. Recently, Spiridonakos and Fassois [11] presented an adaptable functional series TARMA model for nonstationary signal modelling, which employs the B-spline function to model the TARMA coefficients adaptively and structurally. Yang et al. [12] developed a moving Kriging shape function modelling of vector TARMA models, which has the advantage on the sudden change.Furthermore, time-frequency (TF) estimators have been more extensively studied, i.e., by using the Hilbert transform (HT) and wavelet transform (WT), since the response of LTV systems is a typical nonstationary signal which needs to be processed by time-frequency analysis (TFA) tools. Meanwhile, TFA tools become effective in solving time-varying equations with the development of TF theory. Shi et al. [13] firstly addressed the Hilbert transform and empirical mode decomposition (EMD) method to identify time-varying multi-degree-of-freedom (MDOF) systems. Ghanem and Romeo [14] presented a wavelet-based identification method, which projects the governing differential equation on a subspace spanned by a finite number of wavelets and identifies the parameters using the Galerkin approach. Furthermore, a wavelet-based state-space method for LTV system identification was proposed by Xu et al. [15]. Hou et al. [16, 17] proposed a continuous wavelet-based technique for identification of instantaneous modal parameters of time-varying structures. Dziedziech et al. [18] used the wavelet-based frequency response function (FRF) for modal identification of time-variant systems. An identification approach based on the empirical wavelet transform (EWT) was developed for LTV systems by Xin et al. [19].As seen from the above analysis, conventional TFA methods are suitable for processing such nonstationary signals of LTV systems. However, they are insufficient to provide a TF representation with high resolution due to their inherent restrictions. HT can process multicomponent signals only when combined with EMD which always has the shortage of modal aliasing and end swing phenomena. For short-time TFA tools like WT [20, 21], they are established on the assumption of the considered signal being piecewise stationary in a short time. This means that they essentially use horizontal lines to approximate the IF feature. To improve the performance of the TFA methods, many effective methods have been proposed. One relatively effective strategy named synchrosqueezing transform (SST) [22] is to squeeze the TF coefficients into the instantaneous frequency trajectory. However, the strongly modulated frequency component and noise could destroy the original TF representation of the signal heavily. Another strategy named parameterized TFA was inspired by the fact that using adaptive window function (such as adding chirp rate in base function) in WT can achieve a much higher TF resolution. Thus the chirplet transform (CT) was designed by Deng et al. [23]. By using an extra parameter, chirp rate, the CT is able to create a concentrated TF representation for the linear modulated signal. In order to search the best chirp rate according to the TF signature of nonstationary signal, Yu and Zhou [24] developed an adaptive algorithm named GLCT. In terms of the TF characteristics of vibration response, the responses of LTV systems are usually characterized by multicomponent amplitude-modulation (AM) and frequency-modulation (FM) signals. Therefore, the CT is particularly suitable for identifying the instantaneous frequency of time-varying signals. In addition, for nonstationary signals with noise, the Vold–Kalman (VK) filter can more effectively decompose the multicomponent signals [25] compared with EMD.This paper presents a time-frequency decomposition method to extract the individual modes from the free response signal of LTV systems. The IF is extracted based on GLCT as the complex carrier wave of the signal, and then, the VK filter is constructed to perform signal decomposition adaptively. Using this method, the IF of each time-varying component can be identified more accurately, and the amplitude envelope of each signal component can be also obtained. Thus many damping ratio identification methods based on amplitude envelope were proposed such as the time-domain method [26, 27]. The identification of damping ratio is of great value for many engineering applications such as flutter analysis of aircraft. But this task is of challenges since the amplitude envelope is easily disturbed by noise which leads to a much lower identification accuracy of damping ratio than frequency. Considering this shortcoming, the improvement based on energy analysis [28] is explored to make the identification method insensitive to noise. The damping ratio is identified by energy ratio of adjacent half-periodic signals, which can effectively eliminate the effect of noise on amplitude. The aim of this paper is to improve the identification accuracy of instantaneous modal parameters, i.e., instantaneous frequencies and damping ratios.With this prime objective, this paper is organized as follows: the theory of time-frequency decomposition based on GLCT and VK filter is briefly introduced in Section2. The energy analysis for damping ratio identification of LTV systems is presented in Section 3. Section 4 shows a series of comparative numerical simulation examples to illustrate the performance of identification. In these comparative examples, the Morlet wavelet and time-domain peak method are selected for comparison to identify instantaneous frequencies and damping ratios, respectively. Finally, Section 5 provides the discussion and conclusions on the identification results.


## 2. Theoretical Background

Using the adaptive TF method for vibration signal decomposition consists of two main steps: (1) performing TFA of the target vibration signals by using GLCT to extract the TF ridge of each time-varying component and (2) constructing the filtering band based on the extracted TF ridge, then decomposing the target vibration signals into multicomponent AM and FM signals to obtain the amplitude envelope. To apply this TF method for LTV systems, we can identify the instantaneous frequencies based on the TF ridge and identify the damping ratio based on the amplitude envelope, respectively.


### 2.1. General Linear Chirplet Transform (GLCT)

GLCT is an approach to obtain the best demodulated operator for signals with time-varying IF, which has better performance than any individual CT in generating TF representation.According to the theory of Ville, an analytic signal with time-varying IF can be written as(1)st=atej∫φtdt,where at is the amplitude and φt is the IF of the signal.In a short timeτ, φt can be regarded as a linear equation approximately, which can be expanded by the Taylor expansion:(2)φτ=φt+φ′tτ−t,where φt is the IF at time t and φ′t is the first-order derivative of IF at time t.The CT of signalst which considers the time-variant demodulated operator can be written as(3)CTt,ω,ct=∫−∞+∞sτ⋅wτ−t⋅e−jωτ⋅e−jctτ−t2/2dτ,where wτ−t is a window function and e−jctτ−t2/2 is the time-variant demodulated operator.The CT amplitude in IFφt (ω=φt) can be calculated by(4)CTt,ω,ct=∫−∞+∞wτ−t⋅aτ⋅ejφττ+jφ′tτ−t2/2⋅e−jφtτ⋅e−jctτ−t2/2dτ=∫−∞+∞ejφ′t−ctτ−t2/2⋅wτ−t⋅aτ⋅dτ≤∫−∞+∞wτ−taτdτ.From equation (4), it can be seen that, when ct in CT is consistent with φ′t of the analyzed signal, the CT amplitude will achieve the maximum, and then the better TF representation with high energy concentration will be obtained. But for multicomponent signals, the IF feature cannot be known in advance, which leads to that it is hard to determine the demodulated operator accurately.A feasible solution is to use a series of discrete demodulated operators to approximate the best demodulated operator [24]. Actually the demodulated operator is also a rotating operator of signals on the TF result, and the rotating degree is arctan−c. So a parameter α is introduced to characterize the parameter c as(5)α=arctan2TsFsc,α∈−π2,−π2,where Ts is the sampling time and Fs is the sampling frequency.In a discrete way, the parameterα is divided into Nα values on the TF plane from −π/2 to π/2 as(6)α=−π2+πNα+1,−π2+2πNα+1,…,−π2+NαπNα+1.Then, equation (3) can be rewritten as(7)CTt,ω,ct=∫−∞+∞sτ⋅wτ−t⋅e−jωτ⋅e−jtanαFs/2Tsτ−t2/2dτ.From this, the TF spectrum with high energy concentration can be defined as(8)GSpect,ω=CTt,ω,cm2,where cm is a discrete value when the CT amplitude achieves the maximum in TF point t,ω, which can be calculated by(9)cm=argmaxCTt,ω,ct.In order to demonstrate the better TF representation of GLCT, WT and GLCT are used for TF analysis of a time-varying signal, respectively. The results are shown in Figure1. It is shown that GLCT can adaptively track the instantaneous frequency trajectory at different TF points. Compared with WT, GLCT uses adaptive window function instead of rectangular TF window to achieve a much higher TF resolution. Thus, this approach can more adaptively and accurately estimate IF than traditional TFA tools especially for signals with time-varying IF.Figure 1
TF representation for a time-varying signal: (a) using WT and (b) using GLCT.
(a)
(b)


### 2.2. Vold–Kalman Filter

In terms of nonstationary signal decomposition, various interesting approaches have been developed such as EMD [13], empirical wavelet transform (EWT) [19], and Kalman filter [29]. But considering a multicomponent signal with time-varying IF in low SNR scenarios, the VK filter is more adapted for decomposing this signal [30]. Since the VK filter needs an estimation of the IF firstly, a novel time-frequency decomposition method is presented when GLCT is applied for IF extraction. In this section, the decomposition process is expounded in detail.Consider a multicomponent analytic signal composed of component signal in equation (1) with random noise which can be written as(10)st=∑k∈iaktΦkt+μt=∑k∈iaktej∫φktdt+μti=1,2,…,K,where K is the number of component signals and μt is a zero mean white noise.In a discrete way, the signal is rewritten as(11)sn=∑k∈iaknΦkn+μn,i=1,2,…,K,where n=1,2,…,N and N is the total number of sampling points.By using GLCT in TFA of this multicomponent signal, the IF of each component can be estimated. Thus, the complex-valued carrier waveΦkn can be constructed.


#### 2.2.1. Structural Equation

It is assumed that the amplitude envelopeakt is the low-frequency modulation of the complex carrier wave Φkt. In other words, envelope is locally approximated by a low-order polynomial, and it is filtered in some specific way. This condition can be expressed by a structural equation with the nonhomogeneity term εkn. The polynomial order designates the number of the filter poles. The equations for l-pole filter are given by(12)∇lakn=εkn.Takingl=2 as an example, the corresponding equation is given by(13)∇2akn=akn−1−2akn+akn+1=εkn.The system of the structural equations containing all the samples takes the following form, which can be rewritten in the matrix form:(14)−210⋯001−21⋯0001−2⋯00⋮⋮⋮⋱⋮⋮000⋯1−2ak1ak2ak3⋮akN=εk1εk2εk3⋮εkN,(15)GAk=ε.


#### 2.2.2. Data Equation

Taking the measurement noise into account, the measurement equation is written in the matrix form:(16)S−W1W2W3⋯WKA1A2A3⋮AK=U,where S=s1s2s3⋮sN, Wk=Φk10⋯00Φk2⋯0⋮⋮⋱⋮00⋯ΦkN, Ak=ak1ak2⋮akN, U=μ1μ2μ3⋮μN, and N is the number of sampling points.


#### 2.2.3. Global Solution

The system of equations (15) and (16) is an underdetermined system for the unknown amplitude envelope akt. The additional condition for the equation solution is that the variances of the nonhomogeneity terms εkn and noise μn have to be minimal while maintaining the given relationship between them.The weighted sum of equations (15) and (16) gives the loss function:(17)J=∑k∈ir2εTε+UTU=∑k∈ir2AkHGTGAk+ST−∑k∈iAkHWkHS−∑k∈iWkAk,where r is a weighting factor. The choice of a large value for the weighting factor r leads to the highly selective filtration in the frequency domain. In contrast, fast convergence with low-frequency resolution is achieved by choosing the small r value. The limit value of r [30] is shown in Table 1.Table 1
The value of the weighting factor.
Pole number
s
=
1
s
=
2
s
=
3
s
=
4
The max value ofr
7
×
10
6
4
×
10
6
2
×
10
6
1.1
×
10
6The first derivative of equation (17) with respect to the vector Ak gives a condition for the minimum of this function, which is called a normal equation:(18)∂J∂AkH=0.Thus the unknown envelopeAk results from equation (18).


## 2.1. General Linear Chirplet Transform (GLCT)

GLCT is an approach to obtain the best demodulated operator for signals with time-varying IF, which has better performance than any individual CT in generating TF representation.According to the theory of Ville, an analytic signal with time-varying IF can be written as(1)st=atej∫φtdt,where at is the amplitude and φt is the IF of the signal.In a short timeτ, φt can be regarded as a linear equation approximately, which can be expanded by the Taylor expansion:(2)φτ=φt+φ′tτ−t,where φt is the IF at time t and φ′t is the first-order derivative of IF at time t.The CT of signalst which considers the time-variant demodulated operator can be written as(3)CTt,ω,ct=∫−∞+∞sτ⋅wτ−t⋅e−jωτ⋅e−jctτ−t2/2dτ,where wτ−t is a window function and e−jctτ−t2/2 is the time-variant demodulated operator.The CT amplitude in IFφt (ω=φt) can be calculated by(4)CTt,ω,ct=∫−∞+∞wτ−t⋅aτ⋅ejφττ+jφ′tτ−t2/2⋅e−jφtτ⋅e−jctτ−t2/2dτ=∫−∞+∞ejφ′t−ctτ−t2/2⋅wτ−t⋅aτ⋅dτ≤∫−∞+∞wτ−taτdτ.From equation (4), it can be seen that, when ct in CT is consistent with φ′t of the analyzed signal, the CT amplitude will achieve the maximum, and then the better TF representation with high energy concentration will be obtained. But for multicomponent signals, the IF feature cannot be known in advance, which leads to that it is hard to determine the demodulated operator accurately.A feasible solution is to use a series of discrete demodulated operators to approximate the best demodulated operator [24]. Actually the demodulated operator is also a rotating operator of signals on the TF result, and the rotating degree is arctan−c. So a parameter α is introduced to characterize the parameter c as(5)α=arctan2TsFsc,α∈−π2,−π2,where Ts is the sampling time and Fs is the sampling frequency.In a discrete way, the parameterα is divided into Nα values on the TF plane from −π/2 to π/2 as(6)α=−π2+πNα+1,−π2+2πNα+1,…,−π2+NαπNα+1.Then, equation (3) can be rewritten as(7)CTt,ω,ct=∫−∞+∞sτ⋅wτ−t⋅e−jωτ⋅e−jtanαFs/2Tsτ−t2/2dτ.From this, the TF spectrum with high energy concentration can be defined as(8)GSpect,ω=CTt,ω,cm2,where cm is a discrete value when the CT amplitude achieves the maximum in TF point t,ω, which can be calculated by(9)cm=argmaxCTt,ω,ct.In order to demonstrate the better TF representation of GLCT, WT and GLCT are used for TF analysis of a time-varying signal, respectively. The results are shown in Figure1. It is shown that GLCT can adaptively track the instantaneous frequency trajectory at different TF points. Compared with WT, GLCT uses adaptive window function instead of rectangular TF window to achieve a much higher TF resolution. Thus, this approach can more adaptively and accurately estimate IF than traditional TFA tools especially for signals with time-varying IF.Figure 1
TF representation for a time-varying signal: (a) using WT and (b) using GLCT.
(a)
(b)


## 2.2. Vold–Kalman Filter

In terms of nonstationary signal decomposition, various interesting approaches have been developed such as EMD [13], empirical wavelet transform (EWT) [19], and Kalman filter [29]. But considering a multicomponent signal with time-varying IF in low SNR scenarios, the VK filter is more adapted for decomposing this signal [30]. Since the VK filter needs an estimation of the IF firstly, a novel time-frequency decomposition method is presented when GLCT is applied for IF extraction. In this section, the decomposition process is expounded in detail.Consider a multicomponent analytic signal composed of component signal in equation (1) with random noise which can be written as(10)st=∑k∈iaktΦkt+μt=∑k∈iaktej∫φktdt+μti=1,2,…,K,where K is the number of component signals and μt is a zero mean white noise.In a discrete way, the signal is rewritten as(11)sn=∑k∈iaknΦkn+μn,i=1,2,…,K,where n=1,2,…,N and N is the total number of sampling points.By using GLCT in TFA of this multicomponent signal, the IF of each component can be estimated. Thus, the complex-valued carrier waveΦkn can be constructed.


### 2.2.1. Structural Equation

It is assumed that the amplitude envelopeakt is the low-frequency modulation of the complex carrier wave Φkt. In other words, envelope is locally approximated by a low-order polynomial, and it is filtered in some specific way. This condition can be expressed by a structural equation with the nonhomogeneity term εkn. The polynomial order designates the number of the filter poles. The equations for l-pole filter are given by(12)∇lakn=εkn.Takingl=2 as an example, the corresponding equation is given by(13)∇2akn=akn−1−2akn+akn+1=εkn.The system of the structural equations containing all the samples takes the following form, which can be rewritten in the matrix form:(14)−210⋯001−21⋯0001−2⋯00⋮⋮⋮⋱⋮⋮000⋯1−2ak1ak2ak3⋮akN=εk1εk2εk3⋮εkN,(15)GAk=ε.


### 2.2.2. Data Equation

Taking the measurement noise into account, the measurement equation is written in the matrix form:(16)S−W1W2W3⋯WKA1A2A3⋮AK=U,where S=s1s2s3⋮sN, Wk=Φk10⋯00Φk2⋯0⋮⋮⋱⋮00⋯ΦkN, Ak=ak1ak2⋮akN, U=μ1μ2μ3⋮μN, and N is the number of sampling points.


### 2.2.3. Global Solution

The system of equations (15) and (16) is an underdetermined system for the unknown amplitude envelope akt. The additional condition for the equation solution is that the variances of the nonhomogeneity terms εkn and noise μn have to be minimal while maintaining the given relationship between them.The weighted sum of equations (15) and (16) gives the loss function:(17)J=∑k∈ir2εTε+UTU=∑k∈ir2AkHGTGAk+ST−∑k∈iAkHWkHS−∑k∈iWkAk,where r is a weighting factor. The choice of a large value for the weighting factor r leads to the highly selective filtration in the frequency domain. In contrast, fast convergence with low-frequency resolution is achieved by choosing the small r value. The limit value of r [30] is shown in Table 1.Table 1
The value of the weighting factor.
Pole number
s
=
1
s
=
2
s
=
3
s
=
4
The max value ofr
7
×
10
6
4
×
10
6
2
×
10
6
1.1
×
10
6The first derivative of equation (17) with respect to the vector Ak gives a condition for the minimum of this function, which is called a normal equation:(18)∂J∂AkH=0.Thus the unknown envelopeAk results from equation (18).


## 2.2.1. Structural Equation

It is assumed that the amplitude envelopeakt is the low-frequency modulation of the complex carrier wave Φkt. In other words, envelope is locally approximated by a low-order polynomial, and it is filtered in some specific way. This condition can be expressed by a structural equation with the nonhomogeneity term εkn. The polynomial order designates the number of the filter poles. The equations for l-pole filter are given by(12)∇lakn=εkn.Takingl=2 as an example, the corresponding equation is given by(13)∇2akn=akn−1−2akn+akn+1=εkn.The system of the structural equations containing all the samples takes the following form, which can be rewritten in the matrix form:(14)−210⋯001−21⋯0001−2⋯00⋮⋮⋮⋱⋮⋮000⋯1−2ak1ak2ak3⋮akN=εk1εk2εk3⋮εkN,(15)GAk=ε.


## 2.2.2. Data Equation

Taking the measurement noise into account, the measurement equation is written in the matrix form:(16)S−W1W2W3⋯WKA1A2A3⋮AK=U,where S=s1s2s3⋮sN, Wk=Φk10⋯00Φk2⋯0⋮⋮⋱⋮00⋯ΦkN, Ak=ak1ak2⋮akN, U=μ1μ2μ3⋮μN, and N is the number of sampling points.


## 2.2.3. Global Solution

The system of equations (15) and (16) is an underdetermined system for the unknown amplitude envelope akt. The additional condition for the equation solution is that the variances of the nonhomogeneity terms εkn and noise μn have to be minimal while maintaining the given relationship between them.The weighted sum of equations (15) and (16) gives the loss function:(17)J=∑k∈ir2εTε+UTU=∑k∈ir2AkHGTGAk+ST−∑k∈iAkHWkHS−∑k∈iWkAk,where r is a weighting factor. The choice of a large value for the weighting factor r leads to the highly selective filtration in the frequency domain. In contrast, fast convergence with low-frequency resolution is achieved by choosing the small r value. The limit value of r [30] is shown in Table 1.Table 1
The value of the weighting factor.
Pole number
s
=
1
s
=
2
s
=
3
s
=
4
The max value ofr
7
×
10
6
4
×
10
6
2
×
10
6
1.1
×
10
6The first derivative of equation (17) with respect to the vector Ak gives a condition for the minimum of this function, which is called a normal equation:(18)∂J∂AkH=0.Thus the unknown envelopeAk results from equation (18).


## 3. Instantaneous Modal Parameter Identification of Linear Time-Varying Systems

### 3.1. Modal Response of Time-Varying Structures

The motion equation of aP-DOF LTV structure is given as(19)Mtx¨t+Ctx˙t+Ktxt=ft,where Mt, Ct, and Kt are time-varying mass, damping, and stiffness matrices, respectively; xt is the displacement response vector, x˙t and x¨t are the velocity and acceleration response vectors, respectively; and ft is the excitation force vector on the structure.Structural acceleration responses can be obtained as the superposition ofP vibration modal responses [31], and when an impulse force is applied on the p-th dof, the acceleration response of the k-th generalized modal coordinate is given as(20)x¨pt=∑k=1PA0ke−ξktω0kttcosωdktt+φ0k=∑k=1Px¨pkt=∑k=1PAktcosφkt,where ωdkt=ω0kt1−ξk2t and ξkt, ω0kt, and ωdkt denote the k-th modal damping ratio, modal frequency, and natural frequency, respectively. It can be further written as the superposition of multicomponent AM-FM signals.The acceleration impulse response signal of the time-varying structure is similar to multicomponent harmonic signals of equation (10). It can be adaptively decomposed by using the TF decomposition method proposed in this paper, thereby reconstructing each component signal, i.e., Aktcosφkt.Then, the instantaneous modal parameters can be identified by using the following estimator:(21)ωdkt=dφktdt,ηkt=dlnaktdt⋅1ωdkt,ξkt=ηkt2ηkt2+1.


### 3.2. Energy Analysis Method for the Identification of Damping Ratio

As already noted, the higher precision of the identification of the IF can be obtained from the well-established TF representation, but for damping ratio identification based on amplitude envelope in equation (21), it has much lower identification accuracy since the amplitude envelope is easily disturbed by noise. In order to enhance the identification accuracy of the damping ratio, the estimator based on signal energy analysis for LTV systems is presented.It can be assumed that the signal is time-invariant in a short time (that is, the frequency of the signal is not changed in a sampling period). For most signals without fast time-varying IF, this assumption can really be valid. Consider a time-varying signal as shown in Figure2; the zero position of this signal is denoted as tq, and the IF is considered to be constant within a sampling period Tq around tq. In other words, the IF of the signal from time tq−Tq/2 to time tq+Tq/2 is denoted as ωdktq. The signal energy of the first and second half sampling periods is calculated by Equ and Eqd, respectively. They are given by(22)Equ=∫tq−Tq/2tqx¨pk2tdt=A0k2∫tq−Tq/2tqe−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt,(23)Eqd=∫tqtq+Tq/2x¨pk2tdt=A0k2∫tqtq+Tq/2e−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt.Figure 2
A time-varying signal for energy analysis.According tocos2ωdktqt+φ0k=1−sin2ωdktqt+φ0k, equation (22) can be rewritten as(24)Equ=A0k22∫tq−Tq/2tqe−2ξktqω0ktqtdt+A0k22∫tq−Tq/2tqe−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt=A0k22×B1+A0k22×B2,where B1=∫tq−Tq/2tqe−2ξktqω0ktqtdt and B2=∫tq−Tq/2tqe−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt, which can be calculated by the real part of Γ2:(25)Γ2=∫tq−Tq/2tqe−2ξktqω0ktqt⋅ej2ωdktqt+φ0kdt=ej2φ0kj2ωdktq−2ξktqω0ktq⋅ej2ωdktq−2ξktqω0ktq⋅tq−ej2ωdktq−2ξktqω0ktq⋅tq−Tq/2=ej2φ0k⋅−j2ωdktq−2ξktqω0ktq4ωdk2tq+4ξk2tqω0k2tq⋅ej2ωdktq−2ξktqω0ktq⋅tq−ej2ωdktq−2ξktqω0ktq⋅tq−Tq/2.According to the inherent characteristics of the vibration system, the following conditions can be satisfied:(26)ωdktqTq=2π,ωdktq=ω0ktq1−ξk2tq.Then,Γ2 can be rewritten as(27)Γ2=ej2φ0k⋅−j2ωdktq−2ξktqω0ktq4ω0k2tq⋅ej2ωdktq−2ξktqω0ktq⋅tq−Tq/2⋅ej2ωdktq−2ξktqω0ktqTq/2−1=−j2ωdktq−2ξktqω0ktq4ω0k2tq⋅ej2ωdktqtq−Tq/2+φ0k⋅e−2ξktqω0ktq⋅tq−Tq/2⋅e−2πξktq/1−ξk2tq−1.In the same way, equation (23) can be rewritten as(28)Eqd=A0k22∫tqtq+Tq/2e−2ξktqω0ktqtdt+A0k22∫tqtq+Tq/2e−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt=A0k22×B3+A0k22×B4,where B4 can be calculated by the real part of Γ4:(29)Γ4=−j2ωdktq−2ξktqω0ktq4ω0k2tq⋅ej2ωdktqtq−Tq/2+φ0k⋅e−2ξktqω0ktq⋅tq⋅e−2πξktq/1−ξk2tq−1.From the above equations, the relationships can be obtained as follows:(30)B3=e−ξktqω0ktqTqB1,B4=e−ξktqω0ktqTqB2.Therefore, when the IF has already been identified, the instantaneous damping ratio can be identified by the relationship between signal energyEqd and signal energy Eqd:(31)Eqd=e−ξktqω0ktqTqEqu,ξktq=−lnEqdEqu/ω0ktqTq.In such a way, the damping ratioξktq at each point tq can be identified using the same method. Compared with the identification method based on amplitude envelope, the method based on energy analysis can effectively improve the robustness to noise.


## 3.1. Modal Response of Time-Varying Structures

The motion equation of aP-DOF LTV structure is given as(19)Mtx¨t+Ctx˙t+Ktxt=ft,where Mt, Ct, and Kt are time-varying mass, damping, and stiffness matrices, respectively; xt is the displacement response vector, x˙t and x¨t are the velocity and acceleration response vectors, respectively; and ft is the excitation force vector on the structure.Structural acceleration responses can be obtained as the superposition ofP vibration modal responses [31], and when an impulse force is applied on the p-th dof, the acceleration response of the k-th generalized modal coordinate is given as(20)x¨pt=∑k=1PA0ke−ξktω0kttcosωdktt+φ0k=∑k=1Px¨pkt=∑k=1PAktcosφkt,where ωdkt=ω0kt1−ξk2t and ξkt, ω0kt, and ωdkt denote the k-th modal damping ratio, modal frequency, and natural frequency, respectively. It can be further written as the superposition of multicomponent AM-FM signals.The acceleration impulse response signal of the time-varying structure is similar to multicomponent harmonic signals of equation (10). It can be adaptively decomposed by using the TF decomposition method proposed in this paper, thereby reconstructing each component signal, i.e., Aktcosφkt.Then, the instantaneous modal parameters can be identified by using the following estimator:(21)ωdkt=dφktdt,ηkt=dlnaktdt⋅1ωdkt,ξkt=ηkt2ηkt2+1.


## 3.2. Energy Analysis Method for the Identification of Damping Ratio

As already noted, the higher precision of the identification of the IF can be obtained from the well-established TF representation, but for damping ratio identification based on amplitude envelope in equation (21), it has much lower identification accuracy since the amplitude envelope is easily disturbed by noise. In order to enhance the identification accuracy of the damping ratio, the estimator based on signal energy analysis for LTV systems is presented.It can be assumed that the signal is time-invariant in a short time (that is, the frequency of the signal is not changed in a sampling period). For most signals without fast time-varying IF, this assumption can really be valid. Consider a time-varying signal as shown in Figure2; the zero position of this signal is denoted as tq, and the IF is considered to be constant within a sampling period Tq around tq. In other words, the IF of the signal from time tq−Tq/2 to time tq+Tq/2 is denoted as ωdktq. The signal energy of the first and second half sampling periods is calculated by Equ and Eqd, respectively. They are given by(22)Equ=∫tq−Tq/2tqx¨pk2tdt=A0k2∫tq−Tq/2tqe−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt,(23)Eqd=∫tqtq+Tq/2x¨pk2tdt=A0k2∫tqtq+Tq/2e−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt.Figure 2
A time-varying signal for energy analysis.According tocos2ωdktqt+φ0k=1−sin2ωdktqt+φ0k, equation (22) can be rewritten as(24)Equ=A0k22∫tq−Tq/2tqe−2ξktqω0ktqtdt+A0k22∫tq−Tq/2tqe−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt=A0k22×B1+A0k22×B2,where B1=∫tq−Tq/2tqe−2ξktqω0ktqtdt and B2=∫tq−Tq/2tqe−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt, which can be calculated by the real part of Γ2:(25)Γ2=∫tq−Tq/2tqe−2ξktqω0ktqt⋅ej2ωdktqt+φ0kdt=ej2φ0kj2ωdktq−2ξktqω0ktq⋅ej2ωdktq−2ξktqω0ktq⋅tq−ej2ωdktq−2ξktqω0ktq⋅tq−Tq/2=ej2φ0k⋅−j2ωdktq−2ξktqω0ktq4ωdk2tq+4ξk2tqω0k2tq⋅ej2ωdktq−2ξktqω0ktq⋅tq−ej2ωdktq−2ξktqω0ktq⋅tq−Tq/2.According to the inherent characteristics of the vibration system, the following conditions can be satisfied:(26)ωdktqTq=2π,ωdktq=ω0ktq1−ξk2tq.Then,Γ2 can be rewritten as(27)Γ2=ej2φ0k⋅−j2ωdktq−2ξktqω0ktq4ω0k2tq⋅ej2ωdktq−2ξktqω0ktq⋅tq−Tq/2⋅ej2ωdktq−2ξktqω0ktqTq/2−1=−j2ωdktq−2ξktqω0ktq4ω0k2tq⋅ej2ωdktqtq−Tq/2+φ0k⋅e−2ξktqω0ktq⋅tq−Tq/2⋅e−2πξktq/1−ξk2tq−1.In the same way, equation (23) can be rewritten as(28)Eqd=A0k22∫tqtq+Tq/2e−2ξktqω0ktqtdt+A0k22∫tqtq+Tq/2e−2ξktqω0ktqt⋅cos2ωdktqt+φ0kdt=A0k22×B3+A0k22×B4,where B4 can be calculated by the real part of Γ4:(29)Γ4=−j2ωdktq−2ξktqω0ktq4ω0k2tq⋅ej2ωdktqtq−Tq/2+φ0k⋅e−2ξktqω0ktq⋅tq⋅e−2πξktq/1−ξk2tq−1.From the above equations, the relationships can be obtained as follows:(30)B3=e−ξktqω0ktqTqB1,B4=e−ξktqω0ktqTqB2.Therefore, when the IF has already been identified, the instantaneous damping ratio can be identified by the relationship between signal energyEqd and signal energy Eqd:(31)Eqd=e−ξktqω0ktqTqEqu,ξktq=−lnEqdEqu/ω0ktqTq.In such a way, the damping ratioξktq at each point tq can be identified using the same method. Compared with the identification method based on amplitude envelope, the method based on energy analysis can effectively improve the robustness to noise.


## 4. Simulation Example

In order to validate the proposed method, a 3-DOF mass-spring-dashpot system with time-varying stiffness and time-varying damping is simulated as a numerical example. The LTV system is shown in Figure3.Figure 3
3-DOF time-varying system.The initial parameters are given by(32)m1=0.5,m2=m3=1kg,k1=8000,k2=10000,k3=6000,k4=4000N/m,c1=c2=0.2,c3=c4=0.3Ns/m.The stiffness and damping variations for the case are given by(33)Δk1=12000atan10×−22−atan10t−22,Δk3=1200sin1.5πt⋅t,Δk2=Δk4=0,Δc1=Δc2=0,Δc3=−0.3sin0.8πt,Δc4=0.3t.Assuming that an instantaneous impulse excitation is applied on the massm2 at the initial time, the Newmarkβ algorithm is used to compute the structural free response. The sampling time is 3 s, and the acceleration response data of the mass m2 are shown in Figure 4. Then the identification method can be performed on the data.Figure 4
Acceleration response signal.To further validate the adaptiveness and improvement for the TF representation based on GLCT, the Morlet wavelet is selected as a TFA tool to extract the IF as a comparison. The TFA results are shown in Figures5 and 6, respectively. It can be seen that the response signal contains three components. The IF of the first component varies slowly with time, and the IF of the other two components varies rapidly. For the TF representation around point 2s,20Hz, instantaneous frequency trajectory based on GLCT is more clear and obvious and the TF representation based on GLCT has a higher resolution. However, at TF intervals 1−2s,15−25Hz and 1−3s,35−45Hz, the Morlet wavelet appears to be inadequate for TF representation of these signal components with fast time-varying IF.Figure 5
TFA results of GLCT.Figure 6
TFA results of the Morlet wavelet.Therefore, the IF of the time-varying system can be identified by TF ridge extraction. The IF of each mode extracted by using the Morlet and GLCT are shown in Figures7–9. Black lines in these figures are the theoretical values, and red dotted lines and blue dotted lines are the identified results by using GLCT and Morlet, respectively. The results clearly indicate that GLCT conduces to better performance in tracking fast time-varying IF than Morlet.Figure 7
Theoretical results and identified results of the 1st mode.Figure 8
Theoretical results and identified results of the 2nd mode.Figure 9
Theoretical results and identified results of the 3rd mode.Once the IF of each mode is identified, the three time-varying components of the acceleration response can be exactly separated by the VK filter. The decomposition results are shown in Figure10.Figure 10
Amplitude envelope and component signal.Moreover, the amplitude envelope of the three components is also computed, which can be used for the identification of the time-varying damping ratio. To further validate the robustness of the method based on energy analysis, the estimator based on signal amplitude in equation (21) is used to estimate the time-varying damping ratio as a comparison. The identified results are presented in Figures 11–13, and significant fluctuations are observed by using the estimator in equation (21). By comparison, the energy method based on half-period integral can effectively resist noise and get higher identification accuracy. Since the identification of the damping ratio is full of challenge and most modal identification methods are invalid for LTV systems, important improvements have been made in instantaneous modal parameter identification by using energy analysis. In order to investigate the performance of the identification method, a root-mean-square error (RMSE) over the total sampling time is further computed to quantify the accuracy of identified results as(34)RMSE=1Ns∑i=1Nsϑi−ϑ˜i2,where Ns represents the number of sampling points and ϑi and ϑ˜i represent the theoretical and identified instantaneous modal parameters, i.e., instantaneous frequencies and damping ratios, respectively.Therefore, the calculated RMSE values of instantaneous frequencies and damping ratios are listed in Table 2.Figure 11
Theoretical results and identified results of the 1st mode.Figure 12
Theoretical results and identified results of the 2nd mode.Figure 13
Theoretical results and identified results of the 3rd mode.Table 2
Errors of identification results of instantaneous frequency.
RMSE
Mode
Morlet
GLCT
Instantaneous frequency
1st
0.306
0.195
2nd
1.882
0.230
3rd
0.599
0.331
RMSE
Mode
Time-domain algorithm
Energy algorithm
Instantaneous damping ratio
1st
3.006
×
10
−
3
3.689
×
10
−
4
2nd
3.944
×
10
−
3
8.792
×
10
−
4
3rd
2.511
×
10
−
3
6.699
×
10
−
4To further investigate the performance and reliability of the proposed method, different signal-to-noise ratios (SNR) of noise are added to the response of this structure. The same procedure as above is followed to identify the instantaneous modal parameters, and a mean absolute percentage error (MAPE) over the total sampling time is further computed to quantify the accuracy of the identified results as(35)MAPE=1Ns∑i=1Nsϑi−ϑ˜iϑi×100%,where Ns represents the number of sampling points and ϑi and ϑ˜i represent the theoretical and identified instantaneous modal parameters.Therefore, the calculatedMAPE values under the effects of measurement noise are listed in Tables 3 and 4.Table 3
Errors of identification results of instantaneous frequency.
MAPE (%)
Mode
SNR (dB)
Noise-free
100
50
20
Morlet
1st
2.421
2.444
2.444
2.514
2nd
8.141
8.152
8.153
8.157
3rd
1.091
1.076
1.081
1.136
GLCT
1st
1.553
1.556
1.557
1.571
2nd
1.305
1.312
1.314
1.316
3rd
0.911
0.905
0.907
0.927Table 4
Errors of identification results of the instantaneous damping ratio.
MAPE (%)
Mode
SNR (dB)
Noise-free
100
50
20
Time-domain algorithm
1st
>50
>80
>80
>100
2nd
>50
>80
>80
>100
3rd
>50
>80
>80
>100
Energy algorithm
1st
14.135
16.929
16.429
21.008
2nd
20.873
21.560
21.607
21.396
3rd
26.783
29.496
29.531
31.387By comparing theRMSE values and MAPE values using different identification methods, it can be found that the proposed method has better performance in tracking time-varying parameters. The results also show the good adaptive ability and high accuracy even in low SNR scenarios. Especially for the identification of damping ratio, the energy analysis algorithm outperforms the traditional time-domain method.For practical application of the proposed method, instantaneous frequencies and damping ratios can be identified only using one set of data from single-output. However, for the estimation of mode shapes of LTV systems, there are two hard issues in the implementation of this algorithm. Firstly, multiple sets of data from multioutput are needed to establish the mode shape matrix. Secondly, it is necessary to identify the time-varying mode shapes at each time point, which will lead to a large amount of calculation. Therefore, further studies should be done to estimate mode shapes and improve algorithm efficiency.


## 5. Conclusion

This paper proposes a method based on adaptive TF decomposition for instantaneous modal parameter identification of LTV systems. This method consists of two main steps. The first step is the IF extraction algorithm based on GLCT. By applying the algorithm to perform the VK filter, the free response of LTV systems can be decomposed into multiple mode responses. Then the other step is the energy analysis algorithm based on each mode response for the identification of the instantaneous damping ratio. Based on the results of numerical simulations, the corresponding conclusions can be concluded as following:(1)
GLCT is a more effective TFA tool than WT for the response signals of LTV systems. It can obtain the better TF representation with high energy concentration and therefore extract the IF with better accuracy.(2)
Since the wavelet transform is essentially a signal filter, application of a TFA tool with high energy concentration can enhance antinoise capability of parameter identification. By applying GLCT to perform the VK filter, it not only improves the adaptability of decomposition but also has good antinoise ability. It can be seen from the results of error analysis that our method has good immunity to noise.(3)
By adding different SNRs of noise, identification of the instantaneous damping ratio based on the time-domain algorithm and energy analysis algorithm is also discussed. The energy analysis algorithm can eliminate the influence of amplitude fluctuation in the time domain and improve the identification accuracy of the damping ratio. The results show that the energy analysis algorithm has better recognition accuracy, convergence, and noise immunity compared with the traditional method.



---

*Source: 1475981-2019-12-28.xml*