Phase space interrogation of the empirical response modes for seismically excited structures

Mechanical Systems and Signal Processing 91 (2017) 250–265

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Phase space interrogation of the empirical response modes for seismically excited structures Bibhas Paul, Riya C. George, Sudib K. Mishra ⇑ Department of Civil Engineering, Indian Institute of Technology, Kanpur, UP 208016, India

a r t i c l e

i n f o

Article history: Received 3 May 2015 Received in revised form 30 November 2016 Accepted 5 December 2016

Keywords: Phase portrait Embedding dimension Topology Seismic Empirical mode decomposition

a b s t r a c t Conventional Phase Space Interrogation (PSI) for structural damage assessment relies on exciting the structure with low dimensional chaotic waveform, thereby, significantly limiting their applicability to large structures. The PSI technique is presently extended for structure subjected to seismic excitations. The high dimensionality of the phase space for seismic response(s) are overcome by the Empirical Mode Decomposition (EMD), decomposing the responses to a number of intrinsic low dimensional oscillatory modes, referred as Intrinsic Mode Functions (IMFs). Along with their low dimensionality, a few IMFs, retain sufficient information of the system dynamics to reflect the damage induced changes. The mutually conflicting nature of low-dimensionality and the sufficiency of dynamic information are taken care by the optimal choice of the IMF(s), which is shown to be the third/fourth IMFs. The optimal IMF(s) are employed for the reconstruction of the Phase space attractor following Taken’s embedding theorem. The widely referred Changes in Phase Space Topology (CPST) feature is then employed on these Phase portrait(s) to derive the damage sensitive feature, referred as the CPST of the IMFs (CPSTIMF). The legitimacy of the CPST-IMF is established as a damage sensitive feature by assessing its variation with a number of damage scenarios benchmarked in the IASC-ASCE building. The damage localization capability, remarkable tolerance to noise contamination and the robustness under different seismic excitations of the feature are demonstrated. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Interrogation of structures are essential to acquire information about the ‘‘health” of a structure, in view of incipient damages that may occur during their service lifetime. This process is referred as Structural Health Monitoring (SHM). The SHM assess the state of structures (stiffness and/or strength) due to changes in the material properties, geometric dimensions or boundary conditions. Earlier technologies for the SHM are visual inspection and nondestructive testing [1]. However, these techniques require full access to the structure that might not be feasible. Manual inspections are also time and resource consuming. Further, all type of damages cannot be detected by eyes. In fact, hidden damage(s) in structures resulted in a number of failures in the past that have triggered the development of the modern SHM methodologies. Modern SHM is automatic, continuous, robust and reliable in real time performance. The SHM not only provides information about the structural health under aging and degradation from operational environment but also helps in rapid assessment and screening after an extreme event (e.g. earthquake) for the safe occupancy and reuse of a structure. ⇑ Corresponding author. E-mail address: [email protected] (S.K. Mishra). http://dx.doi.org/10.1016/j.ymssp.2016.12.008 0888-3270/Ó 2017 Elsevier Ltd. All rights reserved.

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

251

Modern SHM involves sensor network for continuous monitoring of structural responses, which are subsequently processed for extracting information on structural health [2,3]. Statistical analysis are also taken up to assess the robustness of the feature with respect to severity and localization of damage. The SHM works at several level, Level-I refers to detection of damage in the structure, level II identifies the location of the damage, Level-III would identify the type of damage and the quantification of damage is done in Level-IV. Finally, Level (V) predicts on the remaining lifetime and prognosis [3]. Although, the SHM involves several components, the principal challenge to the structural engineers is the analysis of the sensed response(s) in order to extract effective information on structural damage. The SHM methodologies can be broadly classified into two categories, one is the System Identification (SI) based methods, and the other is based on extraction of damage features indicating the presence of damage(s). The SI methodology directly estimates the structural parameters from the responses measured at several locations of the structure by using an optimization algorithm. Although most of the SI methodologies deal with the linear systems, identification of nonlinear systems have also been tried with limited success [4–8]. However, the system identification approach becomes prohibitively exhaustive for large structures with significant amount of nonlinearities and the problems may suffer from non-uniqueness [9,10]. The feature based techniques extract damage sensitive features from the responses collected at several locations of the structure and relate them to the extent and location of damage(s). Most common damage sensitive features are the change in natural frequencies [11], change in mode shape curvature [12], Auto Regressive and Moving Average (ARMA) coefficients in the predictive time series model [13], Wavelet based techniques [14] and so on. A review of these vibration based methodology can be obtained from Doebling et al. [15] and Sohn et al. [16]. Recently, several damage features based on the phase space representation of the system dynamics are proposed. At motion, the dynamical states of a structure trace out certain geometry in phase space, referred as Phase portrait. It has been demonstrated that the damage(s) in the structures can be correlated with certain properties of the phase portrait [17]. The superiority of these features are established over the conventional modal features [18]. The phase portrait based techniques are also readily extendible to nonlinear structures due to incipient damage. The aspect of nonlinearity are often overruled in the conventional SHM methodologies, partly because of their inability to accommodate nonlinearity. Commonly reported damage features based on the disparities of the phase portrait(s) are, the Local Attractor Variance Ratio (LAVR) proposed by Todd et al. [19], Phase Space Warping (PSW) suggested by Chelidze and Liu [20], change of Phase Space Topology (CPST) by Nie et al. [21,22]. A basic requirement of these techniques is that the structure must be excited by a low dimensional chaotic signal to ensure that the low dimensionality of the phase portrait is maintained. Similar methodology under stochastic excitations have also been reported [23,24]. Improvement in the chaotic interrogation based feature are made by employing the hyper-chaotic excitation by Torkamani et al. [25]. Although phase portrait based damage detection methodologies based on linear/nonlinear dynamic responses received considerable attention, its applicability to large scale structural system is limited by the fact that, exciting a structure by chaotic waveform poses practical difficulty. The suitability of the PSI can be enhanced largely, if instead of chaotic waveform, the natural excitations can be made use of. With this being the eventual goal, in this study, an attempt is made to use the measured responses of the structures, subjected to seismic excitations, for use in the PSI. In this regard, it may be pointed out that similar attempt have been made by Nichols [26] using the hydrodynamic wave excitations for the SHM of offshore structure. However, the wave excitations and the structural responses are low dimensional. The present work extends the chaotic interrogation procedure for the SHM of structures subjected to earthquake excitations, which are known to be high dimensional and difficult to embed in phase space using the Taken’s theorem. The IASC-ASCE benchmark building is selected as the structure to obtain its response(s) under recorded ground accelerations. These responses are taken as the measured responses to assess the incipient damage(s), localization and severity. Five damage scenarios, benchmarked in the IASC-ASCE building are adopted for numerical elucidation. The difficulty associated with the high dimensionality of the seismic responses are overcome with the help of EMD, decomposing a high dimensional signal with a few relatively low dimensional components, referred as Intrinsic mode Functions (IMFs). The IMFs retain significant portion of the dynamics to reflect their sensitivity towards the induced damage scenarios. The CPST features pertaining to individual IMFs (referred as IMF-CPST) are then assessed for their suitability as the damage features. The robustness of this feature is assessed under varying degree of noise contamination as well as alternative set of seismic excitations.

2. Phase portrait and extraction of the damage feature The PSI involves several steps. The measured responses are first decomposed into a number of low dimensional IMFs through EMD, so that the IMFs are amenable to low dimensional embedding, following the Taken’s embedding theorem [27,28]. The Singular System Analysis (SSA) is performed on the so-obtained IMF(s) in order to assess their dimensionality for embedding. The lag coordinate for their phase space reconstruction is estimated from the Average Mutual Information Function (AMIF) of the respective IMF(s). The reconstructed phase portraits from the damaged and undamaged structures are analyzed for the CPST and correlated with the different damage scenarios in the structure. Details of these are expanded in the subsequent sections.

252

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

2.1. Reconstruction of the phase portrait from the measured responses The dynamical states of structure can be geometrically described in phase space, in which, all possible states of a system are represented by a set of points. The phase space include all possible combinations of displacement and velocity because the dynamic evolution can be fully predicted by using these quantities. In a structure in motion, the excitation mechanisms are balanced by the dissipative mechanisms to result into an invariant subspace in the phase space, referred as Phase portrait. Depending on the type of system and excitations, the phase portrait can be a point, a set of points, a curve, a manifold or a more complicated object. The simplest example of a phase portrait is that of an undamped spring mass system under free vibration. The phase portrait takes the form of an ellipse. With viscous damping, the portrait becomes a spiral that gradually reduces its radius to merge into the origin. Nonlinear system under deterministic excitations may result into a limit cycle or more complex form of phase portrait due to chaotic behavior of the system [29,30]. Construction of a phase portrait require measurement of a number of response variables, which might not be practically feasible. However, it is possible to reconstruct the qualitative behavior of the attractor based on limited (even by a single) measurements following the Taken’s embedding theorem [27,28]. A measured time series fxðtÞg (e.g. acceleration) over time steps nt, fx1 ; x2 ; . . . ; xnt g can be embedded in a (pseudo)-phase space of dimension m with appropriate lag ðsÞ as fxi ; xiþs ; . . . ; xiþðm1Þs g. The parameters, embedding dimension ðmÞ and the lag coordinate ðsÞ are determined from the SSA and the AMIF of the measured time series, respectively. The SSA forms a circulant matrix from the measured response as

2

x1

x2

x2

x3

:

:

xNqþ1

xNqþ2

6 6 ½A ¼ 6 4

...

xq

3

. . . xqþ1 7 7 7 : 5 ...

ð1Þ

xN

in which ðqÞ defines a time window ðq ¼ f s sf Þ obtained by the product of the sampling frequency ðf s Þ of the response signal  and another parameter sf that depends on the maximum frequency content ðf Þ of the signal. This is estimated by ensuring that

sf 6 1=f 

ð2Þ 

as proposed by Broomhead and King [31]. The limiting frequency f is the maximum cut-off frequency above which the signal contains no significant power. Matrix ½A can be decomposed by the Singular Value Decomposition (SVD) as

SVDðAÞ ¼ S

X

UT

ð3Þ

where the decomposition consists of a square matrix R, the diagonals of which contain the singular values ðri ; i ¼ 1; 2; . . . ; qÞ. The matrices S and U are the unitary matrices. The singular values are normalized by their sum as

r^ i ¼ ri

.X

ri

ð4Þ

Singular values close to zero are due to noise in the system dynamics, while the dominant ones are associated with the active dynamical degrees of freedom. The singular values also represent the share of a singular vector to the variance of the signal. Further, singular values also decide the degree to which the signal explores various directions in phase space, indicated by the singular vectors [31]. A threshold value based on engineering judgment of the problem is chosen and the number of singular values above this are counted as the dimension ðmÞ for embedding. The other parameter for the reconstruction of the attractor is the lag coordinate ðsÞ, which is based on the AMIF of the signal. Although, the autocorrelation function has also been employed, this has shown to cause misleading results for nonlinear time series [32]. The AMIF is defined as a function of the time lag ðTÞ as [33,34]

IðTÞ ¼

X pðxi ; xiþT Þlog2 ðpðxi ; xiþT Þ=pðxi ÞpðxiþT ÞÞ

ð5Þ

nt

where pðxðiÞÞ and pðxði þ TÞÞ are the estimated probability densities of xðiÞ and xði þ TÞ and pðxðiÞ; xði þ TÞÞ is the joint probability density. The time lag corresponds to the first minima of the AMIF and is chosen as the desired lag for embedding. This is because the minima indicate to reduced level of statistical independency. With this choice of delay, it is ensured that the phase portrait is completely unfolded in phase space. Embedding a phase portrait using too less lag result in trajectories being falsely projected onto the top of another. It is demonstrated that the number of points required to accurately describe an attractor scales in power of its dimension, as suggested by Nerenberg and Essex [35]

N ¼ 2m ðm þ 1Þm

ð6Þ

Therefore, a system with modest dimension will require prohibitively large data to accurately represent the underlying dynamics. It is important to note that the Taken’s embedding theorem [30] is for deterministic excitations and if the system is subjected to stochastic excitation, the embedding theorem, in principle, cannot be used. However, for structure under seismic excitations, the dimension cannot be controlled as in case of chaotic interrogation, in which, the structure is excited by pre-

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

253

scribed chaotic waveforms [18–25]. Seismic responses are high dimensional and the reconstruction of attractor based on embedding theorem are not applicable, in principle. However, the embedding theorem is still employed for stochastically driven system [26], especially if the process is narrow band, so that it can effectively be treated as low dimensional. Trickey et al. [34] have demonstrated that exciting a structure with a narrow band stochastic process can result in practically low dimensional responses that can be embedded in phase space with fewer dimensions. In this context, it may be mentioned that Orstavik and Stark [36] have demonstrated that, a seventy dimensional system can be rendered into a four dimensional through optimal modeling. Stark et al. [37] have generalized Takens embedding theorem for stochastically driven system. 2.2. EMD of high dimensional response signals It is emphasized that a signal must be narrow band in order to rely on its phase space embedment and subsequent analysis of the phase portrait. However, the response of structure subjected to broad band seismic excitations are high dimensional and not qualified for the phase space analysis. In the present study, the measured response signal, which are truly high dimensional are decomposed into a number of oscillator modes, (referred as IMFs) through EMD. The procedure of EMD is described by Huang et al. [38,39]. In EMD, all the local extrema of a response signal ðxðiÞ; i 2 1; ntÞ are identified. These local maxima are then connected by a cubic spline in order to obtain the upper envelope fxu ðtÞg. Similarly all the local minima are connected to get the lower envelope fxl ðtÞg. The mean of the upper and lower envelopes is then obtained as fm1 ðtÞg

fm1 ðtÞg ¼

1 ½fxu ðtÞg þ fxl ðtÞg 2

ð7Þ

The difference between the signal fxðtÞg and the mean fm1 ðtÞg provide a vector fh1 ðtÞg as

fh1 ðtÞg ¼ fxðtÞg  fm1 ðtÞg

ð8Þ

The residual vector fh1 ðtÞg is then treated as a ‘‘pro-IMF” and is further iterated (by treating it as the signal itself) following the k number of similar iterations (referred as ‘‘sifting”) as ðk1Þ

ðkÞ

ðkÞ

 m1 ¼ h1

h1

ð9Þ

The sifting procedure is repeated until a measure, the ‘‘Sum of the Difference” (SD) attains value below a certain threshold. The SD at any iterative step k is defined as ðkÞ

SD

nt  2 X  ðk1Þ ðkÞ  ¼  fh1 ðiÞg  fh1 ðiÞg i¼1

,

nt h X

ðk1Þ

fh1 ðiÞg

i2

ð10Þ

i¼1

the value of SD must be less than a predetermined tolerance to stop the sifting iterations [38,39]. Once the convergence is achieved, the first IMF is obtained as ðkÞ

c 1 ¼ h1

ð11Þ

The first IMF ðc1 Þ contain the finest scale information of the signal. The data becomes coarser in higher IMFs. The first IMF can be separated from the rest of the data to obtain a residual signal as

fr 1 g ¼ fxg  fc1 g

ð12Þ

The residual signal fr1 g can be treated as the new data to extract the subsequent IMFs following the procedure detailed above. The extraction of IMFs stop when frg becomes a monotonic function, from which no more IMF can be extracted. Summing up all the nf number of IMFs fcg and the residue frg provides the original signal as

fxg ¼

nf X fcgj þ frg

ð13Þ

j¼1

An IMF represents a simple oscillatory mode, which can be thought of as counterpart to simple harmonic functions but more generalized. The EMD also assume that any data consists of a number of simple intrinsic modes having same number of extrema and zero-crossings. It is obvious that the IMFs contain increasingly coarser information and each IMF is associated with certain characteristic time scale, dictated by the time lapses between the successive extrema, which also ensures the relative low-dimensionality of the component IMFs over the original broad band signal [39–41]. All these IMFs are zero mean and are normalized to have unit variance as

 ¼ fCg =rfCg fCg j j j

ð14Þ

where rfCgj is the variance of the individual IMFs. For their low dimensionality, some of the IMFs are amenable to phase space embedding of their respective share of attractor and subsequent analysis aiming at extraction of damage features. The EMD algorithm as presented suffers from several drawbacks, which are illustrated in brief. It is stated that the cubic spline fitting is performed to obtain the upper and lower envelopes of the response signal. However, in instances, the data

254

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

might have large swings at the ends and fitting of envelopes lead to propagation of these swings into the original signal, leading to ineffective decomposition. This leads to spurious components and subsequent over estimation of the signal energy. Huang et al. [39] has proposed addition of characteristic waves at the both ends of the signal having capacity to contain the wide swings. There are a number of other alternatives, however, the most common trick is to use a tapering window to the signal to offset any swings. For the seismic response signal, employed in the present study, the amplitude builds up and decays gradually and the signals are virtually free from the effects of end swing, and therefore is not a point of concern in the present study. Another problem associated with the EMD is the difficulty in extracting two intrinsic modes with closely spaced/similar frequencies. This is because the mode extraction is sequential in nature and can only be extracted one after the other. Superior algorithms are proposed to remove this problem using adaptive algorithm [40]. For the present study, however, the structural modes are well separated and the resulting IMFs of the responses are well separated as well. The problem of modal confusion does not arise and the conventional EMD algorithm is employed adequately. It is understood that for using natural excitations, such as earthquake, wind and wave loading, the damage feature must be largely irrespective of the nature of the excitations, employed for generation of the structural responses. In fact, the Chaotic interrogation make use of Chaotic waveform to excite the structure and thereby eliminating this issue. For the present study, it may be noted that other than reducing the high dimensionality, usage of the EMD also helps in reducing the wide variability in the frequency contents that is commonly encountered in earthquake ground motions and the respective structural responses. Further, before estimating the CPST, a normalization on the respective response time histories may be performed. This normalization is performed for each IMF of response time history separately with respect to the respective IMF of the exciting ground motion. This scaling on each narrow band first few IMFs, largely helps to eliminate the dependency of the CPST on the specific waveform of the ground motion. However, it must be kept in mind that the dependency can only be minimized by such scaling, but cannot be eliminated completely. Conventional chaotic interrogation assumes that the baseline structural information are available, so that the responses of the pristine structure can be obtained for comparison of the CPST(s) from the damaged and undamaged structure. In line with this, it may be noted that, even if, the dependency on the specific seismic excitation is not eliminated; the PSI methodology under seismic excitations can still be applied, if the baseline information for the pristine structure is available. 2.3. CPST of the phase portrait as the damage sensitive feature The reconstructed phase portrait is analyzed for certain measure that can be associated with the extent of damage in the structure. This is because the dynamics of the structure is bound to change due to the damage and resulting changes in the stiffness of the structure. The disparities in the attractors show up when contrasted with the undamaged one. This is subsequently quantified as the damage sensitive feature. From a number of features (for example LAVR, CPST, NPE, PSW) suggested in the literature to measure the phase portrait dissimilarity, the CPST is adopted herein because of its simplicity. The proposition of CPST [21,22] was originally made for chaotic interrogation of structure. A brief review of the CPST is presented herein. Let us consider a md dimensional phase space trajectory from the damaged and undamaged structure. A point on the damaged portrait is randomly chosen to be the Fiducial point with its coordinates as

fxgtf ¼ f xðtÞ xðt þ sÞ . . . xðt þ ðmd  1ÞsÞ g

ð15Þ

The same point is also marked on the undamaged portrait by using only spatial information on phase space. It must be reminded that this point may or may not be on the trajectory of the undamaged portrait because of the damage induced changes. The objective here is to extract some metric (as CPST) that can distinguish between the undamaged and damaged portrait. The nearest p neighbors of this fiducial point are selected on the undamaged trajectory as fxgpj , based on the Euclidean distance between the neighbors and the fiducial point as

  R fygpj ; fxgtf ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  T   fxgtf  fygpj fxgtf  fygpj

ð16Þ

It may be noted that the nearest neighbor search is purely based on the spatial coordinate of the points in the phase space and nothing to do with their time localization. All these nearest neighbors of the fiducial point on the undamaged portrait are then predicted ahead of time ðtf þ Dt f Þ using the time evolution information. From this predicted neighborhood on the undamaged portrait, the prediction of the Fiducial point after the time instant Dtf can be obtained as

f^xgtf þDtf ¼

1X fy g p j¼1 j tf þDtf p

ð17Þ

This time interval ðDtf Þ for prediction (‘‘prediction horizon”) is adopted in the range of 1 6 Dt f 6 s=2 [25], where, s is the time lag obtained from the AMIF. The number of nearest neighbors ðpÞ is typically chosen in the range of 104 to 103 of the total number of points on the phase portrait [26]. The predicted fiducial point on the undamaged portrait is then mapped back onto the respective point on the damaged portrait, once again based on only spatial information. It might be noted that this

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

255

mapped point on the damaged portrait might not fall on the trajectory itself. The actual evolution of the damaged fiducial   point at time step Dt f ahead on the damaged trajectory is now obtained from the time evolution information fxgtf þDtf . The CPST is then defined as

CPST ¼

1 p

r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T   f^xgtf þDtf  fxgtf þDtf f^xgtf þDtf  xtf þDtf

ð18Þ

which is the Euclidean norm of the difference among the predicted damaged trajectory and actual damaged trajectory around the fiducial point. The calculations for the CPST are repeated for a number of randomly selected fiducial points on the damaged portrait in order to get a reliable estimate of the CPST as

CPST ¼

nf X CPSTðiÞ=nf

ð20Þ

i¼1

The number of fiducial points are estimated by adaptively taking more and more number of points to obtain a stable estimate. It was suggested that the number of fiducial points can be taken as 5% of the total population of points on the attractor in order to get a reasonable estimate [26]. It might so happen that when an undamaged attractor is contrasted with itself for CPST, instead of zero, some small value may appear, which may be termed as the baseline error. This error may be subtracted from the damaged CPST in order to get the change in CPST caused by the damage only. The appealing feature of the approach is that the damage prediction is entirely data driven and is irrespective of the adopted structural model. Once a structure is damaged then the value of the CPST increase monotonically with the severity of damage [21,22]. A flow chart of the proposed methodology is depicted in Fig. 1 by explaining the various steps involved. The duration of ground motions, adopted herein are typically in the range of 30–40 s. However, ground motions of shorter duration are also common, for instance 20 s or so. In such cases, the length of the data may not be adequate to populate the embedded phase space in order to obtain a reliable estimate of the CPST. In anticipation of small duration earthquakes (such as aftershocks after a main event), the sampling interval of the data acquisition may be reduced so that a larger number of samples can be acquired to meet the need of phase space analysis. The modern Data Acquisition Systems can do this efficiently. However, in case, the data length is short, additional data points can be generated using interpolation between the successive data points.

3. Numerical illustration The approach is illustrated numerically with respect to a benchmark structure with varying degree of damage(s), suggested by the IASC-ASCE task group of SHM. The details of the benchmark building model and the IMF-CPST damage detection methodology is illustrated subsequently.

Fig. 1. Flow chart of the proposed IMF-CPST methodology for damage detection.

256

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

3.1. IASC-ASCE benchmark building The IASC-SCHM benchmark building [42] is a four-story, two-bay by two-bay laboratory scale (one-third scaled) steel frame originally fabricated in the Earthquake Engineering Research Laboratory at the University of British Columbia (UBC) as a part of the IASC-ASCE SHM initiatives. A schematic of the model is shown in Fig. 2. The model building is 2.5 m  2.5m in plan and 3.6 m in elevation. The members are hot rolled 300 W steel with a nominal yield stress of 300 MPa. The sectional details of the members in the model can be obtained from Table 1 of the study reported by Johnson et al. [42]. On the exterior face of each panel, there are two diagonal braces (K-braces) that are progressively removed/modified in order to emulate damage. There is one floor slab in each bay of each floor. The mass of each slab at the first floor is 800 kG, in second and third floor, these are 600 kG each. On fourth floor, either four 400 kG slabs or three 400 kG and one 550 kG slabs are placed to create symmetric and asymmetric structure, respectively. Presently, the symmetric distribution of the floor masses are adopted. Two alternative finite element models of this benchmark structure have been proposed to generate simulated responses as input for damage detection algorithm. The one that is adopted presently is a 12 degree of freedom (DOF) shear building model that constrains all motions, except two translational and one rotational motion at each floor. The beams and columns are modeled as Euler-Bernoulli beam element. The braces are modeled as axially loaded element with no bending stiffness. Details of this 12-DOF model (mass and stiffness matrices) are provided in literature [42]. The natural frequencies of the model in various modes of vibration are provided in Table 6 in Ref. [42] that are employed here for verification. The xdirection is the strong direction due to orientation of the column. The weak direction (+ve y) is referred as south and strong direction (+ve x) as west. Six damage scenarios are simulated in the benchmark building in order to verify the efficacy of any SHM methodology in detecting, localizing and quantifying damage(s). The detection and localization are investigated presently. The six damage scenarios are described. In case 1, the braces are removed from the first story to reduce the respective lateral stiffness. In case 2, the braces in the first and third stories are removed to reduce stiffness in these stories. In case 3, only one brace in the first story (north brace on the west face of the structure) is removed. In case 4, one brace in the first story (north brace on the west face) and one brace in the third story (west brace on the north face) are removed. Case 5, is identical to case 4 but the joint between the north floor beam at the first story on the west face of the structure is partially weakened (by removing the screw) from the north-west column. The beam-column connection can then only transmit forces but cannot transmit moment. In case 6, the stiffness in one brace (north brace on the west face) in the first story is reduced one third ð1=3Þ of its original stiffness.

Y(south)

X (west)

Fig. 2. Schematic of the IASC-ASCE SHM benchmark structure.

257

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265 Table 1 Percentage loss in the lateral storey stiffness corresponding to various damage scenarios in 12-DOFs benchmark building. Element

Damage scenarios (percentage change in lateral stiffness)

Story

DOFs

(1)

(2)

(3)

(4)

(5)

(6)

1 1 1 2 2 2 3 3 3 4 4 4

x y h x y h x y h x y h

45.24% 71.03% 64.96% 0 0 0 0 0 0 0 0 0

45.24% 71.03% 64.96% 0 0 0 45.24% 71.03% 64.96% 0 0 0

0 17.76% 9.87% 0 0 0 0 0 0 0 0 0

0 17.76 9.87 0 0 0 11.31% 0 9.16% 0 0 0

0 17.76% 9.87% 0 0 0 11.31% 0 9.16% 0 0 0

0 5.92% 2.88% 0 0 0 0 0 0 0 0 0

The above six damage scenarios are inflicted to the structure. The respective percentage loss in the horizontal story stiffness, owing to these damage cases are tabulated in reference [39] and are reproduced in Table 1. It may be pointed out that an additional case (Case 0) is added to refer to the pristine structure without any damage. It may be noted that these damage scenarios do not follow particular trend while going from case 1 to case 6. However, a comparison among the reduction in lateral stiffness can provide an idea about the relative severity of these specific scenarios. Detection of a relatively less severe damage case (as indicated in Table 1) by the SHM methodology is more stringent than the more severe cases. The variation of the stiffness reductions for these scenarios are graphically shown in Fig. 3. Due to a number of element damage in each case, the square norm of these stiffness reductions are used in this plot. The values of the norms are further normalized with respect to the maximum reduction, in order to get a maximum value of 1 for most severe damage case, as indicated in Fig. 3. The stiffness variation corresponding to these damage scenarios are shown in sequence, which will be useful to assess the efficacy of the damage feature. This is because the variation in the damage feature value with these changing sequence of damage must corroborate with the observed trend in Fig. 3. A general requirement of a legitimate damage feature is that the parametric trend of variation of damage feature and stiffness variation with varying degree of damage(s) must be identical or similar. Three recorded earthquake ground motions pertaining to real earthquakes are employed as input excitations to excite the benchmark structure. Details of these motions are listed in Table 2. This motions are applied along the weak (south) direction of the structure for exciting the structure. The equation of motion of the structure subjected to ground acceleration f€ xg ðtÞg is written as

€ g þ ½Cfug _ þ ½Kfug ¼ ½Mfrg€xg ½Mfu

ð21Þ

€ g; fug _ and fug are the relative acceleration, where ½M; ½C and ½K are the mass, damping and stiffness matrix, respectively, fu velocity and displacement vectors pertaining to the degrees of freedom of the structure, respectively. The vector frg is the influence coefficient vector and € xg is the ground motion due to earthquake. € g at specific degrees of freedom can be simulated from these responses The noisy sensor measurements of accelerations fv by adding noise. The measured responses can be expressed as

_ þ frg€xg þ fwg fv€ g ¼ ½C 1 fug þ ½C 2 fug

ð22Þ

Fig. 3. Variation in the stiffness reductions pertaining to the damage scenarios in the IASC-ASCE benchmark building (Damage scenario ‘‘0” refers to the pristine structure).

258

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

Table 2 Recorded earthquake ground accelerations used as input excitations. Ground motions

Event

Station

Peak ground accelerations (PGA)

GM-1 GM-2 GM-3

1992 Imperial Valley 1995 Kobe 1989 Loma Prieta

El Centro Takatori Treasure Island

0.31 g 0.61 g 0.16 g

in which fwg is the sensor noise vector, the elements of which are Gaussian pulse processes with root mean square (r.m.s) as percentage of the r.m.s acceleration and ½C 1  and ½C 2  are two matrices to express the physical degrees of freedom for measurements in terms of the active degrees of freedom. Eq. (21) is integrated in time by using step by step integration to obtain the response(s) of the structure at specific degrees of freedom. The time step adopted for integration satisfies the stability and accuracy criteria. The simulated response(s) of the structure with different damage scenarios are employed as the measured responses for subsequent analysis.

3.2. Numerical illustration on the IASC-ASCE benchmark building The seismic phase space interrogation procedure as discussed above is numerically illustrated by taking the IASC-ASCE benchmark building with seismic excitations as input. The sectional properties of the building are adopted from the literature. Time step integration with time step ðDt ¼ 0:001sÞ is adopted to integrate the equation of motion to obtain the responses of interest. It is reminded that in actual practice, these responses are measured by sensors (e.g. accelerometers) installed at various locations of the structure. Typical acceleration time history at a node at the first floor of the structure is shown in Fig. 4a considering GM1 excitation. The acceleration response is visibly wide band and is decomposed into a number of IMFs by the EMD, as shown in Fig. 4b–g. The acceleration time history is high dimensional, whereas, the higher IMFs are increasingly low dimensional. First few IMFs (IMFs 1, 2) are also observed to be high dimensional, whereas, coarser ones (IMFs 3, 4. . .) can sufficiently be treated as low dimensional. It is understood that these IMFs do not retain full evolution of dynamics as in the original signal, rather, they can be treated as reduced order representation of the original response signal. While using the IMFs for damage feature extraction, it is implied that, even with such reduced order information, specific IMFs will still retain significant dynamics so that the damage induced changes will be reflected in the extracted damage feature from them. More precise estimate of the low dimensionality of the IMFs can be obtained from the results of the SSA. The SSA results are presented in terms of the singular values vs. embedding dimension plots in Fig. 5a–d for the acceleration time history and its 1-st, 3-rd and 9-th IMFs, respectively. The SSA spectra clearly indicate how the required embedding dimension for the phase portrait reconstruction are reduced for the higher IMFs. Particularly, the significant reduction in the required dimensions for the third and higher IMFs are noteworthy. It may be reminded that the significant singular values represent the number of active dynamical degrees of freedom in the system that must be retained for the phase space embedding. It is also noted that as one move from finer to coarser (or lower to higher) IMFs, the required dimension for embedment reduces significantly. It might appear that the coarsest IMF might be the best choice as far as the low dimensionality is concerned. However, this is not the sole criteria because the respect IMF should also retain sufficient dynamics so that it is still able to portray the changes in dynamics owing to the inflicted damages. This require a balance between the low-dimensionality and the extent of the retained dynamics. However, these are two mutually conflicting requirements and an optimal choice of the IMF should be made based on such conflicting objectives. Without having a more rigorous analysis of the pareto-optimality, it is generally observed that the 3-rd and/or 4-th IMFs are the best choice to accommodate both the requirements. The damage detection can be based on both of these IMFs to infer on the detection and localization. It worth mention that the optimal choice of IMF may be based on a more rational quantitative criteria. However, this is not done in the present analysis and the choice is rather based on the trend of the numerical results. However, this may be taken up in future studies. The lag coordinate for the embedment is obtained from the first minima of the respective AMIF, which are plotted for the first and third IMFs in Fig. 6a and b respectively. Whereas the minima are clearly indicated from the first IMF, the AMIF for the third IMF is flat, which render the lag coordinate of the embedding for the higher IMFs (e.g. IMF-3) much larger than the lower ones. Once the embedding dimension and the lag coordinate are obtained, the IMFs are embedded in phase space. The reconstructed phase portrait from the third IMF of the acceleration time history is shown in Fig. 7 using a three dimensional embedding. It must be mentioned that, although a three dimensional embedment is shown for convenience, it is not necessary that the computation are always based on such reduced embedment. The CPST value from the reconstructed phase portrait of an IMF for different sensor locations are shown in a three dimensional bar chart in Fig. 8a–c. The value of the prediction horizon adopted in this study is (Dtf ¼ 2Dt) s, which is two times of the sampling interval of the response signal(s). The plot in Fig. 7 is referred as CPST map. The maps are also obtained from the first three IMFs, shown in Fig. 8a– c, respectively. The variations of the CPST with the benchmarked damage scenarios must corroborate with the stiffness variation with respect to the changing damage scenarios, shown in Fig. 2 earlier. For a good damage feature the trends from both should be identical.

259

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

(a)

(b)

(c)

IMF-1

IMF-2

(d)

(e)

IMF-4

IMF-3

(f)

(g)

IMF-8

IMF-9

Fig. 4. (a) Acceleration time history at a node of the first story of the benchmark structure and respective, (b) IMF-1, (c) IMF-2, (d) IMF-3, (e) IMF-4, (f) IMF-8 and (g) IMF-9 obtained from the EMD.

(a)

Acceleration time history

(c)

IMF-3

(b)

(d)

IMF-1

IMF-9

Fig. 5. Singular value spectrum from the SSA of (a) acceleration, (b) 1st, (c) 3-rd, (d) and (e) 9th IMFs.

Out of all the trends portrayed by the CPSTs from all sensor locations and with varying combination of damage scenarios (Fig. 8a–c), the one that most closely mimics the parametric trend of damages (as shown in Fig. 2), is observed to be the IMF3. The CPST map from the other two IMFs show significant disparity from that in Fig. 3 and may not be a good damage feature. Thus, it can be concluded that IMF-3 is the best choice for this purpose. In continuation with the previous discussion, it

260

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

Fig. 6. AMIF for the 1-st and 3-rd IMFs derived from the acceleration time history.

Fig. 7. Phase portrait embedded in three dimension from the third IMF.

can be commented that the IMF-3 is the low dimensional reduced representation of the response signal but it retains sufficient system dynamics to reflect the significant changes due to damage(s). The IMF-3 becomes the optimal choice that perhaps optimally accommodates the conflicting objectives of low-dimensionality and damage sensitivity. It may be noted that such choice may not be universal and might vary slightly depending on the system and excitation. An indication regarding the location of the damage can also be derived from the CPST map. The peak of the CPST map, obtained from a particular sensor and its subsequent decay at locations far away from the site of damage helps in localizing the damage. The damage localization is explained with the help of Figs. 9 and 10. Figs. 9 and 10 indicate the localized damage scenarios in the benchmark building and the respective variations observed in the CPST values from all the sensors, respectively. The damage locations in Fig. 9 are corroborated with the sensor-wise-CPST-variation pattern, shown in Fig. 10. The location of damage at the first storey (damage case 1) is reflected in the respective CPST variation attaining a peak at the sensor location 1, thereby, indicating the damage at the first story level. Damage pattern 2 have damage both at the first as well as third story level, those are also reflected in the respective variation observed in the CPST(s), indicating peak at the first sensor location (corresponds to first story damage) and a subsequent increase in the CPST value at the third story level (indicating damage at the third story). Further, comparison of the CPST(s) among damage scenario 1 and 2 clearly indicates more severe damage for the latter (scenario 2). Similarly, damage location 3 is reflected by the peak around sensor 1, and the less severity

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

261

Fig. 8. Bar chart showing the CPST map for varying combinations of damage scenario and sensor location. The maps are obtained from the (a) first, (b) second and (c) third IMFs of the acceleration time history. (The zero-th scenario refers to the pristine structure).

(over damage 1) is also indicated by the lower value of the CPST. However, the comparison appears critical among damage cases 4, 5 and 6. The less severe nature of the damage pattern 6 is reflected in terms of the lowest value of the CPST peak at sensor 1, which also indicate the location of damage at the first storey. Damage pattern 4 and 5 show almost identical trend of variation in the CPST(s). Both shows peak at sensor 1, indicating the location of damage at the first storey. Accompanying damages at the third story are also indicated by the increasing trend of CPST at sensor locations 2 and 3. However, in the current plot (even with the existing discrepancies among the CPSTs from the damage case 4 and 5), it is almost impossible to localize the loosening of bolts in second story in damage scenario 5. However, such precise localization and more local information can be obtained by increasing the number of sensors, that will provide more finer resolution about the spatial variation of CPSTs. Thus, with limited number of sensors and their locations, the CPST feature is able to correctly point out to the localization of damage, at least globally. In next, the robustness of the reduced CPST is assessed by considering different degrees of noise contamination in the response signals at various location of the sensors. This is because the measured responses are often contaminated by noises from various sources, such as thermal noise, mechanical relaxation and electrical noise from the data acquisition systems. A robust damage feature must be tolerant to moderate degree of noise, so that its efficiency does not get impaired even with the occurrence of slightest noise in the signal. The degree of contamination is expressed in terms of the Signal to Noise ratio (SNR), expressed as

SNR ¼ 10log 10 ðrs =rn Þ2 ¼ 20log 10 ðrs =rn Þ

ð23Þ

where rs and rn are the root mean square (r.m.s) value of the signal ðsðtÞÞ and the noise ðnðtÞÞ. The SNR is typically expressed in deci-bel (dB) unit. Typical values of the SNR considered herein are 20, 25 and 30 dB and the respective variations of CPST under different damage scenarios are shown in Fig. 11. The noise is of Gaussian type with different values of its r.m.s as dictated by the respective dB value. Because of its Gaussian nature, a wide range of frequencies are present in the noise. While the Gaussian nature of the noises take care of the frequencies, the amplitude of the noise is dictated by their r.m.s value. It is observed that the trend of variation in the CPST does not change with varying degree of noise, which implies that the occur-

262

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

Fig. 9. Location of the damage (indicated in red) in the benchmark building and the location of the sensors (yellow star). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Variation in CPST along the location of measurements (sensors), shown in Fig. 9.

rence and localization of damage will not be affected by the presence of noise. However, it is also observed that the CPST value changes remarkably due to addition of noise. This may affect the quantification of damage, which require further study. The efficacy of the CPST feature is also studied using alternative seismic ground motion scenarios listed in Table 2. The variations of the CPSTs with varying damage scenarios, as obtained from different IMFs are shown in Fig. 12a–i. While plotting, the IMFs that offer best performance are chosen. This happens to be mostly the third and fourth IMFs for each of these motions. However, slight variations among the trend of CPST (with varying damage scenario) obtained from the IMFs are

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

263

Fig. 11. Sensitivity of the CPST under varying degree of noise contamination (damage scenario ‘‘0” refers to pristine structure).

observed. For instance, under GM1 and GM 2, the third IMFs perform best, whereas the fourth IMF performs much better than the others in case of GM-3. Needless to mention that, the relative superiority of one IMF on another is based on the similarity of their respective CPST variations with the extent of damages. It can generally be concluded that the third IMF should be the optimal choice in absence of any specific information. Fig. 12a–i also indicate to the robustness of the CPST feature in detecting damages under varying earthquake with widely varying characteristics in their frequencies as well as amplitudes. 3.3. Localization of damage To further illustrate the damage localization capability of the CPST_IMF feature, a 12 storied shear building is considered in Fig. 13a. The 12 storied shear building have uniform mass of 1000 kG lumped at each storey level. The storey stiffness of each story is 15,000 kN/m. With this the fundamental time period comes around 1.29 s. The damage is induced by reducing (20% reduction) the lateral story stiffness along different stories. It is assumed that the acceleration time histories can be obtained from the accelerometers mounted along these stories and the CPST variation from the responses are obtained. The variations in the CPST along the stories are shown for two damage cases. Firstly, the damage is induced at the first story and next when the damage is induced at the eighth story. These are shown in Fig. 13b and c, respectively. The occurrence of damage(s) at the first and eight story are indicated by the CPST peaks at the respective locations, as shown in Fig. 13b and c, respectively. Thus the damage locations corroborate with the location of the CPST peaks and thereby helps in localizing the damage.

Fig. 12. Performance of the CPST for varying damage scenarios as obtained from the (a) IMF-2 of GM-1, (b) IMF-2 of GM-2, (c) IMF-2 of GM-3, (d) IMF-3 of GM-1, (e) IMF-3 of GM-2, (f) IMF-3 of GM-3, (g) IMF-4 of GM-1, (h) IMF-4 of GM-2 and (i) IMF-4 of GM 3.

264

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

Fig. 13. (a) A 12-storied shear building and the variations in CPST along the stories when damage is inflicted at the (b) first and (c) at eighth storey.

4. Conclusion The phase portrait dissimilarity based chaotic interrogation methodology for structural damage detection is extended for earthquake excitations to eliminate the need for exciting the structure with chaotic waveform. The usage of EMD is proposed to overcome the difficulty associated with the high dimensionality of the seismic responses. The EMD allows decomposition of high dimensional response signal into a number of low dimensional IMFs, containing reduced order information of the dynamics, yet retaining sufficient enough information to reflect the changes in the dynamics due to incipient damage(s). The optimal choice of the IMF(s) based on the conflicting objectives of low dimensionality and the sufficiency of information for damage sensitivity are proposed. The suitability of the third/fourth IMFs for this purposes, along with their slight variations from one excitation to another are revealed. The study demonstrates that the CPST of such reduced order portrait still retains the ability to indicate the presence, severity as well as localization of damage. The procedure is numerically demonstrated in reference to the IASC-ASCE benchmark building. The CPST_IMF feature is shown to be tolerant to significant level of noise contamination. The robustness of the feature is also assessed by considering alternative ground motions from different earthquakes.

References [1] C. Yun, K. Choi, S. Kim, Y. Song, Comparative evaluation of nondestructive test methods for in-place strength determination, in: H.S. Lew (Ed.), Nondestructive Testing, ACI SP-112, ACI, Detroit, 1988, pp. 111–136. [2] J.P. Lynch, K.H. Loh, A summary review of wireless sensors and sensor networks for structural health monitoring, Shock Vib. Digest 38 (2) (2006) 91– 128. [3] Charles R. Farrar, Keith Worden, An introduction to structural health monitoring, Philos. Trans. Roy. Soc. London A: Math., Phys. Eng. Sci. 365 (1851) (2007) 303–315. [4] R. Ghanem, M. Shinozuka, Structural-system identification. I: Theory, J. Eng. Mech., ASCE 121 (2) (1995) 255–264. [5] H.G. Natke, System identification techniques, Struct. Saf. 3 (3–4) (1986) 269–316. [6] C.B. Yun, M. Shinozuka, Identification of nonlinear structural dynamic systems, J. Struct. Mech. 8 (1980) 2. [7] H. Katkhuda, A. Haldar, A novel health assessment technique with minimum information, Struct. Contr. Health Mon. 15 (2008) 821–838. [8] W.X. Ren, H.B. Chen, Finite element model updating in structural dynamics by using the response surface method, Eng. Struct. 32 (8) (2010) 2455– 2465. [9] F.E. Udwedia, Some uniqueness results related to building structural identification, SIAM J. Appl. Math. 34 (1) (1978) 104–118. [10] F.E. Udwadia, Some nniqueness results related to soil and building structural identification, SIAM J. Appl. Math. 45 (4) (1985) 674–685. [11] O.S. Salawu, Detection of structural damage through changes in frequency: a review, Eng. Struct. 19 (9) (1997) 718–723. [12] A.K. Pandey, M. Biswas, M.M. Samman, Damage detection from changes in curvature mode shapes, J. Sound Vib. 145 (2) (1991) 321–332. [13] K. Krishnan Nair, A.S. Kiremidjian, K.H. Law, Time series-based damage detection and localization algorithm with application to the ASCE benchmark structure, J. Sound Vib. 291 (1-2) (2006) 349–368. [14] Z. Hou, M. Noori, R. Amand, Wavelet-based approach for structural damage detection, J. Eng. Mech. 126 (7) (2000) 677–683. [15] S.W. Doebling, C.R. Farrar, M.B. Prime, D.W. Shevitz, A review of damage identification methods that examine changes in dynamic properties, Shock Vib. Digest 30 (2) (1998) 91–105. [16] H. Sohn, C.R. Farrar, F.M. Hemez, J.J. Czarnecki, D.D. Shunk, D.W. Stinemates, B.R. Nadler, A review of structural health monitoring literature: 1996– 2001. Los Alamos National Laboratory Report, LA-13976-MS, 2003. [17] J.M. Nichols, M.D. Todd, M. Seaver, Use of chaotic excitation and attractor property analysis in structural health monitoring, Phys. Rev. E 67 (2003) 016209. [18] M.D. Todd, J.M. Nichols, L.M. Pecora, L.N. Virgin, Vibration based damage assessment utilizing state space geometry changes: local attractor variance ratio, Smart Mater. Struct. 10 (2001) 1000–1008. [19] J.M. Nichols, S.T. Trickey, M.D. Todd, L.N. Virgin, Structural health monitoring through chaotic interrogation, Meccanica (2003) 239–250. [20] D. Chelidze, Ming Liu, Multidimensional damage identification based on phase space warping, Nonlinear Dyn. 5 (2006) 9007.

B. Paul et al. / Mechanical Systems and Signal Processing 91 (2017) 250–265

265

[21] Z. Nie, H. Hao, H. Ma, Structural damage detection based on the reconstructed phase space for reinforced concrete slab: experimental study, J. Sound Vib. 332 (2013) 1061–1078. [22] Z. Nie, H. Hao, H. Ma, Using vibration phase space topology changes for structural damage detection, Struct. Health Mon. 11 (5) (2012) 538–557. [23] G. Liu, Z. Mao, M. Todd, Z. Huang, Damage assessment with state-space embedding strategy and singular value decomposition under stochastic excitation, Struct. Health Mon. 13 (2) (2014) 131–142. [24] L.A. Overbey, C.C. Olson, M.D. Todd, A parametric investigation of state-space-based prediction error methods with stochastic excitation for structural health monitoring, Smart Mater. Struct. 16 (2007) 1621–1638. [25] S. Torkamani, E.A. Butcher, M.D. Todd, G. Park, Hyperchaotic probe for damage identification using nonlinear prediction error, Mech. Syst. Signal Process. 29 (2012) 457–473. [26] J.M. Nichols, Structural Health Monitoring of offshore structures using ambient excitation, Appl. Ocean Res. 25 (2003) 101–114. [27] H. Whitney, Differentiable manifolds, Ann. Math. 37 (1936) 645–680. [28] F. Takens, Detecting strange attractors in turbulence, Lecture Notes in Mathematics, vol. 898, Springer, Berlin, 1981, 366. [29] R. Rand, Lecture Notes on Nonlinear Vibrations, . [30] S.H. Strogatz, Nonlinear Dynamics and Chaos, Perseus, New York, 1994. [31] D.S. Broomhead, G.P. King, Extracting qualitative dynamics from experimental data, Physica D 20 (1986) 217–236. [32] H. Tong, Non-Linear Time Series Analysis, Oxford University, UK, 1990. [33] A.M. Fraser, H.L. Swinney, Independent coordinates for strange attractors from mutual information, Phys. Rev. A 33 (2) (1986) 1134–1140. [34] T. Trickey, M.D. Todd, J.M. Nichols, M. Seaver, Considerations for attractor property analysis in vibration based damage detection, Proc. SPIE 5047 (2003) 340–349. [35] M.A.H. Nerenberg, C. Essex, Correlation dimension and systematic geometric effects, Phys. Rev. A 42 (1990) 7065–7074. [36] S. Orstavik, J. Stark, Reconstruction and cross-prediction in coupled map lattices using spatio-temporal embedding techniques, Phys. Lett. A 247 (1998) 145–160. [37] J. Stark, D.S. Broomhead, M.E. Davies, Huke J. Takens, Embedding theorems for forced and stochastic systems, Nonlin. Anal., Theor., Meth. Appl. 30 (8) (1997) 5303–5314. [38] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, The empirical mode decomposition and hilbert spectrum for nonlinear and nonstationary time series analysis, Proc. R Soc. London Ser. A 454 (1998) 903–995. [39] N.E. Huang, Z. Shen, R.L. Steven, A new view of nonlinear water waves: the Hilbert Spectrum, Annu. Rev. Fluid Mech. (1999) 417–457. [40] K. Dragomiretskiy, D. Zosso, Variational mode decomposition, IEEE Trans. Signal Process. 62 (3) (2014) 531–544. [41] A. Kunwar, R. Jha, M. Whelan, K. Janoyan, Damage detection in an experimental bridge model using Hilbert-Huang transform of transient vibrations, Struct. Contr. Health Mon. 20 (2013) 1–15. [42] E.A. Johnson, H.F. Lam, L.S. Katafygiotis, J.L. Beck, Phase I IASC-ASCE structural health monitoring benchmark problem using simulated data, J. Eng. Mech., ASCE 130 (2004) 3–15.