Multivariate empirical mode decomposition–based structural damage localization using limited sensors

In this article, multivariate empirical mode decomposition is proposed for damage localization in structures using limited measurements. Multivariate empirical mode decomposition is first used to decompose the acceleration responses into their mono-component modal responses. The major contributing modal responses are then used to evaluate the modal energy for the respective modes. A damage localization feature is proposed by calculating the percentage difference in the modal energies of damaged and undamaged structures, followed by the determination of the threshold value of the feature. The feature of the specific sensor location exceeding the threshold value is finally used to identify the location of structural damage. The proposed method is validated using a suite of numerical and full-scale studies. The validation is further explored using various limited measurement cases for evaluating the feasibility of using a fewer number of sensors to enable cost-effective structural health monitoring. The results show the capability of the proposed method in identifying as minimal as 2% change in global modal parameters of structures, outperforming the existing time–frequency methods to delineate such minor global damage.


Introduction
Structural health monitoring (SHM) has emerged as a foundational concept in the management and rehabilitation of critical civil infrastructure. The aging infrastructure, varying environmental effects, and other factors, such as increasing traffic conditions, have created a burden on maintaining the desired residual capacity of the structures. The encumbrance on civil structures because of such factors inherently induces damage in the structures, thereby compromising structural integrity and demanding timely maintenance. Damage identification and its localization (Burgos et al., 2020) are one of the critical steps toward efficient monitoring and maintenance of the structures. Through continuous monitoring, damages can be detected in the structures to prevent their major catastrophes at the later stage. In this article, a new and cost-effective damage localization technique is proposed to detect and isolate damage using a limited number of sensors.
Vibration-based SHM techniques (Barthorpe and Worden 2020;Cawley 2018;Erazo et al., 2019;Gatti 2019;Okayasu and Yamasaki 2019) offer viable options for tracking damages in the structures based on the measured data. Existing damage identification techniques involve analysis in time domain, frequency domain, and timefrequency (TF) domain (Sirca and Adeli 2012) along with various artificial intelligence techniques (Salehi and Burgueno 2018;Ying et al., 2012). In a structure, the damage is initiated because of changes in the structural or boundary condition of the system. As these properties change with time, the structure is mathematically modeled as a time-varying system (Chen et al., 2019;Zhou and Li 2017). Damage could be either discrete (because of a sudden shock) or progressive (because of a frequent accumulation of instantaneous damages) that can lead to a catastrophic failure. The time-domain methods only provide information about the time characteristics of the systems. Hence, they are prone to noise contamination and environmental factors, which offer challenges for large structures. On the other hand, time-series methods (Bernal et al., 2002;Johnson et al., 2004;Lei et al., 2003;Nair et al., 2006) require model order selection which is one of the limitations of these methods, especially for large-scale measurements with noise. Frequency-domain methods provide information about the frequency characteristics of the systems. Although frequency changes in the structure can be associated with the presence of damage, discrete changes in natural frequencies may not be sufficient for unique identification of the location of structural damage as a crack at two different locations may cause the same frequency change irrespective of its location (Salawu 1997). Time-frequency method (Neild et al., 2003) can identify the frequency components of the signal and apprehend their time-variant features.
Popular TF methods are short-time Fourier transform (STFT) (Daubechies 1990), Wigner-Ville distribution (WVD) (Andria et al., 1994), and empirical mode decomposition (EMD) (Barbosh et al., 2020;Huang et al., 1998). Short-time Fourier transform is used to increase the sparsity of the signals in TF domain, and the TF resolution of STFT is inversely proportional to window length. Increasing the window length improves the frequency resolution, whereas it decreases its frequency tracking capability. On the other hand, WVD has an unwanted crossproduct disturbance along with the frequency bands, and there is a challenge to reconstruct the signal. Therefore, conventional TF analysis suffers from a trade-off between time and frequency resolution. Moreover, vibration data attribute highly nonlinear and nonstationary behavior due to various natural hazards such as wind and earthquakeinduced excitation, where traditional TF analysis is not adequate. To address these issues, the nonlinear and nonstationary nature of vibration data was analyzed by Hilbert-Huang transform (HHT) (Huang et al., 1998), which is basis-free in nature. The HHT method can decompose any complex signal into a finite number of intrinsic mode functions (IMFs), which can be used to estimate modal parameters. A variant of EMD, multivariate EMD (MEMD) was recently explored as a modal identification tool using multisensor data (Barbosh et al., 2018;Sadhu 2015). Recently, Sony and Sadhu (2020) used a hybrid method combining MEMD and synchrosqueezing transform to study the behavior of time-varying structural systems. However, the important aspects of damage localization were not considered. In this article, a novel MEMD-based damage localization method is proposed that uses a limited number of sensors.
There have been continuous improvements in sensors and sensing technology as they constitute a significant portion of cost during SHM applications and contribute to the overall accuracy of condition assessment. Traditionally, wired sensors are used as a dense array distributed over the structure to acquire long-term SHM data. However, this approach is not a cost-effective and viable option for largespan bridges or tall buildings because of labor-intensive cable installation. The setbacks of wired sensors were envisioned to solve using smart wireless sensors (Abdulkarem et al., 2019;Lynch, 2007). Capabilities and local processing of wireless sensors in the decentralized framework were much later exploited by Spencer et al. (2004). Recently, other modern sensors such as cameras, robotic sensors, smartphones, and drones have been used for SHM through the processing of images and videos (Feng and Feng 2018;Sony et al., 2019). However, such SHM approach is cost-effective only if the hidden structural information and damage characteristics are accurately assessed from the limited measured data.
An optimal number of limited sensors (Li et al., 2014) that can provide the same information as the large array of sensors will reduce the cost of SHM. There has been considerable research on system identification by limited sensors using two popular methods: non-sparse and sparse techniques. Non-sparse methods have gained interest in the operational modal analysis as a nonparametric alternative to the structural identification from output-only measurements (Hazra et al., 2011;Sadhu et al., 2017Sadhu et al., , 2019. A parallel factor decomposition along with Bayesian model updating was used for underdetermined modal identification (Abazarsa et al., 2013). In addition, Guo and Kareem (2016) proposed spatial TF distribution to handle nonstationary response. However, the proposed method highly depends on the quality of the selection of a single auto-spectral component, and its poor selection could lead to inaccurate results.
Unlike non-sparse methods, sparse techniques use TF methods for sparse representation of signals. Yang and Nagarajaiah (2013) used sparse component analysis (SCA) along with L1-minimization to improve underdetermined modal identification. The study by Yu et al. (2014) explored SCA-based methods in the TF domain to estimate time-varying modal parameters and validated the proposed method under thermal effects. Amini and Hedayati (2016) used STFT and SL0 algorithm to perform underdetermined system identification using earthquake and ambient excitation and showed the robustness under noise to separate closely spaced modes. Moreover, in Yao et al. (2018), the authors explored SCA for modal identification using limited sensors and validated the method on both stationary and nonstationary response of structures. Castiglione et al. (2018) proposed frequency banding for largely underdetermined scenarios by decomposing a large underdetermined problem into several overdetermined problems. The method operates directly in the frequency domain and analyzes the cross-spectral matrix of the data. However, user intervention during frequency banding and manual selection of estimated modes is one of the limitations of the proposed method. Yi et al. (2018) proposed a novel SCA method for estimating the number of active modes using statistical properties of normalized single source point vector.
Although the above sparse and non-sparse system identification methods show a wide variety of applications in modal identification using limited sensors, damage localization has not given its due attention under limited sensor measurements (Gomes et al., 2018). Amini and Karami (2012) recently used the system Markov parameters for damage detection. The authors evaluated the effect of noise, number, and location of the sensors. The study was based on finite element study and predefined damage locations and did not identify the unknown damage location. Bagheri et al. (2016) used an optimization algorithm to localize and quantify damage in seismically excited structures using a limited number of sensors. A competitive optimization algorithm combined with moment generating function was used as a damage indicator. The limitation of the study lies in using the moment of the segment as a damage indicator while identifying the instance of the damage. Recently, Karami-Mohammadi et al. (2020) used mode shapes coupled with wavelet transform for damage identification while using a limited number of sensors. The optimal spatial location of the sensors was achieved using the minimization of the nondiagonal entries in the modal assurance criterion matrix. The study is dependent explicitly on mode shapes, and ambiguity arises on interpreting the damage location with respect to different damage scenarios. Moreover, there is a significant challenge in accurate identification of mode shapes using real data with measurement noise.
In this article, MEMD-based damage localization is proposed by taking advantage of the modal energies of individual modes extracted from the measured data. Because of the capability of handling a limited number of sensors, MEMD is further explored to compare the performance of the damage localization using a suite of limited sensors. The article is organized as follows. Section 1 provides a brief introduction to this topic, illustrating the existing literature and the gap areas. Section 2 presents the proposed methodology and the corresponding theoretical background. Section 3 shows the numerical studies using a 10-DOF model subjected to a wide extent of damage and measurements obtained from a suite of limited sensors. A full-scale study comprising of a real bridge with two different damage scenarios is presented in Section 4, followed by key conclusions in Section 5.

Background of MEMD
A brief background of MEMD method is presented first before discussing the details of the proposed algorithm. Multivariate EMD is a variant of EMD (Huang et al., 1998) and EMD works only for a single measurement. The basic idea of EMD is that any complicated dataset can be decomposed into a finite number of "intrinsic mode functions (IMFs)." This decomposition method is adaptive and, therefore, free of any basis function. Because the decomposition is based on the local characteristics of the data, it is applicable to nonlinear and nonstationary processes. In the presence of multiple sensors, EMD faces several challenges (Sadhu, 2015) such as (a) there is no guarantee that the extraction of IMFs from different channels of measurements will match, either in number or in their frequency contents, due to their sole dependencies on a single measurement, and (b) the joint information between multichannel measurements is not considered because signals from multiple sensors are treated individually. For multichannel measurements, a multivariate extension (i.e., MEMD) is proposed by Rehman and Mandic (2010). Because of the geometric projection and averaging tasks, MEMD is free of any linear algebraic operations and allows solving of underdetermined system identification using limited sensors.
Vibration-based monitoring of real-life structures requires use of multiple sensors and structural modes that are closely spaced. An array of sensors for a large-scale structure is expensive and takes a considerable portion of budget of the condition assessment. In this article, MEMD is used to reduce the number of sensors by decomposing the signal into its IMFs. A suitable set of direction vectors are sampled on unit hyperspheres (n-spheres) based on both uniform angular sampling methods and quasi-Monte Carlobased low-discrepancy sequences. To estimate all multidimensional envelopes, the multiple signals are projected onto the chosen direction vectors and the average of all envelopes is considered as a local mean of multiple signals. The algorithm for MEMD (Rehman and Mandic 2010) is presented in Algorithm 1. Algorithm 1: MEMD Input: A signal y(t) Result: Intrinsic mode functions (IMFs) Initialization for y(t) = signal do 1. Evaluate the direction vector, D and its projections along k À th direction; 2. Determine the k-th projection, p k (t) of the input signal y(t) over the k-th direction vector, Y k , for all k (k=1, 2,….., L, where L is the number of direction vectors D). 3. Estimate the corresponding time t k i of maximum p k (t) for all k. 4. Interpolate ½t k i , yðt k i Þ to extract multi-dimensional envelopes, ϵ k (t). 5. Determine the mean of envelope W ðtÞ as, 6. Estimate the residual R(t) using RðtÞ ¼ yðtÞ À W ðtÞ.
If R(t) satisfies the stopping criteria of multivariate IMF, repeat the above steps to [y(t) À R(t)] until the next IMF is diminished. Otherwise, repeat it to R(t). End

Proposed algorithm
Consider a linear dynamical system with n degree-offreedom (DOF) subjected to a broadband excitation X(t), the equation of motion is expressed as

My
:: ðtÞ þ C _ yðtÞ þ KyðtÞ ¼ ðtÞ where y(t) is the displacement vector. Using the classical modal superposition theorem, the solution to equation (2) for those of broadband X(t) can be written in terms of an expansion of vibration modes where y and ν are the response and modal response matrix, respectively. Φ m×n is the modal transformation matrix. n and m are the number of modal responses and measurements, respectively. The measurement at the i-th DOF (i = 1, 2, …, m) can be expressed as where ν j is the j-th modal response and f ij represents the mode shape ordinate of i-th DOF and j-th mode. Because MEMD results in IMFs containing mono-component modal responses, each signal of y(t) can be expressed in terms of its IMFs where I ij is the j-th IMF of y i (t). Comparing the equality of equations (4) and (5), we get Each IMF I ij is then analyzed to obtain Fourier spectrum (E ij (f)), as shown in equation (7) E ij ðf Þ ¼ If E ij (f) is the Fourier transform of I ij (t), then |E ij (τ)| 2 can be interpreted as energy density at the frequency τ, which means that the total energy contained in a small frequency interval [τ À ϵ, τ + ϵ] around τ is approximately given by 2ϵ|E ij (τ)| 2 . However, as per mathematical axiom, if |x| > | y|, then, for every x≠0, y≠0, and x, y > 0, x 2 > y 2 is true, and thus, the peak Fourier amplitude can be interpreted as modal energy. The damage localization feature is then extracted using percentage change in peak amplitude of Fourier spectrum of modal responses as shown in equation (8) where ΔE ij represents absolute percentage change in modal energy of the corresponding damaged (f D ) and undamaged (f U ) modal frequencies. Then, ΔE ai is calculated as the average of ΔE ij for all selected IMFs of i-th sensor Finally, the threshold ΔE is proposed by taking mean value of ΔE ai for all the selected sensor locations The proposed algorithm is applied by first decomposing the multi-sensor data into their IMFs, and then, their Fourier spectra are used to identify ΔE ij of the corresponding frequencies and are used as damage localization feature. The damage is deemed present for any sensor location, if ΔE ai is higher than ΔE. The flowchart of the proposed algorithm is presented in Algorithm 2.

Numerical studies
A 10-DOF model is used to illustrate the performance of the proposed method. The model is subjected to a broadband earthquake (i.e., Imperial Valley earthquake) with a peak ground acceleration (PGA) of 0.1 g to study the damage localization using the earthquake-induced measurements. Considering the fact that the first floor columns are primarily damaged during the earthquakes, three damage scenarios are used to evaluate the performance of the proposed method under varying level of damage in the first floor. For example, C1 indicates undamaged model with first five natural frequencies as 0.78, Table 2. Modal energies (E ij (f)) of IMF-1 and IMF-2 of the 10-DOF model subjected to Imperial Valley earthquake.       frequencies are less than 5% in C2 and C3 (i.e., 1.9% and 4.5%, respectively), where the traditional TF methods  were unable to delineate such minor change in the modal frequencies. The models are subjected to above earthquake and the earthquake-induced responses are analyzed to identify the location of the simulated damage.
The 10-DOF model is excited using Imperial Valley earthquake and damage is induced at the first floor (i = 1) by  reducing the stiffness, and the proposed algorithm is evaluated to localize damage using a combination of all sensors. MEMD is implemented to decompose the simulated responses into their IMFs. The first two modes are illustrated to show the change in the natural frequencies and ΔE ai is used as a damage indicator. The values of ΔE ij for various floors of I i1 (IMF-1) and I i2 (IMF-2) are tabulated in Tables 2 and 3, respectively. For example, Table 2 provides the peak Fourier amplitude of first floor response as 8.2 unit for C1, as shown in Figure 1(a). In Figure 1-3, the first row shows IMF-1 and second row shows IMF-2, respectively. The corresponding modal energies for C1, C2, and C3 are shown in Figures 1-3, respectively.
To localize the damage, equation (8) is used to calculate ΔE ij , followed by evaluating ΔE ai using equation (8) for all selected IMFs of a sensor. Finally, the average of ΔE ai for all sensors is evaluated using equation (9) and is used as a damage localization feature, as presented in Table 3. The percentage change between C1 and C2 is 6.1% for IMF-1, similarly, for C1 and C3, its 51.2%. The damage localization feature is then evaluated by taking mean of ΔE ij of all IMFs. For example, C2 yields ΔE ai as 32.5 unit, which is equal to average of corresponding values of IMF-1 and IMF-2. ΔE ai values for various floors are shown in Figure 4(a). For example, the first floor has a value of 32.5 unit, the fifth floor has a value of 22.7 unit, and 10th floor has a value of 20.7 units. In this sub-figure, ΔE is calculated as 25.3 unit, which is shown in the dotted line and it can be clearly seen that the proposed damage index exceeds ΔE at i = 1, indicating accurate damage localization at the first floor.
The study is further enunciated by considering a limited number of sensors, whereas keeping sensor locations 1, 5, and 10 intact for all the cases to have consistent reference sensors for the comparison purposes. For example, Figure 4 shows the identification results of C2 for a limited sensor cases using 10, 9, 7, and five sensors, respectively. The damage location can be clearly identified in floor 1; however, in case of eight and six sensors, the proposed method could not classify accurately. In case of C3, all the sensors classify the damaged floor accurately, as shown in Figure 5. Therefore, the sensitivity of the proposed method improves with the severity of the damage.   Table 4, where 0 and 1 indicate index for misclassification and accurate classification, respectively. It can be concluded that the proposed method is sensitive to the severity of damage, the higher the severity, the higher is the accuracy for damage localization. However, the proposed method has nearly 92% accuracy across all limited sensor cases in C2 and C3 subjected to the example earthquake, indicating its efficacy to accurately identify less than 5% global change in the modal frequencies.

Full-scale study
A benchmark study on Z24 bridge is used to evaluate the performance of the proposed method. The bridge is located in Canton Bern, Switzerland, connecting Koppigen and Utzenstorf. The bridge is a highway overpass linking Bern and Zurich and it is a prestressed bridge, with three spans, two lanes, and 60 m overall length (Maeck and Roeck De (2003); Teughels and Roeck (2004)). The bridge was demolished at the end of 1998 because a new railway adjacent to the highway required a new bridge with a larger side span. The bridge was excited by two shakers, one at the mid-span of the bridge and another at the side span. Because of the size of the bridge, response was measured in nine setups with each setup containing up to 15 sensors each. In all setups, three accelerometers and the two force sensors were common. The data were sampled at 100 Hz, and a total of 65,536 samples were acquired.
The study has multiple damage cases; however, for the purpose of the damage localization, only three cases are considered with a total number of eight sensors, as shown in Figure 6. The pier settlement case was used to localize the damage and differentiate the damaged pier from the undamaged pier. The settlement was simulated by cutting the Koppigen pier and removing about 0.4 m of concrete. Lowering and lifting was applied by six hydraulic jacks. During the tests, the pier rested on steel sections with similar stiffness as the uncut concrete section. The lowering of the  pier was carried out by supporting the structure with scaffolding. The pier was cut to support its dead weight, test equipment, and the impact of a vehicle with and without normal force. The base plates were located and connected using shear connectors. One of the piers (Koppigen) among two piers, as shown in Figure 6, was damaged by lowering the piers. In this study, first, an undamaged case was used as a baseline and then a pier settlement of 20 mm, 40 mm, and 95 mm are used as the damaged cases. Schematic and sensors of the undamaged and damaged piers used in this study are shown in Figure 6. These data were made publicly available by researchers at the Katholieke Universiteit Leuven and is available at: https://bwk.kuleuven.be/ bwm/z24.
The damaged pier is localized using the proposed algorithm. Various cases of limited sensors are used to illustrate the performance of the proposed algorithm. There are four sensors each on both of the piers (a total of eight sensors), and for limited sensor case, the number of sensors was reduced to six and four, respectively. ΔE of each sensor is used to present the damage localization feature for various cases. Three different damage localization scenarios of pier settlement are considered, namely 20 mm, 40 mm, and 95 mm of settlement to include the severity of damage in the analysis. As shown in Figure 7, for 20 mm pier settlement, it can be observed that ΔE ai of the damaged pier is consistently higher than the ΔE in all limited sensor cases. Similar analysis is conducted for 40 mm and 95 mm lowering of pier as shown in Figures 8 and 9, respectively. The results show that the proposed method is able to detect the damage with various severities as well as identify the location of damage using a suite of limited sensor cases across various levels of damage.
The variability for 95 mm pier settlement in ΔE ai for each mode of every sensor is presented as box plots and for various limited sensors damage cases, as shown in Figures  10 and 11(a) and (b), respectively. It can be observed that along with mean of ΔE ai for each sensor, the variability in the modal ΔE ai for the damaged cases is larger than ΔE, indicating clear delineation of the damage location at the damaged pier. It can be observed that high variation is in the sensors from the damaged pier (i.e., 511, 521, 531, and 541) and low variation in the undamaged pier (i.e., 411, 421, 431, and 441). Moreover, the damage index of the damaged pier is always higher than the threshold value irrespective of the number of limited sensors, indicating the efficacy of the proposed method to identify the damage location.

Conclusions
In this article, a novel damage localization method is proposed using a damage index obtained from MEMD. Multivariate EMD is used to decompose multi-sensor data into their mono-components. The mono-component modal responses are further used to evaluate the modal energy to derive the damage localization feature using limited sensors. The threshold for the proposed damage index is adaptive and automated in nature and free of any user discretion or any predefined value, which is one of the merits of the proposed method. To the best of the authors' knowledge, the capabilities of MEMD in performing damage localization using limited sensors have not been undertaken in the literature. The proposed method is demonstrated using a suite of numerical studies and the benchmark data of the Z24 bridge. It is concluded that the proposed method works well on all the studies and is effective in localizing the damage. The limited measurement aspect of damage localization is explored by selecting a fewer number of sensors and it is shown that with limited measurements, the proposed method is as effective as total number of measurements equals to the DOF of the model. The results show the capability of the proposed method in identifying as minimal as 2% change in global modal parameters of structures, outperforming the existing TF methods to delineate the minor global damage.