# Structural Damage Detection by Power Spectral Density Estimation Using Output-Only Measurement

Authors: Eun-Taik Lee; Hee-Chang Eun
Journal: Shock and Vibration
 (2016)

---

## Abstract

Most damage detection methods have difficulty in detecting damage using only measurement data due to the existence of external noise. It is necessary to reduce the noise effect to obtain accurate information and to detect damage by the output-only measurement without baseline data at intact state and input data. This work imported the power spectral density estimation (PSE) of a signal to reduce the noise effect. By estimating the PSE to characterize the frequency content of the signal, this study proposes a damage detection method to trace the damage by the curvature of the PSE. Two numerical applications examine the applicability of the proposed method depending on a window function, frequency resolution, and the number of overlapping data in the PSE method. The knowledge obtained from the numerical applications leads to a series of experiments that substantiate the potential of the proposed method.


---

## Body

## 1. Introduction

Early detection of defects in civil structures is a critical process in assisting structural maintenance and management plans. With a robust damage detection methodology, it becomes possible to fix the structure during the early stages of damage.The existing damage detection methods are performed using measurement data in the time-domain or frequency-domain. The time-domain methods usually go through stochastic processes. Cattarius and Inman [1] provided a nondestructive time-domain approach to examining structural damage using time histories of the vibration response of the structure. Majumder and Manohar [2] developed a time-domain approach to detecting damage in bridge structures by analyzing the combined system of the bridge and vehicle. Sandesh and Shankar [3] presented a time-domain damage detection scheme based on a substructure system identification method using Genetic Algorithms and Particle Swarm Optimization to filter out the updated parameters. Lu and Gao [4] presented a time series model for the diagnosis of structural damage considering a damage sensitive feature without input excitation.Many damage detection methods require information about the baseline data and the input data. However, the data collection is not always practical because they cannot be readily obtained. Output-only methods use only the vibration response signals and may be classified into nonparametric methods based on corresponding time series representations and parametric methods based on scalar or vector parametric time series representations. Yao and Pakzad [5] presented an autoregressive method with exogenous input modeling for measurements and investigated an application to detect damage in a space truss structure. Power spectral density (PSD), defined as the squared value of the signal, describes the power of a signal or time series distributed over different frequencies. The PSD is the Fourier transform of the autocorrelation function, which provides the transformation from the time-domain to the frequency-domain. By estimating the PSD localized in bandwidth regions near resonance, Liberatore and Carman [6] presented a method of structural damage identification. Beskhyroun et al. [7] proposed a damage identification method to detect damage and its location based on changes in the PSD curvature before and after damage. Bayissa and Haritos [8] presented a damage identification method based on bending moment response PSD. Zhou et al. [9] provided a damage identification method based on the PSD transmissibility using output-only response measurements. From the sensitivities of PSD with respect to the structural damage parameters and finite element model updating, Chen et al. [10] presented a method to identify structural damage. Zheng et al. [11] considered a structural damage detection method from the finite element model, which is updated using the measured PSD.PSE estimates the spectral density of a random signal from a sequence of time sample. The techniques may be based on time-domain or frequency-domain analysis. Welch [12] provided a method for the application of the fast Fourier transform (FFT) algorithm for the estimation of power spectra. The Welch method is a nonparametric method, including the periodogram, with the advantage of possible implementation using FFT. The Welch method is used to find the PSD of a signal and to reduce the effect of noise. Vamvoudakis-Stefanou et al. [13] analyzed six well-known, output-only statistical time series methods and compared the methods through the detection and identification of various types of damage. Kopsaftopoulos and Fassois [14, 15] investigated the potential and effectiveness of the statistical time series methods through experiments on a laboratory aluminum truss structure. They assessed several scalar and vector statistical time series methods for vibration based on structural health monitoring. Gupta et al. [16] compared several types of window functions and observed that the rectangular and Hamming windows gave better results than the Bartlett, Hamming, and Blackman windows. Zimin and Zimmerman [17] compared Time-Domain Periodogram Analysis (TDPA) with Frequency-Domain Periodogram Analysis (FDPA) for synthetic signals and suggested that TDPA can be utilized as an index to evaluate the existence of structural damage. Giles and Spencer Jr. [18] provided a damage detection algorithm using the changes in the PSD of a structure between undamaged and damaged systems and only output the measurement data. Fang and Perera [19] introduced power mode shape curvature and power flexibility, and they proposed the damage detection method using their variation between undamaged and damaged states. Masciotta et al. [20] presented a spectrum-driven method to detect the damage using a proper combination of the eigenparameters extracted from the PSD matrix.This work uses only the output data without the baseline data, unlike the existing methods, to compare the difference in damage indices such as PSE, mode shapes, and flexibility between the undamaged and damaged states. There are damage detection methods that are very sensitive to external noise so that the damage cannot be traced. This study provides a damage detection method that utilizes the response data transformed to the frequency-domain from the time-domain and the Welch method to find the PSE of a signal. Damage is evaluated using the curvature of the PSE. This work investigates the validity of the damage detection method depending on the rectangular, Hamming, Bartlett, Hann, and Blackman windows with the Welch method, overlapping of 40% and 50%, and frequency resolution of 0.5 Hz and 1.0 Hz. Two numerical applications for detecting damage in a beam structure examine and compare the applicability of the proposed method depending on the window function, the frequency resolution, and the number of overlapping data. The knowledge obtained from the numerical applications leads to a series of experiments that substantiate the potential of the proposed method.


## 2. Damage Detection Scheme Based on the Welch Method

The Welch method divides the time series data into segments, computes a modified periodogram of each segment, and averages the PSE. A portion of the data stream near the boundaries of the window function is ignored in the analysis, and its situation can be improved by letting the segments overlap. The PSD represents the strength of the variations as a function of frequency. The spectral density characterizes the frequency content of the signal and its estimation detects any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities. The PSESxxf is the discrete Fourier transform (DFT) of the autocorrelation estimate Rxxk or(1)Sxxf=∑k=-∞∞Rxxke-iωkT,where i=-1. This assumes that N point data sequence Xn is a discrete time, random process with an autocorrelation function Rxxk. The number of overlapping data, window function, and frequency resolution affect the PSE in (1) and are defined in the following.


### 2.1. Sample Overlapping

DividingXn, 0≤n≤N-1, from a stationary, second-order stochastic sequence into L segments of M samples (N≥LM), the data segment can be expressed as(2)X1n=Xn,n=0,1,2,…,M-1.Similarly(3)X2n=Xn+D,n=0,1,2,…,M-1,Xkn=Xn+k-1D,k=0,1,2,…,L-1,where kD is the starting point for the kth sequence. If D=M, the segment does not overlap.Welch’s method to modify Bartlett’s method applies the windowWn directly to the data segment before computing the individual periodograms. The L modified or windowed periodogram to form X1nWn,…,XknWn, n=0,1,…,M-1, can be defined as(4)S^xxkf=1UM∑n=0M-1XknWne-i2πfn2,k=0,1,…,L-1,where i=-1 and U is a normalization factor for the power written as(5)U=1M∑n=0M-1W2n.The Welch PSESxxwf is the average of modified periodogram defined as(6)Sxxwf=1L∑k=0L-1S^xxkf.And its expected value is given by(7)ESxxwf=SxxfWf,where(8)Sxxf=1NXf2,Wf=1UM∑n=0M-1Wne-i2πfn.The window functions Wn considered in this study are as follows.


### 2.2. Window Functions

Rectangular window:(9)Wn=1,n<N-12,0,otherwise.Bartlett window:(10)Wn=1-2nN,n<N-12,0,otherwise.Hann window:(11)Wn=0.5+0.5cos⁡2πnN,n<N-12,0,otherwise.Hamming window:(12)Wn=0.54+0.46cos⁡2πnN,n<N-12,0,otherwise.Blackman window:(13)Wn=0.42-0.5cos⁡2πnN-1+0.08cos⁡4πnN-1,n<N-12,0,otherwise.


### 2.3. Frequency Resolution

Frequency resolution is the distance in Hz between two adjacent data points in DFT. It is defined as(14)Fr=FsN,where Fs denotes the sampling rate and N is the number of data. From the relationship of the sampling rate and time, N=Fst, it can also be expressed by(15)Fr=1t.This equation indicates that the frequency resolution is affected by the prescribed time.The normalizing factorU is required so that the PSE of the modified periodogram will be asymptotically unbiased. The PSE is determined by the spectral resolution of each segment of length L, and it depends on the window function to minimize the effect of leakage to better represent the frequency spectrum of the data.This work considers the optimal PSE data corresponding to the frequency to display the maximum energy. The curvature is utilized as an index to evaluate the damage because the beam is characterized by the flexural response. The damage present in the region represents the abrupt variation in the curvature. The curvature at each locationl, Sxx,lw′′, is numerically obtained by a central difference approximation:(16)Sxx,lw′′=Sxx,l-1w-2Sxx,lw+Sxx,l+1wh2,l=1,2,…,N-2,where h is the distance between two successive nodes. The validity of the proposed method is illustrated by a numerical application and a beam test.


## 2.1. Sample Overlapping

DividingXn, 0≤n≤N-1, from a stationary, second-order stochastic sequence into L segments of M samples (N≥LM), the data segment can be expressed as(2)X1n=Xn,n=0,1,2,…,M-1.Similarly(3)X2n=Xn+D,n=0,1,2,…,M-1,Xkn=Xn+k-1D,k=0,1,2,…,L-1,where kD is the starting point for the kth sequence. If D=M, the segment does not overlap.Welch’s method to modify Bartlett’s method applies the windowWn directly to the data segment before computing the individual periodograms. The L modified or windowed periodogram to form X1nWn,…,XknWn, n=0,1,…,M-1, can be defined as(4)S^xxkf=1UM∑n=0M-1XknWne-i2πfn2,k=0,1,…,L-1,where i=-1 and U is a normalization factor for the power written as(5)U=1M∑n=0M-1W2n.The Welch PSESxxwf is the average of modified periodogram defined as(6)Sxxwf=1L∑k=0L-1S^xxkf.And its expected value is given by(7)ESxxwf=SxxfWf,where(8)Sxxf=1NXf2,Wf=1UM∑n=0M-1Wne-i2πfn.The window functions Wn considered in this study are as follows.


## 2.2. Window Functions

Rectangular window:(9)Wn=1,n<N-12,0,otherwise.Bartlett window:(10)Wn=1-2nN,n<N-12,0,otherwise.Hann window:(11)Wn=0.5+0.5cos⁡2πnN,n<N-12,0,otherwise.Hamming window:(12)Wn=0.54+0.46cos⁡2πnN,n<N-12,0,otherwise.Blackman window:(13)Wn=0.42-0.5cos⁡2πnN-1+0.08cos⁡4πnN-1,n<N-12,0,otherwise.


## 2.3. Frequency Resolution

Frequency resolution is the distance in Hz between two adjacent data points in DFT. It is defined as(14)Fr=FsN,where Fs denotes the sampling rate and N is the number of data. From the relationship of the sampling rate and time, N=Fst, it can also be expressed by(15)Fr=1t.This equation indicates that the frequency resolution is affected by the prescribed time.The normalizing factorU is required so that the PSE of the modified periodogram will be asymptotically unbiased. The PSE is determined by the spectral resolution of each segment of length L, and it depends on the window function to minimize the effect of leakage to better represent the frequency spectrum of the data.This work considers the optimal PSE data corresponding to the frequency to display the maximum energy. The curvature is utilized as an index to evaluate the damage because the beam is characterized by the flexural response. The damage present in the region represents the abrupt variation in the curvature. The curvature at each locationl, Sxx,lw′′, is numerically obtained by a central difference approximation:(16)Sxx,lw′′=Sxx,l-1w-2Sxx,lw+Sxx,l+1wh2,l=1,2,…,N-2,where h is the distance between two successive nodes. The validity of the proposed method is illustrated by a numerical application and a beam test.


## 3. Numerical Applications

A numerical application was performed to detect damage of a finite element model of a fixed-end beam, as shown in Figure1. The nodal points and the elements are numbered as shown in the figure. A beam with a length of 1 m is modeled using 50 beam elements. The beam has an elastic modulus of 2.0×105 MPa and a unit mass of 7,860 kg/m3. The beam’s gross cross section is 75 mm × 9 mm, and its damage section is established as 75 mm × 8.5 mm. The damping matrix is assumed as Rayleigh damping and is expressed by the stiffness matrix and a proportional constant of 0.0001. This application considers a beam with multiple damage at elements 19 and 41.Figure 1
A fixed-end beam structure model.Assume the beam structure is subjected to the white noise excitation at the left support with a step size of 0.001 sec and 0.002 sec for one second, as shown in Figure2. The dynamic responses in the time-domain are simulated by solving the dynamic equation of the finite element model subject to the excitation, and the response data of the same length as the input force are collected. We assume that the information about the excitation is unknown for this numerical application.Figure 2
Acceleration of external excitation at the left support: (a) time step of 0.001 sec, (b) time step of 0.002 sec.
(a)
(b)The dynamic responses of the practical system are contaminated by external noise. In this study, the measured response data are simulated by adding a series of random numbers to the calculated displacement responses. The displacementsur for describing the practical responses at node r can be calculated from the simulated noise-free displacements, ur0, by the following equation:(17)ur=ur01+ασr,where α denotes the relative magnitude of the error and σj is a random number variant in the range -1,1. The measured displacement responses, including the contaminated external noise, are established by inserting α=3% into (17).Most damage detection methods are sensitive to external noise contained in the measurement data. The measurement data sets collected from repeated numerical simulations are used to reduce the effect of the external noise. For this numerical application, ten data sets were utilized. Taking the FFT, the ten displacement response data sets in the time-domain were transformed into responses in the frequency-domain. Figure3(a) represents the absolute amplitude curves of the displacement responses in the frequency-domain after transforming 1,000 displacement response data points at all nodes for one second with a time step of 0.001 seconds. The first resonance frequency is located at approximately 36 Hz/sample. For the maximum values corresponding to all nodes in the neighborhood of 36 Hz/sample, the curve to connect them is exhibited in Figure 3(b). It is difficult to trace the damage location from the plot. Considering that the beam is a flexural member, the flexural curvature is predicted by a central difference method. Figure 3(c) displays the curvature curves on ten data sets, and abrupt changes in all cases are found near nodes 19 and 41, including the damage. The FFT data set can be utilized to detect damage. The concept is expanded to the periodogram and PSE to detect damage because the PSE Welch method originated from the FFT.Figure 3
Numerical results in the frequency-domain: (a) FFT curves, (b) maximum FFT curve corresponding to the first resonance frequency, and (c) its curvature curves.
(a)
(b)
(c)This work estimates the Welch PSE by dividing the response into eight segments in length. The numerical results are compared according to the following test parameters: the rectangular, Hamming, Bartlett, Hann, and Blackman window functions, sample overlaps of 40% and 50%, and frequency resolutions of 0.5 Hz and 1.0 Hz. Figures4(a) and 4(b) represent the PSE curves for 50% overlapping samples between two adjacent segments of eight segments with 1000 and 500 sampling record lengths, respectively, with a rectangular window. The signal power is concentrated at approximately 36 Hz/sample, and the curves are characterized by the frequency resolution.Figure 4
PSE curves: (a) Fr = 1.0 Hz, (b) Fr = 0.5 Hz.
(a)
(b)Figure5 represents the PSE and the curvature curves according to the window functions. The curve in Figure 5(a) consists of the maximum PSE values at all nodes at approximately 36 Hz/sample and displays the maximum energy in the PSE curves of Figure 4(a), which looks similar to the fundamental mode of the beam structure. Figures 5(b)–5(f) display the curvature curves according to the window functions on the ten data sets with the identical condition of 50% overlap and Fr = 0.5 Hz. It is shown that the damage is located in the region representing the abrupt change in the curvature. All plots exhibit more abrupt changes in the curvature at the damage location than in any other regions, which indicate that the window functions introduced in this study can be used to detect damage.Figure 5
PSE curve and its curvature under the condition of Fr = 0.5 Hz: (a) PSE curve (rectangular window, 50% overlap) by maximum PSE values around 36 Hz/sample, (b) PSE curvature (rectangular window), (c) PSE curvature (Hamming window), (d) PSE curvature (Bartlett window), (e) PSE curvature (Hann window), and (f) PSE curvature (Blackman window).
(a)
(b)
(c)
(d)
(e)
(f)Figure6 compares the PSE curvatures of ten data sets depending on the overlaps of 50% and 40% under the same Fr = 0.5 Hz using rectangular and Hamming windows. The 50% overlap leads to a larger change in the damage region than the 40% overlap. This indicates that more data overlap provides more conservative results. Figure 7 compares the applicability of damage detection by the PSE curvature depending on frequency resolutions of 0.5 Hz and 1.0 Hz using rectangular and Hamming windows. In the plots, the low frequency resolution leads to a more abrupt change in the curvature, and the damage region can be more explicitly found. This numerical application demonstrates that the damage detection method using the data samples with 50% sample overlap and 0.5 Hz frequency resolution with the window functions considered in this study is more effective.Figure 6
PSE curvature depending on the overlap: (a) rectangular window, (b) Hamming window.
(a)
(b)Figure 7
PSE curvature: (a) rectangular window, 50% overlap, (b) Hamming window, 50% overlap.
(a)
(b)As another numerical application, a crossed-beam structure with multiple damages in Figure8 is considered to detect the damage location. Fixed-end beam 1 is jointed on fixed-end beam 2. The nodal points and the members are numbered as shown in the figure. Assuming a Bernoulli-Euler plane beam element, the beam finite elements are obtained by subdividing the beam members longitudinally. Beams 1 and 2 are then modeled using 40 and 20 elements. Each node has two DOFs of transverse displacement and slope. The two beams have an identical elastic modulus E=200 GPa and a unit mass of 7,860 kg/m3. Undamaged beams 1 and 2 are 2 m and 0.8 m in length, and the cross sections are 100 × 12 mm and 100 × 10 mm, respectively. Two damages of beam 1 are located at elements 8 and 30 with the same cross section of 100 × 11.5 mm, and damage to beam 2 is located at element 5 with a cross section of 100 × 9 mm. This work assumes a Rayleigh damping of the stiffness matrix and a proportional constant of 0.0002 and assumes the beam structure is subjected to the base excitation in Figure 2(b). Based on these variables, this numerical application considers the validity of the proposed method originating from the displacement responses contaminated by 3% external noise in (17).Figure 8
A crossed-beam structure: (a) a beam structure, (b) finite elements of beam 1, and (c) finite elements of beam 2.
(a)
(b)
(c)Utilizing the 50% overlap samples of eight segments with 500 sampling record lengths with a rectangular window, the PSEs were extracted from the displacement responses in the time-domain. Figures9(a) and 9(b) represent the PSE curves of beams 1 and 2 transformed from the displacement responses in the time-domain without the external noise, respectively. The first resonance frequency of beams 1 and 2 exists at the frequency of 16 Hz/sample. The PSE curvatures are estimated by the second-order difference method. Taking the maximum PSE values at the first resonance frequency and calculating the curvatures, the resulting curvature curves are plotted as shown in Figures 9(c) and 9(d). It is observed that the damaged elements correspond with the locations to represent the abrupt curvature change, except the joint nodes used to connect the two beams.Figure 9
PSE and the curvature curves of noise-free structure: (a) PSE curves of beam 1, (b) PSE curves of beam 2, (c) PSE curvature of beam 1 using the maximum PSEs at the first resonance frequency, and (d) PSE curvature of beam 2 using the maximum PSEs at the first resonance frequency.
(a)
(b)
(c)
(d)Utilizing the numerically simulated data contaminated by 3% external noise and taking the same process as the previous case, the PSE curves in Figures10(a) and 10(b) display more complicated relations than Figures 9(a) and 9(b). The first resonance frequency of 16 Hz/sample coincides with the PSE curve without the external noise. The curvature curves represent the explicit curvature variations at elements 8 and 30 of beam 1 in Figure 10(c) and element 5 of beam 2 in Figure 10(d). It is observed that the proposed method can be utilized in detecting beam structure damage despite the existence of external noise.Figure 10
PSE and the curvature curves containing 3% noise: (a) PSE curves of beam 1, (b) PSE curves of beam 2, (c) PSE curvature of beam 1 using the maximum PSEs at the first resonance frequency, and (d) PSE curvature of beam 2 using the maximum PSEs at the first resonance frequency.
(a)
(b)
(c)
(d)


## 4. Beam Test

A beam test was performed to detect the damage of a simply supported beam, as shown in Figure11. The gross cross section of the beam is b×t=100 mm × 11 mm, and the net length of the span between end supports is 1200 mm. The damage is located at 350 mm from the left support, which is between nodes 2 and 3. The damaged cross section was established as b×t=100 mm × 8 mm. The eight points measured by accelerometers are numbered, as shown in Figure 11. The beam is excited by hitting an impact hammer with supersoft tip on the support to collect response data within the low frequency range. The experiment was conducted using DYTRAN model 3055B1 uniaxial accelerometers and a miniature transducer hammer (Brüel and Kjaer model 8204) for the excitation of the system. The data acquisition system was a DEWETRON model DEWE-43. The acceleration response data in the time-domain were measured by the DEWETRON. The input acceleration excited by the impact hammer was not measured.Figure 11
Test beam: (a) plan view, (b) section, (c) photo of test beam, and (d) impact hammer.
(a)
(b)
(c)
(d)Figure12(a) represents the acceleration responses in the time-domain at all nodes for one second after hitting. Figures 12(b) and 12(c) display the FFT and PSE curves at all nodes in the frequency-domain to transform the acceleration responses. The maximum energy is revealed at 17.1 Hz. Figures 12(d) and 12(e) indicate the maximum values of FFT at all nodes corresponding to the first resonance frequency and their curvature, respectively. Observing the abrupt change in the curve taken from the maximum absolute values of the FFT, the damage can be detected by the response curve itself only, unlike the results of the numerical application. However, the curvature curve more clearly indicates the damage region near node 3, despite the existence of external noise. Figures 12(f) and 12(g) represent the PSE curve and its curvature, respectively, using the data samples with 50% sample overlap, 0.55 Hz frequency resolution, and the rectangular window function. The plots are similarly interpreted as the results of using the FFT. The numerical and beam tests indicate that the PSE curvature method can be effectively utilized for detecting damage without the baseline or other input data, despite the existence of external noise.Figure 12
Experimental results: (a) acceleration responses in time-domain, (b) absolute values of FFT in frequency-domain, (c) PSE curves, (d) maximum values of FFT nearby the first resonance frequency, (e) their curvature, (f) curve of maximum PSE, and (g) its curvature.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 5. Conclusions

Practical damage detection was performed using output-only response data without baseline data from the intact state and without other input data. This study proposed a damage detection method to trace damage by the curvature of the PSE using the Welch method to reduce the noise effect. The validity of the proposed method was evaluated according to the rectangular, Hamming, Bartlett, Hann, and Blackman window functions with the Welch method, overlapping of 40% and 50%, and frequency resolution of 0.5 Hz and 1.0 Hz through a numerical application and a beam test. The proposed method can be effectively utilized for detecting damage without the baseline or other input data, despite the existence of external noise.



---

*Source: 8761249-2016-01-06.xml*