# A Study on Instantaneous Time-Frequency Methods for Damage Detection of Nonlinear Moment-Resisting Frames

Authors: Ehsan Darvishan; Gholamreza Ghodrati Amiri; Pedram Ghaderi
Journal: Shock and Vibration
 (2014)

---

## Abstract

Most of the civil structures exhibit nonlinear hysteresis behavior during earthquakes. However, detection of damage in these structures is a challenging issue due to successive change in structural characteristics during a seismic excitation. The current paper presents a promising approach for damage detection of nonlinear moment frames. First, several instantaneous time-frequency methods including Hilbert-Huang transform, direct quadrature, Teager energy operator, and higher-order energy operator are investigated as signal processing tools and the most appropriate method is selected using an outlier analysis. Next, a procedure is proposed based on time-frequency analysis in conjunction with clustering to find damage extension in moment frames under a seismic excitation using frequency, amplitude, and energy damage measures. A probabilistic approach is implemented to investigate capability of the procedure for different ground motion records using incremental dynamic analysis. Results show that frequency is not an appropriate feature to detect damage in nonlinear structures.


---

## Body

## 1. Introduction

During an earthquake many civil structures undergo damages. Therefore, it is important to ensure serviceability and safety of the structures after a seismic excitation. Today, nondestructive evaluation techniques (NDE) are widely utilized for health monitoring and damage detection of structures. At present, current practical methods are based on visual inspection, CT scanning, ultrasonic, stress waves, acoustic emission, and so forth [1]. However, these local damage detection methods are effective only for small structures or structural members. Recently, global vibration-based techniques have attracted researchers since they can solve the problem of large and complicated structures. These methods include power spectrum, Fourier transform, spectrum analysis, and cepstrum analysis [2]. The main drawback of these procedures is that they cannot trace time-varying nature of signals due to their fixed functions. Instead, time-frequency methods are effective tools for analyzing nonstationary signals which makes them appropriate for health monitoring and fault detection of civil, mechanical, and aerospace systems [3].Doebling et al. [4] and also Fan and Qiao [5] gave a detailed review of vibration-based damage detection methods. An experimental investigation was conducted on a three-story shear building by Xu and Chen [6]. They studied applicability of empirical mode decomposition (EMD) to detect sudden change in structural stiffness. Yang et al. [7] tried to extract damage spikes from vibration response of the ASCE benchmark building using EMD. Rezaei and Taheri [8] proposed a new damage index for damage detection of beams using EMD. Chen [9] experimentally investigated the capability of EMD for structural health monitoring of large civil structures. Hilbert-Huang transform (HHT) was utilized as a signal processing tool to detect existing damages in an actual benchmark problem by Liu et al. [10]. Lin and Chu [11] characterized the acoustic emission signals due to cracks in welded parts of an offshore structure by application of HHT. Hsu et al. [12] proposed a new damage index to detect damage in a steel structure at early state. Xu et al. [13] employed HHT to identify modal properties of Di Wang tower in China. Bouchikhi et al. [14] compared EMD and B-spline EMD associated with HHT and THT for multicomponent nonstationary signals.However, despite extensive investigations on structural damage detection, most researchers have focused on damage detection of elastic structures by monitoring stiffness change due to damage and researches on nonlinear hysteretic structures are scarce in the literature. Saadat et al. [15] proposed an intelligent parameter varying procedure to identify damage in a nonlinear mass-spring system. Xia [16] proposed a new approach to detect time-varying characteristics of a nonlinear structure with hysteresis behavior.This study focuses on detection of damage in nonlinear moment frames due to seismic excitations. This paper is organized into two parts. In the first part, capability of several well-known time-frequency procedures to extract signals is investigated using acceleration response of the 3-story SAC steel structure [17] under a seismic excitation. An outlier analysis is performed to investigate the adverse effect of outliers on statistical parameters of extracted frequencies. The second part presents a promising approach to detect damage due to hysteretic behavior of moment connections of a structure under seismic excitations. Three damage measures including frequency, amplitude, and energy are compared. A probabilistic approach is implemented to quantify damage under various ground motion records.


## 2. Signal Processing Procedures

### 2.1. Empirical Mode Decomposition

In 1996, Huang et al. [18] proposed empirical mode decomposition (EMD) for nonlinear and nonstationary data. EMD is an adaptive procedure which relies on two basic assumptions: (
1
) the number of extrema and zero-crossings must be equal or differ at most by one for all the data; (
2
) mean value of envelope defined by local maxima and minima must be zero at each point. Extraction of intrinsic mode functions (IMFs) is performed via sifting process described below.First, determine local extrema of the data. Plot an upper envelope by using a cubic spline passing all the extrema,e
u
(
t
). Repeat the same procedure for minima to find e
l
(
t
). Calculate their mean, m
1
(
t
), as
(1)
m
1
(
t
)
=
e
u
(
t
)
+
e
l
(
t
)
2
.The difference between signal andm
1
(
t
) is designated as the first component, h
1
(
t
):
(2)
h
1
(
t
)
=
x
(
t
)
-
m
1
(
t
)
.Now, the first component is considered as the signal and the above procedure is repeatedk times to satisfy the basic conditions:
(3)
h
1
k
(
t
)
=
h
1
(
k
-
1
)
(
t
)
-
m
1
k
(
t
)
.Thus,h
1
k
(
t
) will be the first IMF:
(4)
c
=
h
1
k
(
t
)
.A stopping criterion is computed at each repeat as follows:(5)
SD
=
∑
t
=
0
T
[
|
h
(
k
-
1
)
(
t
)
-
h
(
k
)
(
t
)
|
2
h
2
(
k
-
1
)
(
t
)
]
.When the stopping criterion is satisfied the first IMF,c
1, is obtained. To obtain other IMFs the residue, r
1, is considered as the signal and the procedure is repeated:
(6)
r
1
=
x
(
t
)
-
c
1
.When the residue is constant or a monotonic signal or no more sifting is possible the procedure is stopped. Finally, we can reconstruct the original data by summing up(7)
x
(
t
)
=
∑
j
=
1
n
c
j
+
r
n
.


### 2.2. Hilbert-Huang Transform

Hilbert-Huang transform (HHT) implements the benefits of EMD and Hilbert spectral analysis (HSA) to extract a signal to its intrinsic mode functions [19]. Hilbert transform of a signal, x
(
t
), is an integral transform that can be thought of as the convolution with the Cauchy kernel
(8)
H
[
x
(
t
)
]
=
1
π
PV
∫
-
∞
+
∞
x
(
t
)
t
-
τ
d
τ
,
where PV is the Cauchy principal value. Frequency content of a signal, x
(
t
), can be obtained by providing the conjugate part to give an analytic complex signal as
(9)
z
(
t
)
=
x
(
t
)
+
i
y
(
t
)
=
a
(
t
)
e
i
θ
(
t
)
.Thus(10)
a
(
t
)
=
x
2
(
t
)
+
y
2
(
t
)
θ
(
t
)
=
arctan
(
y
x
)
,
wherea(t) is the instantaneous amplitude and θ
(
t
) is the instantaneous phase. The instantaneous frequency is the time derivative of phase function
(11)
ω
=
d
θ
d
t
.


### 2.3. Teager Energy Operator

Teager energy operator (TEO) is a powerful technique for time-frequency analysis of nonstationary signals. TEO tracks the modulation energy and determines the instantaneous frequency (FM) and instantaneous amplitude (AM) of a signal [20].For a continuous TEO,ψ
(
·
) is defined as
(12)
ψ
[
x
(
t
)
]
=
[
x
˙
(
t
)
]
2
-
x
(
t
)
x
¨
(
t
)
,
where x
˙
(
t
) and x
¨
(
t
) are the first and second time derivatives of x
(
t
), respectively. In the discrete case, the time derivatives are substituted by time differences. Thus, the discrete time TEO is given by
(13)
ψ
[
x
(
n
)
]
=
x
2
(
n
)
-
x
(
n
+
1
)
·
x
(
n
-
1
)
.According to (13), TEO only requires three samples to calculate energy. Therefore, it provides high resolution and is able to capture energy variations.Maragos et al. [20] proposed discrete energy separation algorithm (DESA-1) for the instantaneous frequency and amplitude of an AM-FM signal as follows:
(14)
y
(
n
)
=
x
(
n
)
-
x
(
n
-
1
)
,
f
(
n
)
=
arcos
(
1
-
ψ
[
y
(
n
)
]
+
ψ
[
y
(
n
+
1
)
]
4
ψ
[
x
(
n
)
]
)
|
a
(
n
)
|
=
ψ
[
x
(
n
)
]
1
-
(
1
-
(
(
ψ
[
y
(
n
)
]
+
ψ
[
y
(
n
+
1
)
]
)
/
4
ψ
[
x
(
n
)
]
)
)
2
.Or(15)
|
a
(
n
)
|
=
ψ
[
x
(
n
)
]
1
-
(
1
-
(
(
ψ
[
x
(
n
)
]
-
ψ
[
x
(
n
-
1
)
]
)
/
2
ψ
[
x
(
n
)
]
)
)
2
f
(
n
)
=
arcos
(
1
-
ψ
[
x
(
n
)
]
-
ψ
[
x
(
n
-
1
)
]
2
ψ
[
x
(
n
)
]
)
.The main drawback of TEO is its sensitivity. If the value of data points changes abruptly the discrete time differentiation yields unrealistic results. Thus, TEO is appropriate for bandwidth signals.


### 2.4. Higher-Order Energy Operator

Higher-order energy operator (HEO) is an extension of TEO to a large class of operators [21]. These generalized operators can improve the sampling rate precision and computation effort and show finer frequency resolution in the presence of noise. HEO, ψ
H
d
k
,
p
,
m
(
·
), implements four partial pth, qth, mth, and lth derivatives of a signal as follows:
(16)
ψ
H
d
k
,
p
,
m
[
x
(
t
)
]
=
x
(
p
)
(
t
)
x
(
q
)
(
t
)
-
x
(
m
)
(
t
)
x
(
l
)
(
t
)
,in which conditions ofk
=
p
+
q
=
m
+
l and (
p
,
q
)
≠
(
m
,
l
) must be satisfied. The combination of (
p
=
1, q
=
k
-
1, m
=
0, l
=
k
) yields TEO.In discrete case, continuous time derivatives are obtained by substituting two sample symmetric differences(17)
x
(
p
)
(
n
T
s
)
=
x
(
p
-
1
)
[
(
n
+
1
)
T
s
]
-
x
(
p
-
1
)
[
(
n
-
1
)
T
s
]
2
T
s
.Thus, the symmetric higher-order discrete operator is obtained byk
even:(18)
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
=
1
2
[
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
ii
-
(
x
(
n
+
m
)
x
(
n
+
l
)
+
x
(
n
-
m
)
x
(
n
-
l
)
)
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
]
,
k odd:
(19)
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
=
1
2
[
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
i
-
(
x
(
n
+
m
)
x
(
n
+
l
)
+
x
(
n
-
m
)
x
(
n
-
l
)
)
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
]
.By including higher-order derivatives, more sampling points are included which can lead to less sensitivity of energy operator against frequency fluctuations but less instantaneous results:(20)
x
1
(
n
)
=
x
(
n
)
-
x
(
n
-
1
)
,
Ω
^
2
=
ψ
H
d
k
,
p
,
m
[
x
1
(
n
)
]
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
,
A
^
=
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
sin
[
(
m
-
p
)
Ω
]
sin
[
(
m
+
p
-
k
)
Ω
]
.


### 2.5. Direct Quadrature Method

Direct quadrature (DQ) method is defined based on a normalization scheme. An iterative normalization is applied on the original signal to separate AM and FM parts [22]. First, local maxima of absolute values of data, x
(
t
), are determined. Next, a cubic spline is passed from the local maxima points to create envelope of the signal, e
1
(
t
). The first normalized data is obtained by
(21)
y
1
(
t
)
=
x
(
t
)
e
1
(
t
)
.The normalized data is considered as the original data and the same procedure is used repeatedly to obtain thenth normalized data as
(22)
y
n
(
t
)
=
y
n
-
1
(
t
)
e
n
(
t
)
=
F
(
t
)
.The normalization process stops when all the values ofy
n
(
t
) are equal to or less than unity. y
n
(
t
) is now considered as the empirical FM part of the data
(23)
y
n
(
t
)
=
cos
⁡
ϕ
(
t
)
.The empirical AM part is defined by(24)
A
(
t
)
=
x
(
t
)
F
(
t
)
=
e
1
(
t
)
×
e
2
(
t
)
×
⋯
×
e
n
(
t
)
.Also, the original signal can be reconstructed by(25)
x
(
t
)
=
A
(
t
)
×
F
(
t
)
=
A
(
t
)
cos
⁡
ϕ
(
t
)
.


## 2.1. Empirical Mode Decomposition

In 1996, Huang et al. [18] proposed empirical mode decomposition (EMD) for nonlinear and nonstationary data. EMD is an adaptive procedure which relies on two basic assumptions: (
1
) the number of extrema and zero-crossings must be equal or differ at most by one for all the data; (
2
) mean value of envelope defined by local maxima and minima must be zero at each point. Extraction of intrinsic mode functions (IMFs) is performed via sifting process described below.First, determine local extrema of the data. Plot an upper envelope by using a cubic spline passing all the extrema,e
u
(
t
). Repeat the same procedure for minima to find e
l
(
t
). Calculate their mean, m
1
(
t
), as
(1)
m
1
(
t
)
=
e
u
(
t
)
+
e
l
(
t
)
2
.The difference between signal andm
1
(
t
) is designated as the first component, h
1
(
t
):
(2)
h
1
(
t
)
=
x
(
t
)
-
m
1
(
t
)
.Now, the first component is considered as the signal and the above procedure is repeatedk times to satisfy the basic conditions:
(3)
h
1
k
(
t
)
=
h
1
(
k
-
1
)
(
t
)
-
m
1
k
(
t
)
.Thus,h
1
k
(
t
) will be the first IMF:
(4)
c
=
h
1
k
(
t
)
.A stopping criterion is computed at each repeat as follows:(5)
SD
=
∑
t
=
0
T
[
|
h
(
k
-
1
)
(
t
)
-
h
(
k
)
(
t
)
|
2
h
2
(
k
-
1
)
(
t
)
]
.When the stopping criterion is satisfied the first IMF,c
1, is obtained. To obtain other IMFs the residue, r
1, is considered as the signal and the procedure is repeated:
(6)
r
1
=
x
(
t
)
-
c
1
.When the residue is constant or a monotonic signal or no more sifting is possible the procedure is stopped. Finally, we can reconstruct the original data by summing up(7)
x
(
t
)
=
∑
j
=
1
n
c
j
+
r
n
.


## 2.2. Hilbert-Huang Transform

Hilbert-Huang transform (HHT) implements the benefits of EMD and Hilbert spectral analysis (HSA) to extract a signal to its intrinsic mode functions [19]. Hilbert transform of a signal, x
(
t
), is an integral transform that can be thought of as the convolution with the Cauchy kernel
(8)
H
[
x
(
t
)
]
=
1
π
PV
∫
-
∞
+
∞
x
(
t
)
t
-
τ
d
τ
,
where PV is the Cauchy principal value. Frequency content of a signal, x
(
t
), can be obtained by providing the conjugate part to give an analytic complex signal as
(9)
z
(
t
)
=
x
(
t
)
+
i
y
(
t
)
=
a
(
t
)
e
i
θ
(
t
)
.Thus(10)
a
(
t
)
=
x
2
(
t
)
+
y
2
(
t
)
θ
(
t
)
=
arctan
(
y
x
)
,
wherea(t) is the instantaneous amplitude and θ
(
t
) is the instantaneous phase. The instantaneous frequency is the time derivative of phase function
(11)
ω
=
d
θ
d
t
.


## 2.3. Teager Energy Operator

Teager energy operator (TEO) is a powerful technique for time-frequency analysis of nonstationary signals. TEO tracks the modulation energy and determines the instantaneous frequency (FM) and instantaneous amplitude (AM) of a signal [20].For a continuous TEO,ψ
(
·
) is defined as
(12)
ψ
[
x
(
t
)
]
=
[
x
˙
(
t
)
]
2
-
x
(
t
)
x
¨
(
t
)
,
where x
˙
(
t
) and x
¨
(
t
) are the first and second time derivatives of x
(
t
), respectively. In the discrete case, the time derivatives are substituted by time differences. Thus, the discrete time TEO is given by
(13)
ψ
[
x
(
n
)
]
=
x
2
(
n
)
-
x
(
n
+
1
)
·
x
(
n
-
1
)
.According to (13), TEO only requires three samples to calculate energy. Therefore, it provides high resolution and is able to capture energy variations.Maragos et al. [20] proposed discrete energy separation algorithm (DESA-1) for the instantaneous frequency and amplitude of an AM-FM signal as follows:
(14)
y
(
n
)
=
x
(
n
)
-
x
(
n
-
1
)
,
f
(
n
)
=
arcos
(
1
-
ψ
[
y
(
n
)
]
+
ψ
[
y
(
n
+
1
)
]
4
ψ
[
x
(
n
)
]
)
|
a
(
n
)
|
=
ψ
[
x
(
n
)
]
1
-
(
1
-
(
(
ψ
[
y
(
n
)
]
+
ψ
[
y
(
n
+
1
)
]
)
/
4
ψ
[
x
(
n
)
]
)
)
2
.Or(15)
|
a
(
n
)
|
=
ψ
[
x
(
n
)
]
1
-
(
1
-
(
(
ψ
[
x
(
n
)
]
-
ψ
[
x
(
n
-
1
)
]
)
/
2
ψ
[
x
(
n
)
]
)
)
2
f
(
n
)
=
arcos
(
1
-
ψ
[
x
(
n
)
]
-
ψ
[
x
(
n
-
1
)
]
2
ψ
[
x
(
n
)
]
)
.The main drawback of TEO is its sensitivity. If the value of data points changes abruptly the discrete time differentiation yields unrealistic results. Thus, TEO is appropriate for bandwidth signals.


## 2.4. Higher-Order Energy Operator

Higher-order energy operator (HEO) is an extension of TEO to a large class of operators [21]. These generalized operators can improve the sampling rate precision and computation effort and show finer frequency resolution in the presence of noise. HEO, ψ
H
d
k
,
p
,
m
(
·
), implements four partial pth, qth, mth, and lth derivatives of a signal as follows:
(16)
ψ
H
d
k
,
p
,
m
[
x
(
t
)
]
=
x
(
p
)
(
t
)
x
(
q
)
(
t
)
-
x
(
m
)
(
t
)
x
(
l
)
(
t
)
,in which conditions ofk
=
p
+
q
=
m
+
l and (
p
,
q
)
≠
(
m
,
l
) must be satisfied. The combination of (
p
=
1, q
=
k
-
1, m
=
0, l
=
k
) yields TEO.In discrete case, continuous time derivatives are obtained by substituting two sample symmetric differences(17)
x
(
p
)
(
n
T
s
)
=
x
(
p
-
1
)
[
(
n
+
1
)
T
s
]
-
x
(
p
-
1
)
[
(
n
-
1
)
T
s
]
2
T
s
.Thus, the symmetric higher-order discrete operator is obtained byk
even:(18)
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
=
1
2
[
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
ii
-
(
x
(
n
+
m
)
x
(
n
+
l
)
+
x
(
n
-
m
)
x
(
n
-
l
)
)
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
]
,
k odd:
(19)
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
=
1
2
[
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
i
-
(
x
(
n
+
m
)
x
(
n
+
l
)
+
x
(
n
-
m
)
x
(
n
-
l
)
)
(
x
(
n
+
p
)
x
(
n
+
q
)
+
x
(
n
-
p
)
x
(
n
-
q
)
)
]
.By including higher-order derivatives, more sampling points are included which can lead to less sensitivity of energy operator against frequency fluctuations but less instantaneous results:(20)
x
1
(
n
)
=
x
(
n
)
-
x
(
n
-
1
)
,
Ω
^
2
=
ψ
H
d
k
,
p
,
m
[
x
1
(
n
)
]
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
,
A
^
=
ψ
H
d
k
,
p
,
m
[
x
(
n
)
]
sin
[
(
m
-
p
)
Ω
]
sin
[
(
m
+
p
-
k
)
Ω
]
.


## 2.5. Direct Quadrature Method

Direct quadrature (DQ) method is defined based on a normalization scheme. An iterative normalization is applied on the original signal to separate AM and FM parts [22]. First, local maxima of absolute values of data, x
(
t
), are determined. Next, a cubic spline is passed from the local maxima points to create envelope of the signal, e
1
(
t
). The first normalized data is obtained by
(21)
y
1
(
t
)
=
x
(
t
)
e
1
(
t
)
.The normalized data is considered as the original data and the same procedure is used repeatedly to obtain thenth normalized data as
(22)
y
n
(
t
)
=
y
n
-
1
(
t
)
e
n
(
t
)
=
F
(
t
)
.The normalization process stops when all the values ofy
n
(
t
) are equal to or less than unity. y
n
(
t
) is now considered as the empirical FM part of the data
(23)
y
n
(
t
)
=
cos
⁡
ϕ
(
t
)
.The empirical AM part is defined by(24)
A
(
t
)
=
x
(
t
)
F
(
t
)
=
e
1
(
t
)
×
e
2
(
t
)
×
⋯
×
e
n
(
t
)
.Also, the original signal can be reconstructed by(25)
x
(
t
)
=
A
(
t
)
×
F
(
t
)
=
A
(
t
)
cos
⁡
ϕ
(
t
)
.


## 3. Modeling and Analysis

The SAC 3-story building in Los Angeles was selected for analysis. This building is designed as a standard office building as a part of SAC task 5.4.3 [17]. Figure 1(a) depicts the building in plan and elevation. Perimeter moment frames are indicated in Figure 1(b) with thick lines. Design of the structure is performed based on UBC-94 [23] and provisions of FEMA 267 [24] are also considered. Grade A36 steel is utilized in conjunction with Gr. 50 steel for beams and girders. Seismic mass of 65.5 kips-sec2/ft is assumed for stories and 70.9 kips-sec2/ft for the roof. It is assumed that the building is located on stiff soil [23].Configuration of the building. (a) Plan and elevation view and (b) perimeter moment frames.
(a)
(b)The open source finite element program, Opensees [25], is used for modeling. The building is modeled as a centerline 3D moment frame. Rigid diaphragm is defined for each story. No cyclic degradation is defined. “Steel01” material is utilized with “fiber” sections to model moment-rotation relationship of moment connections. 5% strain hardening is defined. “Displacement-based beam-column” element is used to model nonlinear behavior of elements. P-Δ effect is considered to account for second-order displacements during the analysis.To improve the robustness of the procedure several sensors are installed on the structure. Four sensors are mounted on the corners of each story for each direction. Figure2 illustrates the orientation of principal directions (i.e., x and y) of the building and position of the sensors.Figure 2
Location of the sensors on the building.Two horizontal components of “EL Centro Array number 1” record of “Imperial Valley” earthquake are utilized for nonlinear dynamic analyses. Both components are simultaneously imposed on the structure in two orthogonal directions for analysis.


## 4. Selection of Signal Processing Method

Time-frequency methods are robust procedures to detect instantaneous changes in frequency of a system. However, time dependence of extracted frequencies can increase complexity in identification of the structure characteristic. To clarify this issue, the frequency of vibration of the structure is investigated in healthy state. Four instantaneous time-frequency methods including two well-known HHT and TEO procedures and two newer DQ and HEO methods are implemented and the results are compared.For this reason, a nonlinear dynamic analysis is performed. A low scale factor is selected for the acceleration record to ensure elastic behavior of the structure (maximum drift ratio below 0.01). Acceleration output of sensor 17 is recorded during the nonlinear dynamic analysis (Figure3). Sampling frequency of 800 Hz is assigned for the sensor. EMD is used to extract the acceleration records of the sensor to its IMFs, as depicted in Figure 4. Only the first IMF is used for further process since it contains most of the frequency content. Finally, the aforementioned time-frequency methods are applied to extract instantaneous frequencies and amplitudes (Figure 5). Statistical variables of the frequencies are listed in Table 1. All the calculations are performed using MATLAB software package [26].Table 1
Statistical variables of extracted frequencies.
HHT
TEO
HEO
DQ
Mean
3.764
3.417
1.992
2.804
Median
2.846
2.323
1.393
0.738
Std.
15.571
13.963
3.346
89.918
RMS
16.019
14.375
3.894
89.960
Range
527.4
1942.3
340.2
8759.7
K-factor
8448.6
27920.8
1324.9
731797.6Components of ground motion acceleration and response acceleration of sensors 17 and 18.
(a)
(b)
(c)
(d)Extracted IMFs of the acceleration response.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Extracted time-frequencies of acceleration response. (a) HHT, (b) TEO, (c) HEO, and (d) DQ.
(a)
(b)
(c)
(d)As can be easily seen from Figure5, large fluctuation in frequencies is observed which is associated with numerous spikes. An outlier analysis is carried out to further investigate the outlier frequencies.


## 5. Outlier Analysis

An outlier is an unusual point that is markedly inconsistent from rest of the data. Outliers can appear accidentally in a set of data, but in most cases this deviated data is due to measurement error or natural deviations in the dataset [27]. A common practice for the former case is to discard outliers. However, statistical tools should be cautiously implemented for the latter case. Outliers can influence the estimation of statistical parameters of the observation. Therefore, it is important to distinguish the origin of outliers.


### 5.1. Univariate Analysis

This section describes the results of a univariate analysis of frequencies. The following statistical features were considered to compare the procedures including root mean square (RMS), range. andK
-factor:
(26)
RMS
=
1
n
∑
i
=
1
n
x
i
2
,
Range
=
max
⁡
(
x
i
)
-
min
⁡
(
x
i
)
,
K
-factor
=
max
⁡
(
x
i
)
×
RMS,
where x
i is data sample. Table 1 compares the aforementioned extracted frequencies.Comparison of mean and median values shows the deviation of mean values due to presence of outliers that are significantly deviant from the rest of measurements. Therefore, median values are more appropriate for comparison.HHT yields the largest median frequency equal to 2.84. Standard deviation of the HHT is large due to high values of outliers. Although the range of frequencies for HHT is less than TEO and DQ, the number of outliers with large values is high which leads to large RMS andk-factor.In contrast to HHT, HEO has low median frequency and the lowest range, RMS, andK-factor. The extremely low RMS andk-factors show that although the number of outliers is not less than the others, the outliers with very large values are less. It is notable that although HEO shows poor accuracy to predict the median frequency in comparison to TEO, the low standard deviation of the results shows that using more data samples for numerical differentiation (17) can reduce the fluctuations of the results. However, this can adversely affect the instantaneous nature of the HEO and leads to poor estimations. DQ has the lowest median frequency among all. Also, DQ has the highest range, RMS, and K
-factor. The combination of numerous outliers with considerable deviation from the rest of data set has resulted in large difference between median and mean values and also large standard deviation. This can be due to normalization process (22) which may disassemble the frequency content.


### 5.2. Multivariate Analysis

To study the effect of outliers for all the sensors simultaneously, a multivariate outlier analysis is performed using a discordancy test. Mahalanobis squared distance (MSD) is designated as a discordancy test variant. MSD for ap-dimensional problem withn samples is defined as
(27)
D
ξ
=
(
{
x
ξ
}
-
{
x
-
}
)
T
·
[
K
]
-
1
·
(
{
x
ξ
}
-
{
x
-
}
)
,
where {
x
-
} is the mean vector of data samples, {
x
ξ
} is the observed outlier, and [
K
] is the covariance matrix. {
x
-
} and {
x
ξ
} arep-dimensional vectors, and [
K
] is a p
×
p square matrix. Figure 6 plots MSD values for the assumed signal processing procedures.MSD values for frequencies of all sensors. (a) HHT, (b) TEO, (c) HEO, and (d) DQ.
(a)
(b)
(c)
(d)It is clear that HHT has the lowest values of MSD. Also the number of MSDs with large values is the least among all. The most outspread results belong to TEO. HEO shows a slight improvement in comparison to TEO. However, the number of MSDs with large values is still remarkable. Although DQ has lower values of MSD compared to TEO and HEO, large discordance is observed specially for the first 5000 samples (first 6 second of seismic analysis). According to the above results, HHT is selected for further postprocessing since it provides more appropriate results.Since the outliers are observed in the undamaged case of the structure, they cannot be considered as damage spikes as a result of sudden change in structural stiffness. Therefore, we assumed that these outliers cannot be useful to detect damage and eliminated them since they adversely affect the postprocessing results.


### 5.3. Sampling Rate Effect

A sensitivity analysis is carried out to study the sampling rate effect for HHT procedure. The previous statistical variables are used to measure dispersion of data. Table2 refers to statistical characteristics of extracted frequency for sampling rate of 25 to 800 Hz.Table 2
Statistical variables for different sampling frequencies.
Sampling freq. (Hz)
25
50
100
200
400
800
Mean
3.278
3.206
3.614
3.646
3.764
3.764
Median
2.840
2.827
2.837
2.835
2.844
2.846
Std.
12.211
12.052
14.823
15.267
15.568
15.571
RMS
12.608
12.442
15.253
15.698
16.016
16.019
Range
523.7
523.7
527.4
527.4
527.4
527.4
K-factor
6603.1
6516.2
8044.4
8279.5
8446.8
8448.6As shown in Table2, the median frequencies are almost constant. But standard deviation increases by increase in sampling rate which affects the mean frequency. Increase in dispersion due to sampling rate increase is also confirmed by RMS, range, andk-factor which is due to presence of spikes (outliers) in frequencies. These spikes can appear as a result of numerical differentiation of phase function which is sensitive to sampling time step (11). However, difference in results of 400 Hz and 800 Hz is not sensible. Generally speaking, frequency results are less sensitive to sampling rates over 400 Hz.On the other hand, the lower values of spikes for lower sampling frequencies (i.e., 25 and 50 Hz) are due to lowerNyquist frequency that limits the maximum frequency representable in the data sample, not as a result of more precise extraction. Thus, 800 Hz sampling frequency was designated for postprocessing. However, the outliers were eliminated as mentioned before.


## 5.1. Univariate Analysis

This section describes the results of a univariate analysis of frequencies. The following statistical features were considered to compare the procedures including root mean square (RMS), range. andK
-factor:
(26)
RMS
=
1
n
∑
i
=
1
n
x
i
2
,
Range
=
max
⁡
(
x
i
)
-
min
⁡
(
x
i
)
,
K
-factor
=
max
⁡
(
x
i
)
×
RMS,
where x
i is data sample. Table 1 compares the aforementioned extracted frequencies.Comparison of mean and median values shows the deviation of mean values due to presence of outliers that are significantly deviant from the rest of measurements. Therefore, median values are more appropriate for comparison.HHT yields the largest median frequency equal to 2.84. Standard deviation of the HHT is large due to high values of outliers. Although the range of frequencies for HHT is less than TEO and DQ, the number of outliers with large values is high which leads to large RMS andk-factor.In contrast to HHT, HEO has low median frequency and the lowest range, RMS, andK-factor. The extremely low RMS andk-factors show that although the number of outliers is not less than the others, the outliers with very large values are less. It is notable that although HEO shows poor accuracy to predict the median frequency in comparison to TEO, the low standard deviation of the results shows that using more data samples for numerical differentiation (17) can reduce the fluctuations of the results. However, this can adversely affect the instantaneous nature of the HEO and leads to poor estimations. DQ has the lowest median frequency among all. Also, DQ has the highest range, RMS, and K
-factor. The combination of numerous outliers with considerable deviation from the rest of data set has resulted in large difference between median and mean values and also large standard deviation. This can be due to normalization process (22) which may disassemble the frequency content.


## 5.2. Multivariate Analysis

To study the effect of outliers for all the sensors simultaneously, a multivariate outlier analysis is performed using a discordancy test. Mahalanobis squared distance (MSD) is designated as a discordancy test variant. MSD for ap-dimensional problem withn samples is defined as
(27)
D
ξ
=
(
{
x
ξ
}
-
{
x
-
}
)
T
·
[
K
]
-
1
·
(
{
x
ξ
}
-
{
x
-
}
)
,
where {
x
-
} is the mean vector of data samples, {
x
ξ
} is the observed outlier, and [
K
] is the covariance matrix. {
x
-
} and {
x
ξ
} arep-dimensional vectors, and [
K
] is a p
×
p square matrix. Figure 6 plots MSD values for the assumed signal processing procedures.MSD values for frequencies of all sensors. (a) HHT, (b) TEO, (c) HEO, and (d) DQ.
(a)
(b)
(c)
(d)It is clear that HHT has the lowest values of MSD. Also the number of MSDs with large values is the least among all. The most outspread results belong to TEO. HEO shows a slight improvement in comparison to TEO. However, the number of MSDs with large values is still remarkable. Although DQ has lower values of MSD compared to TEO and HEO, large discordance is observed specially for the first 5000 samples (first 6 second of seismic analysis). According to the above results, HHT is selected for further postprocessing since it provides more appropriate results.Since the outliers are observed in the undamaged case of the structure, they cannot be considered as damage spikes as a result of sudden change in structural stiffness. Therefore, we assumed that these outliers cannot be useful to detect damage and eliminated them since they adversely affect the postprocessing results.


## 5.3. Sampling Rate Effect

A sensitivity analysis is carried out to study the sampling rate effect for HHT procedure. The previous statistical variables are used to measure dispersion of data. Table2 refers to statistical characteristics of extracted frequency for sampling rate of 25 to 800 Hz.Table 2
Statistical variables for different sampling frequencies.
Sampling freq. (Hz)
25
50
100
200
400
800
Mean
3.278
3.206
3.614
3.646
3.764
3.764
Median
2.840
2.827
2.837
2.835
2.844
2.846
Std.
12.211
12.052
14.823
15.267
15.568
15.571
RMS
12.608
12.442
15.253
15.698
16.016
16.019
Range
523.7
523.7
527.4
527.4
527.4
527.4
K-factor
6603.1
6516.2
8044.4
8279.5
8446.8
8448.6As shown in Table2, the median frequencies are almost constant. But standard deviation increases by increase in sampling rate which affects the mean frequency. Increase in dispersion due to sampling rate increase is also confirmed by RMS, range, andk-factor which is due to presence of spikes (outliers) in frequencies. These spikes can appear as a result of numerical differentiation of phase function which is sensitive to sampling time step (11). However, difference in results of 400 Hz and 800 Hz is not sensible. Generally speaking, frequency results are less sensitive to sampling rates over 400 Hz.On the other hand, the lower values of spikes for lower sampling frequencies (i.e., 25 and 50 Hz) are due to lowerNyquist frequency that limits the maximum frequency representable in the data sample, not as a result of more precise extraction. Thus, 800 Hz sampling frequency was designated for postprocessing. However, the outliers were eliminated as mentioned before.


## 6. Damage Detection Methodology

This section investigates the capability of HHT to detect damage in a nonlinear moment frame. First, the moment connections of the structures are modeled to present nonlinear behavior. Second, they are modeled to behave elastically. An IDA-based [28] approach is utilized for analysis. Seven intensity levels are defined for analysis covering the range of elastic (maximum drift ratio of 0.01) to near collapse of the structure (maximum drift ratio of 0.1) in accordance with FEMA350 [29]. Table 3 lists the details of the intensity levels and the corresponding response of the nonlinear structure.Table 3
Response of the structure for different intensity levels.
Intensity level
Level 1
Level 2
Level 3
Level 4
Level 5
Level 6
Level 7
Max driftx-dir.
0.0007
0.0066
0.0110
0.0201
0.0315
0.0546
0.0844
Max drifty-dir.
0.0016
0.0098
0.0182
0.0356
0.0485
0.0658
0.0891
Max disp.x-dir.
0.0065
0.0342
0.0669
0.1496
0.2181
0.4473
0.7350
Max disp.y-dir.
0.0116
0.0501
0.0741
0.2009
0.4241
0.6866
0.8712It is worth noting that no damage occurs in intensity level 1. Response acceleration of the structure under the aforementioned intensity levels is recorded by all the sensors for both elastic and nonlinear cases.Three features are implemented to detect damage including frequency, amplitude, and energy. Frequency and amplitude are directly extracted from the HHT. Energy index is defined as sum of square of the first IMF [30]:
(28)
E
=
∫
0
t
0
(
IMF
)
2
d
t
,
wheret
0 is the time length of the signal.We did not filter the frequencies since filtering may disassemble the frequency content of the response. But the outlier frequencies are discarded, as aforementioned. According to Figure7, Fourier spectrum of response of the healthy state shows that vibration frequency of the structure is mainly between 1 and 7 Hz. Thus, the frequencies outside this range are cut off.Figure 7
Fourier spectrum of acceleration response.To quantify changes in the frequencies and amplitudes, fuzzyc-means clustering [31] is employed. Fuzzyc-means allows each of the data to belong to several clusters by using rationale of fuzzy sets. The aim is to minimize the least-squared error of data,x
i, and center of cluster,c
j:
(29)
J
m
=
∑
i
=
1
N
∑
j
=
1
C
u
i
j
m
∥
x
i
-
c
j
∥
,
where m is weighting exponent (
1
≤
m
<
∞
) and u
i
j is degree of membership of x
i for cluster j. Partitioning is an iterative process. In each iteration, the memberships and center of clusters are updated as follows:
(30)
u
i
j
=
1
∑
k
=
1
C
(
∥
x
i
-
c
j
∥
/
∥
x
i
-
c
k
∥
)
2
/
(
m
-
1
)
,
c
j
=
∑
i
=
1
N
u
i
j
m
·
x
i
∑
i
=
1
N
u
i
j
m
.When the iteration thresholdma
x
i
j
{
|
u
i
j
(
k
+
1
)
-
u
i
j
(
k
)
|
}
<
ε is reached, the iterations are terminated.Before clustering data, a silhouette statistic [32] is employed to find the appropriate number of clusters. For a given sampling pointi, d
(
i
,
c
) is the average dissimilarity of i to all other data points for any cluster c. If b
i is defined by the minimum of all d
(
i
,
c
), the silhouette width for theith point is
(31)
s
w
i
=
(
b
i
-
a
i
)
max
⁡
(
b
i
,
a
i
)
.By finding the average of silhouette widths for all sampling points the average silhouette width is obtained by(32)
s
w
ave
=
1
n
∑
i
=
1
n
s
w
i
.Larger silhouette width indicates better clustered data. For this case using two clusters yields the best results. For instance, the average silhouette width of frequencies for different numbers of clusters is plotted in Figure8.Figure 8
The average silhouette widths.In this study, the value of center of the first cluster for frequencies and amplitudes is designated as frequency and amplitude feature indexes, respectively. Also, the energy value is considered as energy feature index. Figure9 shows the major components of the procedure.Figure 9
Overview of the algorithm.To detect damage, the sum of residues between nonlinear and elastic feature indexes for all the sensors is considered as damage index. Table4 lists the damage indexes for different intensity levels. As shown in the table, damage indexes of all features are zero for intensity level 1 implying the undamaged state for the nonlinear structure. Frequency index rapidly increases for intensity levels 2 and 3 since the structure does not tolerate severe nonlinearity. The frequency index remains approximately constant for higher intensity levels. It can be due to hysteresis behavior of moment connections which results in successive changes in stiffness of members and consequently successive changes in vibration frequency of the structure. When large numbers of beam-column connections undergo successive frequency changes under a moderate or severe seismic action (e.g., intensity levels 4 to 7) the overall frequency change of the structure will be complicated. In this study no drastic change in frequency index was observed for high-intensity levels. Thus, frequency may not be an appropriate damage index to estimate extent of damage. Amplitude index shows more uniform behavior. No residue is observed for intensity level 1. However, intensity level 2 shows a nonnegligible residue since the structure behaves slightly nonlinearly. The rate of increase in amplitude index increases for higher levels of damage. Energy index shows the most rational results. The energy index increases more rapidly after 0.01 maximum drift.Table 4
Feature indexes for intensity levels.
Damage
Level 1
Level 2
Level 3
Level 4
Level 5
Level 6
Level 7
x-direction
Frequency
0
1.940
4.193
5.643
6.019
6.292
6.189
Amplitude
0
4.982
16.044
44.701
75.828
109.658
143.224
Energy
0
1.676
e
+
06
8.642
e
+
06
3.957
e
+
07
9.317
e
+
07
1.687
e
+
08
2.650
e
+
08
y-direction
Frequency
0
2.904
5.648
8.658
9.434
9.434
9.468
Amplitude
0
6.550
21.424
60.805
104.836
154.087
202.991
Energy
0
2.182
e
+
06
1.151
e
+
07
5.548
e
+
07
1.385
e
+
08
2.683
e
+
08
4.481
e
+
08Figure10(a) shows the behavior of normalized values of the indexes versus maximum drift ratios for x
- and y-directions. It can be easily seen that amplitude and energy indexes show better correlation with maximum drift ratios than frequencies.Correlation of feature indexes with (a) maximum drift ratio and (b) spectral acceleration.
(a)
(b)Also, feature indexes are plotted versus spectral acceleration (S
a) as an intensity measure in Figure 10(b) to study the ability of feature indexes to be implemented as a damage measure (DM) [28] in IDA analysis. Acceptable slope change of the amplitude and energy indexes is observed for both directions according to the figure. In contrast, frequency curve shows a negative slope at higher intensity levels. This ability is investigated in more detail in the next section.


## 7. Ground Motion Variability

The procedure in the previous section was tested for only one ground motion record. However, large variability is often observed in characteristics of different ground motions. In attempt to investigate the damage of the structure under different seismic actions the above procedure was assessed for a set of ground motion records. FEMA-P695 [33] far-field record set was selected for analysis. This record set includes twenty-two pairs of horizontal components of ground motions that is claimed to be adequate to consider record-to-record variability. An IDA analysis is carried out. The analysis is stopped when 10% maximum story drift ratio is reached for each ground motion record. Figure 11 refers to IDA curves of the features in x-direction. The measures are normalized to provide an ability to compare the features.IDA curves for feature indexes. (a) Frequency, (b) amplitude, and (c) energy.
(a)
(b)
(c)Slope of frequency IDA curves significantly increases for high level of spectral acceleration implying the lower rate of frequency change for these intensity levels. Nevertheless, IDA curves show relatively erratic behavior in some cases. Abrupt slope changes occur even in low damage measures. Regarding the figure, amplitudes present more smooth behavior. Slope of amplitude IDA curves gradually decreases to reach the maximum damage measure. However, no sensible slope change is observed after amplitude index of about 0.2, on average, making estimation of extent of damage difficult. Slope of IDA curves for energy index shows a smooth transition from elastic to collapse state. In contrast to amplitude feature, one can observe a moderate slope reduction for energy feature over 0.2. These results are consistent with the previous results.As expected, large record-to-record variability exists in IDA curves of all the feathers. For interpretation of randomness included in the records, a probabilistic approach is employed. IDA curves are summarized to form 16%, 50%, and 84% cross-sectional fractiles [28], as depicted in Figure 12.IDA curves for feature indexes. (a) Frequency, (b) amplitude, and (c) energy.
(a)
(b)
(c)Although the IDA curves of frequency behave irregularly for some records, the median curve shows more uniform behavior. However, it is notable that the median curve cannot show the frequency tendency properly for normalized values of more than about 0.9 since the frequency indexes stay constant or decrease. Still, energy seems to be a more reliable feature to evaluate extent of damage.To quantify the damage extent, 5 limit states are defined. Table5 compares median values of DM in x
- and y-directions. All three median values follow almost similar paths in both directions. This is more obvious for amplitude and energy index. Frequency median values in y-direction show more rapid increase inS
a from 0.5 to 1.0 which may be a result of dispersion of IDA curves in this range since the other features do not show this deviation in y-direction.Table 5
Median values of damage measures.
x-direction
y-direction
DM
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
Frequency
0.56
1.27
1.77
2.62
5.71
0.63
1.30
2.15
3.62
5.91
Amplitude
2.05
3.15
4.18
5.20
6.14
2.00
3.16
4.21
5.10
5.97
Energy
2.72
3.80
4.74
5.46
6.14
2.74
3.87
4.69
5.37
5.97To evaluate sufficiency of the damage measures, dispersion of estimations for different limit states is investigated. It is clear that more dispersion in the results reduces the confidence in estimation of the features. Dispersion of the damage measures is estimated using the conditional factional standard deviation defined by half the distance between 16% and 84% values for a certain damage measure in logarithmic scale [34]. Table 6 shows that frequency presents the largest dispersion. Amplitude and energy indexes show less dispersion which is almost constant for different limit states.Table 6
Dispersion of intensity level for different damage measures.
x-direction
y-direction
DM
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
Frequency
0.56
0.61
0.64
0.45
0.45
0.57
0.66
0.55
0.52
0.52
Amplitude
0.36
0.35
0.37
0.37
0.37
0.36
0.34
0.33
0.30
0.30
Energy
0.36
0.38
0.38
0.38
0.38
0.30
0.30
0.31
0.30
0.30


## 8. Conclusion

This study aims to find a promising method to detect damage in nonlinear moment frames. The capability of several time-frequency procedures is evaluated to find a suitable signal processing method. Then a probabilistic approach is employed to find the extension of damage under various seismic actions. Based on the results, the following major conclusions can be drawn.(i)
Despite the fact that several time-frequency procedures are proposed recently, HHT is still a robust procedure to extract frequencies of a signal.
(ii)
Large dispersion was observed in estimation of damage features due to large record-to-record variability. Thus, a probabilistic approach can be useful to quantify damage.
(iii)
Energy variable yields more appropriate results than frequency and amplitude for damage detection of nonlinear moment frames.



---

*Source: 523675-2014-03-13.xml*