Numerical and Experimental Investigation on Parameter Identification of Time-Varying Dynamical System Using Hilbert Transform and Empirical Mode Decomposition

This paper proposes an approach to identifying time-varying structural modal parameters using the Hilbert transform and empiricalmode decomposition.Definition of instantaneous frequency and instantaneous damping ratio based onHilbert transform for single-degree-of-freedom (SDOF) system is first introduced.The following is the definition of Hilbert damping spectrum from which the time-varying damping ratio of multi-degree-of-freedom (MDOF) system can be calculated. Identification procedures for both instantaneous frequency and damping ratios based on their definition are then introduced. Applicability of the proposed identification algorithm has been validated through several numerical examples. The instantaneous frequency and damping ratios of SDOF system under free vibration and under sinusoidal and white noise excitation have been identified.The proposedmethod is also applied toMDOF system with slow and sudden changing structural parameters.The results demonstrate that when the system modal parameters are slowly changing, the instantaneous frequency could be easily andwell identifiedwith satisfied accuracy for all cases. However, the instantaneous damping ratio could be extracted only when the system is lightly damped. The damping results are better for free vibration situation than for the forced vibration cases. It is also shown that the suggested method can easily track the abrupt change of system modal parameter under free vibration. The proposed method is then applied to a 12-story short-lag shear wall structure model tested on a shaking table. The instantaneous dynamic properties of the structure were identified and were then introduced as known parameters into a finite element model. Comparisons with the numerical results using constant structural parameters demonstrate that the calculated structural responses using the identified time-varying parameters are much closer to the experimental results.


Introduction
Quite some structures exhibit intrinsically time-varying behavior when subjected to strong excitations as earthquake or wind or when they are damaged.However, most existing structural identification techniques rest on the assumption that the structures are linear or time invariant.They are not readily applicable to examine the possible time-varying behavior during online monitoring of some civil structures.It is therefore necessary to develop identification techniques that can quantify the time-varying structural behavior for the purpose of damage detection, health monitoring, reliability, and safety evaluation of structures.
A mount of research has been carried out regarding the online estimation (tracking) technique of systems with timevarying parameters.Smyth et al. [1] had obtained online estimation of the parameters of multi-degree-of-freedom (MDOF) nonlinear hysteretic systems based on the measurement of restoring forces.Zhang et al. [2] identified dynamic properties of linear single-degree-of-freedom (SDOF) system with gradual changing stiffness based on empirical mode decomposition (EMD) plus Hilbert transform (HT) method and pointed out that the steady part and transient part of structural dynamic responses could be separated through Hilbert spectrum analysis when the system is under sine wave excitation.Pines and Salvino [3] proposed Hilbert damping 2 Mathematical Problems in Engineering spectrum based on intrinsic mode functions decomposed by the EMD method.The damping loss factor was first proposed as a joint distribution function of time and frequency; afterwards Hilbert marginal spectrum was developed and applied to identify the damage of a 3-story shearing building in a shaking table test.Shi and Law [4] addressed the identification of linear time-varying multi-degrees-of-freedom systems by HT+EMD.They used several numerical examples to demonstrate the effectiveness and accuracy of the proposed method.Basu et al. [5] developed online identification of linear time-varying stiffness of structural systems based on wavelet analysis.L-P wavelet was adopted in their study to identify MDOF system with slowly varying stiffness and sudden changing stiffness.Pai et al. [6] compared EMD plus Hilbert transform method with sliding window in timefrequency analysis of nonlinear signal.They pointed out that EMD plus Hilbert transform method as a self-adaptive method needs no basis function and has high frequency resolution dealing with nonstationary and nonlinear signal.However, they also pointed out that Gibbs phenomenon could lead to inaccuracy at the ends.Wang and Genda [7] suggested a recursive Hilbert transform method for the timevarying property identification of shear-type buildings under base excitation.With known floor masses, the stiffness and damping coefficients of each floor were identified one by one from the top to bottom.It is worth noting that the Hilbert transform plus EMD method has also been widely used in parameter identification of linear structures [8,9].The author also has applied the HT+EMD approach for linear identification problems as parameters identification of structures with closed space modal properties [10], structural damage detection [11], and parameter identification of an existing longspan bridge [12].We further demonstrated in Chen et al. [13] that EMD is a very powerful tool in processing nonstationary signal which is frequently encountered in time-varying identification problem.Though successful in numerical examples, experimental verifications of the abovementioned techniques, especially for time-varying system, are still rare.
Inspired by the previous researches regarding the timevarying structural parameter identification problem, this paper presents an approach for identifying time-varying structural parameters using HT+EMD.Definitions of instantaneous frequency and damping ratios based on Hilbert transform in SDOF with time-varying parameters are first introduced.For SDOF, the time-varying parameters can be easily obtained by its definition.For MDOF system, the structural responses will be first processed by EMD to obtain the modal response based on which the time-varying parameters can be identified using the same procedure for SDOF system.
Extensive numerical examples and experimental measurements of a large-scale 12-story building model, which was tested on a shaking table, have been employed to demonstrate the applicability and effectiveness of the proposed method.

Definitions of Instantaneous Parameters by Hilbert Transform
Traditionally frequency is defined as the number of occurrences of a repeating event per unit time.In the classical Fourier transform analysis, the repeating event is regarded as sine or cosine wave that has constant amplitude.According to this definition, any sine or cosine wave whose duration is less than one period cannot lead to a meaningful frequency value.However, as we know, lots of random processes in nature are nonstationary such as earthquake wave or wind speed, it is difficult for us to capture the local modulation properties of nonstationary signals.
In order to depict local properties of nonstationary signals, Leon Cohen proposed the instantaneous frequency systematically and defined it as the first derivative of phase angle.For an arbitrary time series, (), we can always have its Hilbert transform, (), as where ..indicates the Cauchy principal value integral.
According to this definition, () and () are actually the complex conjugate pair, so we can have an analytic signal, (), as in which Equation ( 1) defines the Hilbert transform as the convolution integral of () with 1/, and it therefore emphasizes the local properties of ().In (2) the polar coordinate expression further clarifies the local nature of this representation: it is the best local fit of amplitude and phase varying trigonometric function to ().
With the above Hilbert transform in (1), the instantaneous frequency   is defined as Equation ( 4) can be further deduced as follows: In order to clarify the meaning of (4) physically, Cohen introduced the term "monocomponent function" as a limitation on the original signal (), which means that only when () is a narrow band signal could the instantaneous frequency reflect the local vibration property.
Based on the Hilbert transform and definition of instantaneous frequency, instantaneous damping ratio of SDOF system could also be deduced from the equation of motion.
Consider the equation of motion of a SDOF system with time-varying parameters as follows:  () ü () +  () u () +  ()  () = 0. (6) After normalization by the mass property, (6) could be rewritten as where  0 () and () indicate the time-dependent frequency and damping ratio; through Hilbert transform we know that the analytical signal of () could be formed as where () and () are regarded as the instantaneous amplitude and instantaneous phase angle, respectively.The first and second derivatives of analytical signal () could be easily got as follows: Apply HT to (7) leading to Introducing ( 9) into (10) we get Let the real part and imaginary part in (11) be equal to zero, respectively; we can finally determine the time-dependent damping ratio () as follows: where  and  are instantaneous amplitude and frequency, respectively.Therefore, the instantaneous damping ratio of SDOF system under zero excitation could be defined as a joint function of amplitude and frequency.

Identification Methodology
With definitions of the instantaneous frequency and damping ratio for SDOF system in the previous section, identification algorithms are proposed based on empirical mode decomposition method in this part.The empirical mode decomposition (EMD) method was proposed by Huang et al. in 1998 [14].EMD can decompose any measured signal into intrinsic mode function (IMF) that admits a well-behaved Hilbert transform.EMD is a self-adaptive data processing Time (s) method and has higher frequency resolution compared with traditional signal processing tools in dealing with nonstationary signals.
Firstly, in order to obtain the "monocomponent function" as required by HT, it is necessary to apply EMD to the original signal () to decompose it into summation of several IMF components plus the final residue.That is, where   () is intrinsic mode function, which is a narrow band signal and   () is residual.Therefore, the multicomponent signal () is represented by several monocomponent signals IMFs.Each IMF corresponds to one mode; in this way, dynamic responses of MDOF system could be decoupled into several intrinsic mode responses like   ().Details of the implementation of EMD can be found in Huang et al. [14] and Chen and Xu [10,12].
Having   () obtained, HT can be used to calculate the instantaneous phase angle as defined in (3a)-(3b).Then, according to (4) three difference algorithms are introduced here to estimate   , and the first is forward difference: The second is backward difference: The third is central difference: These three algorithms almost have the same accuracy, and they are all sensitive to noise.Therefore, necessary smoothing or fitting techniques could be taken on f () to eliminate numerical errors.It should be noted that the identification results are also sensitive to calculating parameters adopted in EMD algorithm.In general, a stricter decomposing standard should be taken when the input signal is highly contaminated.Take a nonlinear signal () as example; Figure 1 shows the time history of ():  Equation (17) indicates that () has two modal components, and the theoretical value of instantaneous frequency of each component is given in the following: After empirical mode decomposition we got seven IMFs as depicted in Figure 2. The amplitude of the first two IMFs is significantly larger than the last five IMFs, which are actually due to the numerical calculation error of EMD.By applying Hilbert transform to the first two IMFs, we can get the instantaneous phase angle from which the instantaneous frequency can be estimated by the forward difference algorithm given (14).Figure 3 shows the identified results compared with theoretical ones.It could be observed that identified results fit well with the theoretical ones.Identification of instantaneous damping ratios for SDOF system has been derived in (12); however, it could be used only in free vibration cases; when the system is forced vibration, EMD in conjunction with random decrement technique should be used first to extract the free responses.In order to extend this method to MDOF systems Pines and Salvino, 2006 [3] proposed Hilbert damping spectrum to evaluate damping; he defined damping lose factor () = 2() and demonstrated that  is a function of both frequency and time as follows: According to (19), each IMF has a time-dependent damping loss factor (), and as already mentioned each IMF corresponds to each single mode.In this sense, the instantaneous damping ratios can be estimated from (19) as a marginal distribution of time or frequency from (19).The timedependent damping ratio, for instance, will be In summary, the procedure of time-varying parameter identification of MDOF system consists of the following three steps.First, decompose the measured response of the system by EMD into summation monocomponent signal.Then, for each monocomponent signal, the instantaneous frequency can be estimated by ( 14), (15), or (16).The instantaneous damping ratio can be estimated by (20).

Numerical Simulations
In order to illustrate the proposed algorithm for identification of SDOF and MDOF systems with time-varying parameters, several numerical examples are discussed in this section including SDOF system with slow varying and fast varying parameters and MDOF system with slow varying and sudden change parameters.
The initial velocity and displacement of the system are 10 mm/s and 10 mm, respectively, and time interval is set as 0.001 sec.The free response time histories of displacement, velocity, and acceleration, as shown in Figure 4, were calculated by fourth-order Runge-Kutta method.Take the acceleration response as the measured signal (); the instantaneous phase angle () and amplitude () can be obtained by Hilbert transform.The instantaneous frequency and damping ratios can then be identified by (14) and (12).Second-order five-point smoothing was applied to the identified instantaneous damping ratios to reduce the numerical noise due to difference calculation.Figures 5(a) and 5(b) compare, respectively, the identified instantaneous frequency and damping ratio with theoretical value.
It is found that the identified instantaneous frequency agrees very well with the theoretical value, excluding the results at endpoints.The difference at endpoints is called end effect that is caused by Hilbert transform itself, since the signal is always cut limited while the convolution integral in theory is from −∞ to +∞.Therefore, inaccuracy will occur at endpoints according to numerical difference methods.Figure 5(b) shows that the identified damping ratios are fluctuating around the theoretical value.Table 1 further gives the relative error of identified and theoretical parameters at different time instants.Note that the average identification error of frequency and damping ratio is less than 0.5% and 3%, respectively.

SDOF: Excited Vibration.
The equation of motion of a SDOF system under external excitation can be expressed as follows: where () is system input and () and () are timedependent damping ratios and frequency, respectively.For the first case, we assume the external input as () = 2 −0.04 sin(20+0.1 2 ) and the initial velocity and displacement are zero.The time-varying parameters were set as 0.01 times that in Example 1 and the time step was 0.001 sec.The dynamics responses of the system are shown in Figure 7. Figure 8 gives the Hilbert energy spectrum of the acceleration response.It can be easily observed from Figure 8 that there are two time scales in the response; one is in the frequency range from 10 to 12 Hz and the other is from 3 to 6 Hz.The second one is in accordance with the natural frequency of the SDOF system.Therefore, we take the second IMF as the free vibration response of SDOF system, and following the above steps in the first example, instantaneous frequency and damping ratios could be easily identified.The results are compared with theoretical values as in Figures 9(a Table 2 shows relative errors at several calculation time instants.It is found in Table 2 that identified results for frequency are relatively stable in the middle part of the response with relative error around 1%. Significant errors occur at the end due to the end effect the same as Example 1.For damping ratios, it should be mentioned that the results are fitted by second-order 5 points smoothing method, although    the identified results are not ideal; with average relative error 5%, the trend of time-varying damping ratios could be clearly tracked by this method.White noise excitation was also considered in this example.The identification accuracies for frequency and damping were similar.

MDOF System with Slowly
Varying Parameter.To illustrate the application of the identifying methodology for MDOF systems, an example of a 2DOF system has been considered.The system considered is a shear-building model (as shown in Figure 10) with 2 mass lumped at the two nodes, mm), dynamic responses of system can be calculated by fourth-order Runge-Kutta method, Δ = 0.001 s, and the results of the first floor are shown in Figure 11.Apply EMD to the acceleration response of the first floor leading to three IMF components as shown in Figure 12.The first IMF corresponds to the modal response of the second vibration mode and the second IMF corresponds to the first vibration mode.Following the identification procedures of SDOF system, taking the two IMFs decomposed as free responses of SDOF system, the instantaneous frequency and damping ratios could be easily identified.Figures 13(a), 13(b),  and 13(c) compare, respectively, the identified time-varying frequencies, damping ration of the first, and second vibration mode with the theoretical values.Table 3 further shows the relative identification error at several time instants.Note that the frequency results fit very well with the theoretical values, while the damping ratio results are less ideal, especially at the curve peak or valley.The variation trend of the instantaneous damping can be tracked by the identification results.

MDOF System with Sudden-Change Parameter.
To check the effectiveness of the method for structures with suddenchange parameter, the same 2DOF system above is adopted   but assume that an abrupt change of the first floor's stiffness happens, which is Mathematical Problems in Engineering with the initial conditions {(0)} = {10 0 10 0}, (unit: mm); then responses of the structure are obtained by fourthorder Runge-Kutta method, time interval Δ = 0.001 s. Figure 14 shows the first floor's displacement, velocity, and acceleration response, from which three IMF components and one residual can be obtained by EMD (Figure 15). Figure 16 compares the identified instantaneous frequency with the theoretical values.Note that the change of frequency can be well identified before and after the sudden change of stiffness.17).The structural model was 3.51 m high, 2.81 m wide, and 3 m long as indicated in Figure 18(a).Schematic plan of standard floor of the model with wall numbering and measurement sensor arrangement is illustrated in Figures 18(a) and 18(b).The sensors adopted in this test were piezoelectric accelerometers and displacement  Using the acceleration records of wall Q7, the constant natural frequency and damping ratio can be identified by method previously suggested by the author of [12].The identification results are given in Table 4, in which the amplitude of input  earthquake increased from Case 1 to Case 3 (low, moderate, and high earthquake intensity).Note that the natural frequency has a sharp decrease from 5.127 Hz of the intact structure to 2.441 Hz after Case 3, indicating occurrence and accumulation of damage from test Case 1 to Case 3. It is interesting that from Case 1 to Case 3 the damping ratio increases from 3.5% to 11.06%, implying that the structure has increased energy dissipating capacity when damaged.

Experimental Investigation
For test Case 1 the instantaneous frequencies and damping ratios are identified by the proposed method using displacement responses at the top of walls Q7 and Q8, and the identified results are given in Figure 20.It is seen that, at most part of the time duration, except near the two ends of the record, the identified instantaneous frequency varies slightly around 5 Hz.Test Case 2, the instantaneous frequency, and damping ratio are identified using response of wall Q7, and the results are shown in Figure 21.It is seen that the instantaneous frequency drops quickly from 6 Hz to 2 Hz at the beginning and then varies around 2 Hz for the rest period.As for the damping ratio, it goes down from 8% at the very beginning and varies in a range from 2% to 6% and then in a range from 4% to 6%. Figure 22 compares the marginal damping-frequency spectrum (by integrating (19) over time) of response of Q7 at Case 1 (low level earthquake intensity) with that at Case 2 (high level earthquake intensity).It is observed that, for low level earthquake excitation, the damping ratio changes not much from 2 Hz to 5 Hz.For high level earthquake excitation the damping ratio is higher in lower frequency region as 1 Hz to 3 Hz than other regions, indicating that the input energy is mainly dissipated in vibration modes with lower frequency.The conclusion agrees well with the experimental observations that when the structure model was damaged more energy will be dissipated during opening and closing of existing cracks and development of new cracks.
Time-varying frequency and damping ratio of a structural model under various earthquake excitations are identified with the method proposed, in order to testify the identified results; 3D finite element simulation of the shaking table test is performed.First, a 3D FE model is established, with different inputs of earthquake waves, the displacements  at different shearing walls on top floor are calculated, then, with the identified real time frequency and damping ratios under different earthquake cases, the original FE model is modified, and the displacements are recalculated as   ; then  and   are compared, respectively, with real experimental values; it is easily found that   is closer in amplitude to experimental results, particularly in free vibration piece after main shock.

FE Simulation.
To further compare the effect of constant parameters with the time-varying parameters on structural response prediction, a three-dimensional finite element model of the building was established as shown in Figure 23.Beam 188 element, Shell63 element, and Shell43 element were employed to simulate beam, floor, and shearing wall, respectively.Concrete material properties were determined from control specimens cast during the construction procedure of the model.Floor density was adjusted to account for the additional mass added to the building model in the test.Model analysis on the FE model was performed.The calculated natural frequencies of the first three vibration modes are 5.075 Hz, 6.204 Hz, and 6.785 Hz in the  translational direction, torsional direction, and  translational direction, respectively.The relative error of the first frequency between FE model (5.075Hz) and real model (5.127Hz) is only 1.2% which means the mass/stiffness ratio of the FE model is very close to the real one.
Taking the actual record of earthquake wave in test Case 1 as input to the FE model, the structural responses can be calculated for two situations: first, a constant damping ratio 5% is adopted, and second, the identified time-varying parameters (Figure 20 For the second situation the identified instantaneous parameters are introduced into the above two equations to calculate the instantaneous damping coefficient. Figures 24 and 25 compare, for the constant damping situation and time-varying damping, respectively, the calculated displacements of walls Q7, Q8, and Q10 with the measured responses.It is found that when constant damping is adopted the response amplitude calculated by FE model is much smaller than the measured value, which indicates that the constant damping ratio of 5% is overestimated.The significant difference occurs in duration of 2-4 second where the calculated amplitude is only about half of the measured ones.When time-varying damping is considered, on the other hand, the calculated responses are much closer to the measured responses especially for the duration 6-10 seconds when the structural responses could be broadly treated as free vibration response.The root-mean-square value (RMS) of the calculated responses and the measured responses are computed and given in Table 5.The relative prediction error for time-varying damping situation is much smaller compared to that of the constant damping situation.For response of wall Q7 and Q10 the relative error is no more than 10%.This comparison shows that time-dependent parameters by the proposed method are reasonable.

Figure 8 :
Figure 8: Hilbert energy spectrum for signal in Figure 7.

2 )Figure 11 :
Figure 11: Dynamic responses of the first floor.

Figure 17 :
Figure 17: Panoramic view of the building model.

Figure 19 :
Figure 19: Measured acceleration and displacement responses of wall Q7 Q8 and Q10.

5. 2 .
Parameter Identification.White noise sweep frequency tests were carried out in the test after each loading case.

Figure 21 :
Figure 21: The displacement record of Q7 (a) and identified instantaneous frequency (b) and damping ratio (c).

Figure 22 :
Figure 22: Comparison of Hilbert marginal spectrum of high and low level earthquake excitation.

Figure 23 :
Figure 23: 3D finite element model of the test structure.
(a)) are used.Since Rayleigh damping model is used in ANSYS, for the first situation the two damping coefficients alpha and beta are determined as follows;

Table 1 :
Identification error of frequency and damping ratio at different time instant (SDOF Example 1).

Table 2 :
Identification error of frequency and damping ratio at different time instant (SDOF Example 2).

Table 3 :
Identification error of frequency and damping ratio at different time instant (MDOF Example 1).

Table 4 :
Identified natural frequency damping ratio of the 12-story building model.

Table 5 :
Root mean square value of measured and calculated responses.