Laplace approximation in sparse Bayesian learning for structural damage detection

Mechanical Systems and Signal Processing 140 (2020) 106701

Contents lists available at ScienceDirect

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

Laplace approximation in sparse Bayesian learning for structural damage detection Xiaoyou Wang a, Rongrong Hou a, Yong Xia a, Xiaoqing Zhou b,⇑ a b

Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Kowloon, Hong Kong, PR China College of Civil Engineering, Shenzhen University, Shenzhen, PR China

a r t i c l e

i n f o

Article history: Received 8 June 2019 Received in revised form 12 November 2019 Accepted 31 January 2020

Keywords: Sparse Bayesian learning Structural damage detection Laplace approximation Vibration based methods

a b s t r a c t The Bayesian theorem has been demonstrated as a rigorous method for uncertainty assessment and system identification. Given that damage usually occurs at limited positions in the preliminary stage of structural failure, the sparse Bayesian learning has been developed for solving the structural damage detection problem. However, in most cases an analytical posterior probability density function (PDF) of the damage index is not available due to the nonlinear relationship between the measured modal data and structural parameters. This study tackles the nonlinear problem using the Laplace approximation technique. By assuming that the item in the integration follows a Gaussian distribution, the asymptotic solution of the evidence is obtained. Consequently the most probable values of the damage index and hyper-parameters are expressed in a coupled closed form, and then solved sequentially through iterations. The effectiveness of the proposed algorithm is validated using a laboratory tested frame. As compared with other techniques, the present technique results in the analytical solutions of the damage index and hyper-parameters without using hierarchical models or numerical sampling. Consequently, the computation is more efficient. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction Vibration-based damage detection techniques have received considerable attention in the past decades [1–4]. However, several challenges still exist such as the ill-posed inverse problem, limited measurements and unavoidable uncertainties [5–10]. Specifically, a small perturbation in the measurement data may cause a significant variation in the solution [11]. The incomplete measurements may lead to underdetermined or nonunique solutions [12]. The uncertainties coming from the estimation errors, measurement noise, varying ambient conditions and so on may result in inaccurate identification results [7,9]. The Bayesian theorem has been demonstrated as a rational technique to handle the aforementioned problems [13–16]. The method possesses a rigorous framework in which all unknowns and uncertainties are regarded as stochastic quantities, and represented by variables or parameters following designated probability distributions [14,15]. Rather than treating the problem deterministically, the Bayesian formulation provides a posterior probability distribution over all plausible values by multiplying the likelihood function with the prior. In recent years, the sparse Bayesian learning has been developed for structural damage detection [12,17–19], on the basis that damage usually occurs at limited positions or members in the ⇑ Corresponding author. E-mail addresses: [email protected] (X. Wang), [email protected] (R. Hou), [email protected] (Y. Xia), [email protected] (X. Zhou). https://doi.org/10.1016/j.ymssp.2020.106701 0888-3270/Ó 2020 Elsevier Ltd. All rights reserved.

2

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

preliminary stage of structural failure. The technique treats this sparse feature as the prior. All unknowns are solved automatically by employing a full Bayesian analysis which marginalises over all possible values [17], or some special approaches that identify the most probable values (MPVs) such as the maximum a posterior (MAP) [12] and the maximum likelihood (ML) estimation [20]. For example, Huang et al. [17] proposed a hierarchical sparse Bayesian learning for structural health monitoring using incomplete modal data, in which a full Bayesian analysis is conducted by considering the posterior uncertainty based on the Gibbs sampling. Hou et al. [12] proposed a sparse Bayesian learning method for probabilistic structural damage detection and adopted an iterative procedure to derive the MPVs of unknowns. However, the application of the Bayesian framework to the damage detection of real structures suffers from common limitations. For most problems an analytical posterior probability density function (PDF) is not available due to the computationally prohibited integration in the evidence, which is resulted by the high-dimension and the nonlinear relationship between the measured modal data and structural parameters [21]. In the case wherein the prior PDF does not conjugate the likelihood function [22,23], it is intractable to directly obtain an analytical posterior distribution either. Several techniques have been developed to circumvent these difficulties. One technique is the hierarchical Bayesian model (HBM) [24]. HBM possesses a multi-level structure in which the parameter in the prior PDF relies on other introduced parameters. These parameters are called hyper-parameters and are assumed to follow designated distributions. There has been a lot of applications of HBM in diverse fields, e.g., compressive sensing, signal detection, model updating and damage detection [24–27]. With this method, the nonlinear problem can be extended as the product of a series of coupled hierarchical linear PDFs, which simplify the calculation and make the posterior PDF analytically tractable [18]. However, the HBM is of low efficiency given the need to estimate a series of additional parameters introduced by hyper-priors. Moreover, the successful application of HBM requires proper assumptions of all hyper-priors, which is not guaranteed owing to the demand for sufficient information about individual hyper-parameter. Given that the uncertainties of the posterior are directly related to the variance of the prior [24], the inappropriate assumptions of hyper-priors tend to increase the uncertainties of the prior, further reduce the robustness of the model and lead to a decrease in the sensitivity of the posterior PDF to the prior information. This phenomenon worsens with the increase of hierarchical layers. Another technique is the Markov Chain Monte Carlo (MCMC) simulation. With this method, samples consistent with the posterior probability distribution can be generated [28,29]. Then the asymptotic solution of the posterior PDF is obtained from these samples. One advantage of the MCMC approach is that it can provide a full characterization of the posterior PDF instead of merely focusing on the MPVs [28]. Hence, the MCMC is convenient to calculate the complicated integral and evaluate the posterior uncertainties. However, the basic MCMC technique is extraordinarily time-consuming for models containing a large number of parameters or involving partial differential equations [29]. To improve the computational efficiency, the Metropolis–Hastings algorithm [28] and the Gibbs sampling technique [17], as special cases of MCMC, have been developed to deal with multi-dimensional problems. In the Metropolis–Hastings algorithm, the proposal distribution and the stopping criterion for sampling largely affect the accuracy of the uncertainty evaluations [28], whereas they are difficult to ascertain. Huang et al. applied the Gibbs sampling technique to a hierarchical model for structural damage assessment [17]. While the Gibbs sampling is only applicable to standard distributions and converges slowly for many hierarchical models [24]. Moreover, the ultimate solutions may be unstable and nonunique as a result of stochastic sampling. The Laplace approximation is another technique that can be employed to solve the computationally prohibited integral in the evidence [15]. In contrast to the purely numerical methodologies such as MCMC, the Laplace approximation requires less computational cost and possesses a unique solution [30]. This technique leads to the asymptotic solution of the problem directly without the introduction of additional calculations, therefore can effectively avoid the risk of increased uncertainties caused by HBM. The Laplace approximation technique has been utilized by many researchers to develop asymptotic expansions of posterior interests and expectations, especially in the Bayesian framework [14–16,19,31–33]. In this regard, this study extends the sparse Bayesian model for structural damage detection using the Laplace approximation. The complicated equation to be integrated in the evidence is approximated as a Gaussian PDF. The solutions of the damage index and hyper-parameters are then expressed in closed form and solved iteratively. An experimental frame is utilized to verify the effectiveness of the proposed algorithm. The identification error and computational cost of the present algorithm and the previous proposed expectation–maximization (EM) algorithm [12] are compared. 2. Bayesian probabilistic framework Suppose that Nm vibration modes of a structure are available such that we have a vector of eigenvalues h i h i NP N m b¼ / b ;/ b ;;/ b b 2 RNp (r = 1,. . .,Nm) is b k2 ;    ; b k Nm 2 R1Nm and a matrix of mode shapes w , where / k¼ b k1; b 1 2 Nm 2 R r the rth mode shape vector at N p measured degrees of freedom (DOFs). 2.1. Definition of a model class The model class M is based on a set of linear structural models, each having the same known mass matrix M and an uncertain stiffness matrix K, which is parameterized as a combination of n elemental stiffness matrix Ki (i = 1,2,. . .,n) and calculated using the equation [11]

3

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

K¼

n X

ð1Þ

si Ki

i¼1

where n is the number of elements, Ki is the ith elemental stiffness matrix, and si is the stiffness parameter. Suppose that the structural mass remains constant in the damaged state, while the stiffness matrix is reduced to K

K¼

n X

ð2Þ

si Ki

i¼1

where si denotes the stiffness parameter in the damaged condition. The stiffness reduction factor (SRF) is then defined as [11]

hi ¼

si  si si

ð3Þ

where the value of hi is greater than 1 and no more than 0, particularly, a value of 0 denotes no damage of the element and 1 represents a complete damage. Thus hi can be regarded as the damage index that indicates both damage position and severity. The analytical eigenvalues and mode shapes are governed by the eigenequation defined as



 K  kr M /r ¼ 0

ð4Þ

where kr and /r refer to the rth analytical eigenvalue and mode shape, respectively. 2.2. Bayesian model updating framework With the selected model class and available measurements, the Bayesian formulation can be expressed as [14]

pðhjD; MÞ ¼ c1 pðDjh; MÞpðhjMÞ where D denotes the measured modal parameters including

ð5Þ D

E b b ; h denotes the damage index defined as above; kr ; / r

pðhjD; MÞ is the posterior PDF of damage index given the measured data D and identified model M; pðDjh; MÞ is the likelihood function representing the PDF of modal data given the damage index; pðhjMÞ is the prior PDF; and c ¼ pðDjMÞ is the normalization constant termed evidence. 2.3. Likelihood function for structural modal parameters In order to simplify the notation, the dependence on M is herein dropped in the notation of likelihood function, expressed as pðDjhÞ. Given the damage index, the modal parameters are deemed independent mode by mode. Therefore, Nm      Y    b b jh pðDjhÞ ¼ p b kjh p wjh ¼ p b k r jh p / r

ð6Þ

r¼1

The measurement errors, er and er , are assumed to follow the Gaussian distribution with a zero mean and diagonal variance matrix, which is respectively given by

er ¼

b   k r  kr  N 0; b1 b kr

  b  /  N 0; c1 I er ¼ / r r

ð7Þ ð8Þ

where hyper-parameters b and c equal the reciprocal of the variance, representing the precision of the measured eigenvalues and mode shapes, respectively. b are given by The resulting likelihood functions of h based on the measurement b k and /

8 " #2 9 Nm < bX =    b N2m b k  k ð h Þ r r p b kjh; b ¼ exp  b : 2 r¼1 ; 2p kr

( ) Nm    Np N2 m cX b  / ðhÞk2 b c ¼ c exp  k/ p /jh; r r 2 2p 2 r¼1

ð9Þ

ð10Þ

4

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

2.4. Prior PDF of damage index In the Bayesian framework, the damage index h is assigned with a prior distribution pðhjMÞ according to engineering judgements and knowledge. Considering that the l1 regularization technique is not differentiable at hi ¼ 0 ði ¼ 1; 2; . . . ; nÞ which may increase the difficulty of the calculation, hence, the automatic relevance determination (ARD) prior [13] following the Gaussian distribution is utilized in this study and serves as the regularization item for the ill-posed problem. In the ARD model, each damage index is assigned with an individual hyper-parameter. Therefore, the prior PDF is expressed as the product of n independent Gaussian distributions, each corresponding to one damage index, namely,

pðhjaÞ ¼

n Y i¼1

 pðhi jai Þ ¼

1 2p

n2 Y n

i¼1

1



1 2



a2i exp  ai h2i

ð11Þ

where hyper-parameter ai represents the precision of hi . 2.5. Posterior PDF of damage index The posterior PDF, updated with the likelihood PDF and prior PDF according to Eq. (5), is expressed as

      b a; b; c ¼ c1 p b b c pðhjaÞ p hjb k; w; kjh; b p wjh;

ð12Þ

  b hja; b; c with respect to h, that is, Here the evidence is calculated by integrating p b k; w;

Z    Z      b a; b; c ¼ p b b hja; b; c dh ¼ p b b c pðhjaÞdh c ¼ pðDjMÞ ¼ p b k; wj k; w; kjh; b p wjh; 8 ) ! " #2 Z  N2m  Np Nm  n2 Y Nm Nm n n < bX b X 1  2 k r  kr ðhÞ b c 2 1 cX 2 b  / ðhÞk2  1 ¼ ai exp   k/ a h dh i i r r 2 b : 2 r¼1 2p 2p 2 i¼1 2p 2 r¼1 kr i¼1

ð13Þ

3. Bayesian inference based on Laplace approximation   b a; b; c , that is The damage index h can be solved through the MAP estimate of p hjb k; w;

      b a; b; c ¼ arg max c1 p b b c pðhjaÞ h ¼ arg max p hjb k; w; kjh; b p wjh; !  N2m  Np Nm  n2 Y n 1 b c 2 1 1 2 ai ¼ arg max c 2p 2p 2p i¼1 8 9 " #2 Nm Nm < bX n b  2 = 2 k r  kr ðhÞ cX 1X b  k / r  /r ðhÞk2  a i hi  exp  b : 2 r¼1 ; 2 i¼1 2 r¼1 kr

ð14Þ

For the convenience of optimization, the negative natural logarithm of the right side in Eq. (14) is used, which is written as

" #2 Np h Nm Nm X n i2 X b X X  2 k r  kr ðhÞ b  / ðhÞ þ / h ¼ arg min J ðhÞ ¼ arg min b þc ai hi j;r j;r b kr r¼1 r¼1 j¼1 i¼1

ð15Þ

where the items unrelated with h are omitted for simplification. However, the hyper-parameters fa; b; cg contained in Eq. (15) are unknown and need to be determined first. The MPVs of   b . In particular, k; w the hyper-parameters fa; b; cg can be calculated by maximizing p a; b; cjb

  b a; b; c pða; b; cÞ   p b   k; wj b a; b; c b ¼   /p b k; wj p a; b; cjb k; w b p b k; w

ð16Þ

where we assume that pða; b; cÞ is uniformly distributed as hyper-prior information is unknown. Consequently, the hyper  b a; b; c , or the evidence in Eq. (13). parameters fa; b; cg can be estimated by maximizing p b k; wj Due to the nonlinear relationship between hkr ; /r i and h, it is intractable to obtain a closed-form solution of the evidence. Therefore, an asymptotic or numerical estimation is required, as explained in the Introduction. The Laplace approximation is a desirable technique that can be employed to solve the complicated integral in the evidence, provided that the model is globally or locally identifiable based on available observations [15]. The difference between global and local identifiability

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

5

lies in the number of optimal modes. Beck et al. have proposed that for a locally identifiable system, the PDF should be approximated as a weighted average of the second order Taylor’s series expansion at all optimal modes [15]. By contrast, for a globally identifiable system, the PDF can be expanded merely at the unique optimal mode [14–16]. In Contrast to the traditional deterministic system identification technique, the Bayesian framework in the present study consider uncertainties and quantifies these uncertainties as parameters. In other words, the difference is that the hyper-parameters in the objective function shown in Eq. (15) are stochastic. If the uncertainties are at an acceptable level, the system can be regarded as globally identifiable as long as the measurements are sufficient. Under this condition, the PDF to be evaluated is sharply peaked and can be fully characterized by the unique optimal mode (where the second order Taylor’s series expansion is conducted) and the covariance matrix. Assuming that the system is globally identifiable based on the measurements. By employing the Laplace approximation,   b b hja; b; c , which is to be integrated in Eq. (13), can be approximated as a Gaussian distribution [34], that is p k; w;

         b hja; b; c ffi p b b hja; b; c exp  1 h  h T A h  h p b k; w; k; w; 2

ð17Þ

  b hja; b; c and is equivalent to the MAP estimate calculated by Eq. (15). The matrix A where h represents the mode of p b k; w; denotes the inverse of the variance, which equals the Hessian matrix at h and is defined by

  b h j a; b; c A ¼  rr ln p b k; w;

¼ h¼h

@ J ðhÞ ¼ W þ bH þ cP @hi @hj h¼h

ð18Þ

In the above equation the Hessian matrix is decomposed into three items, i.e., W is a diagonal matrix with entries Wii ¼ ai , and H and P are given by

Hij ¼

"  #2     Nm "  #   Nm Nm b X X b k r  kr h k r  kr h @ 1X 1 @kr h @kr h 1 @kr h  ¼     2 b b b @hi @hj 2 r¼1 @hi @hj kr kr k r @hi @hj r¼1 b r¼1 kr

    Nm N p h   Np h Np Nm X Nm X XX  i2 X  i @/k;r h @/k;r h @/k;r h @ 1X b b /  /k;r h / k;r  /k;r h  Pij ¼ ¼   @hi @hj 2 r¼1 k¼1 k;r @hi @hj @hi @hj r¼1 k¼1 r¼1 k¼1

ð19Þ

ð20Þ

The derivatives of the eigenvalues and eigenvectors with respect to the stiffness parameter can be calculated using the Nelson’s method [35] or the substructural approach [36,37].   b hja; b; c is expressed as A normalized Gaussian distribution proportional to p b k; w;

qðhÞ ¼

   T   1 h  exp  h A h  h ¼ N hjh; A1 n=2 2 ð2pÞ jAj1=2

ð21Þ

where jAj denotes the determinant of the matrix A. On the basis that the integral of a normalized Gaussian distribution   b a; b; c can be solved as equals 1, the evidence p b k; wj

   Z    Z T   1 b a; b; c ¼ p b b hja; b; c dh ffi p b b hja; b; c p b k; wj k; w; k; w; exp  h  h A h  h dh 2   ð2pÞn=2 b hja; b; c ¼p b k; w; 1=2 jAj

ð22Þ

In this way, the asymptotic analytical expression of the evidence is obtained. For the convenience of calculation, the logarithm form is used and shown as

    b a; b; c ¼ ln p b b hja; b; c þ n lnð2pÞ  1 lnjAj ln p b k; wj k; w; 2 2       b h; c þ ln p hja þ n lnð2pÞ  1 lnjAj ¼ ln p b kjh; b þ ln p wj 2 2 ¼

    n n   2 Nm b Np Nm  c  n 1 1X 1X þ þ lnai  ai hi þ ln ln ln 2p 2 2p 2 i¼1 2 i¼1 2 2 2p



"  #2 Nm Nm b k r  kr h bX cX b  / ðhÞk2 þ n lnð2pÞ  1 lnjAj  k/ r r 2 b 2 r¼1 2 2 2 r¼1 kr

ð23Þ

6

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

  b a; b; c . From Eq. (23), setting the As stated previously, the MPVs of fa; b; cg can be estimated by maximizing p b k; wj   b a; b; c with respect to a; b; c equal to zero, respectively, one has derivative of lnp b k; wj

ai ¼

b¼

2

h

1

hi þ ðWþbH þ cPÞ1

i

ð24Þ ii

Nm

2 i Pn h PNm bk r kr ðhÞ þ i¼1 ðWþbH þ cPÞ1 H r¼1 bk r ii

ð25Þ

N N

m p c¼P P h i  i2 Pn h Np Nm b þ i¼1 ðWþbH þ cPÞ1 P r¼1 k¼1 / k;r  /k;r h

ð26Þ ii

The above derivation is detailed in Appendix A. Note that the solutions of h (Eq. (15)) and fa; b; cg are coupled. Therefore, they can be determined from an iterative process of Eqs. (15) and (24)–(26). The above formulations can be summarized below: Algorithm based on Laplace Approximation 1. 2. 3. 4. 5. 6.

Initialize að0Þ ; bð0Þ ; cð0Þ and hð0Þ ; For j ¼ 1; 2;    , update hð jÞ using Eq. (15) with að j1Þ ; bð j1Þ ; cð j1Þ ; Update ai ð jÞ using Eq. (24) with hð jÞ ; bð j1Þ ; cð j1Þ ; Update bð jÞ using Eq. (25) with hð jÞ ; ai ð jÞ ; cð j1Þ ; Update cð jÞ using Eq. (26) with hð jÞ ; ai ð jÞ ; bð jÞ ; Let j ¼ j þ 1, repeat Steps 2 ~ 5 until the convergence criterion is satisfied, that is, khðjÞ  hðj1Þ k2 =khðjÞ k2  Tol.

4. An experimental example 4.1. Model description An experimental three-story steel frame [38] is employed here to validate the proposed algorithm. The structure shown in Fig. 1 was mounted on the strong floor through a steel plate. Therefore the boundary condition can be considered as the fixed constraints. The height of the frame is 1.5 m with each story 0.5 m high, and the span also equals 0.5 m. The cross section of the beams and columns is 75.0  5.0 mm2. The Young’s modulus of the steel is 2.0  1011N/m2, and the mass density is 7.92  103 kg/m3. In this study, the frame is divided into 225 Euler-Bernoulli beam elements, each with a length of 20 mm. A total of 39 measurement points was arranged every 100 mm to obtain the mode shapes of the whole structure. Accelerometers were used to record the acceleration signals induced by an instrument hammer with a rubber tip. The first seven frequencies and mode shapes were extracted using the rational fraction polynomial method [39]. Then three cuts in different crucial locations were sequentially introduced to the structure, as shown in Fig. 1. Cuts 1–3 were located at elements 1, 176, and 12, with the same length b = 20 mm and different depths of d = 22.5, 22.5 and 30 mm, respectively. As the length of the cuts is identical to that of the beam elements, the damage severity, quantified by the SRF (Eq. (3)), equals to the reduction in the moment of inertia of the cross section. That is, the SRFs of the three damaged elements are 60%, 60% and 80%, respectively. Table 1 lists the detailed information of the three damage scenarios (DSs). The number of damaged elements in all three DSs is extremely less than the total number of elements. Therefore, the SRF vector is very sparse with several nonzero items at the damaged elements whereas most zero items at others. The modal tests mentioned previously were conducted in each DS and the corresponding modal data are compared in Table 2. It can be found that the natural frequencies of the frame decrease with the increase of damage severity. The modal assurance criterion (MAC) is employed for mode matching such that the measured modes match the corresponding ones in the finite element model [38]. 4.2. Damage identification The iterative process requires initializing the damage index and hyper-parameters first. Modal experiment experiences show that natural frequencies may comprise 1% noise and mode shapes are usually not as accurate as the frequencies [40]. Therefore, a noise level of 1% and 5% are assigned to the frequencies and mode shapes, respectively, that is,

7

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

Fig. 1. Configuration of the experimental frame (unit: mm).

Table 1 Damage locations and severities for three damage scenarios. Damage Scenario

Cut No.

Element No.

Cut depth (mm)

Damaged component

SRF (h)

DS1 DS2

Cut Cut Cut Cut Cut Cut

1 1 176 1 176 12

22.5 22.5 22.5 22.5 22.5 30

Column Column Beam Column Beam Column

h1 ¼ 60% h1 ¼ 60% h176 ¼ 60% h1 ¼ 60% h176 ¼ 60% h12 ¼ 80%

DS3

1 1 2 1 2 3

Table 2 Frequencies of the frame in the undamaged and damaged states (units: Hz). Mode

Undamaged

DS1

DS2

DS3

1 2 3 4 5 6 7 Average (%)

4.23 14.03 25.45 44.81 58.12 68.36 72.27

4.13 (2.31) 13.75(1.96) 25.14(1.19) 44.70(0.23) 57.39(1.24) 67.34(1.49) 72.06(0.28) (1.44)

4.08 (3.53) 13.45(4.11) 25.13(1.23) 44.69(0.27) 57.28(1.44) 66.11(3.29) 71.42(1.18) (2.32)

4.06 (3.88) 13.42(4.32) 25.09(1.41) 44.62(0.42) 56.55(2.69) 65.31(4.46) 70.74(2.12) (3.37)

Note: Values in parentheses are the frequency change ratios (%) between the damaged and undamaged states.

bð0Þ ¼ 1 = ð0:01Þ2 ¼ 1  104 and cð0Þ ¼ 1 = ð0:05Þ2 ¼ 400. The noise level of the damage index in the prior is estimated as ð0Þ

10% of the exact value, namely, ai

¼ 1 = ð10%Þ2 ¼ 100 ði ¼ 1; 2; . . . ; 225Þ. The damage location and severity are unknown

beforehand, thus the initial damage indexes are assumed as 0, that is, hð0Þ ¼ f0; :::; 0gT . The convergence criterion for the iterðj1Þ k = k h ð jÞ k 6 0:01. ðjÞ  h ation process is set as k h 2

2

8

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

The proposed algorithm is applied to detect damage in the three DSs. An identification error d is defined to quantify the damage severity shown as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kh  hk2 d¼ n

ð27Þ

where h is the identified damage index. In DS1, only one element was damaged. The iterative identification results are shown in Fig. 2. For brevity, only the results in the first two and the last two iterations are shown. In the first two iterations, several elements are falsely identified as damaged. After six iterations, the results converge and the actual sparse damage is accurately detected. The identification error upon convergence is 0.37%. In DS2, elements 1 and 176 were damaged. The convergence is achieved after five iterations. The results in the first two and last two iterations are shown in Fig. 3. Upon convergence, both two damaged elements are accurately detected with the identification error of 1.29%. Three elements (1, 12 and 176) were damaged in the last damage scenario (DS3). The convergence is achieved after four iterations and the damage identification results in each iteration are shown in Fig. 4. Again all three damages are accurately detected with the identification error of 1.36%.

Fig. 2. Damage identification results of DS1.

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

9

Fig. 3. Damage identification results of DS2.

It can be seen that all damages in the three DSs are located accurately and no false identification occurs. All the identification errors are less than 2%. Therefore, the proposed algorithm is able to detect the damage location and severity accurately. 4.3. Variation of hyper-parameters The hyper-parameters alter during the iteration process. Fig. 5 shows the variation of some hyper-parameters of DS3. Although ai of all elements are initialized identically as 100, their values upon convergence vary remarkably. Specifically, ai associated with the undamaged elements (a100 ; a150 and a225 Þ converge to a larger value in comparision with ai corresponding to the damaged elements (a1 ; a12 and a176 ). This feature is consistent with the sparse mechanism of the ARD model, in which each variable is assigned with an individual hyper-parameter. In particular, ai corresponding to the undamaged elements become significantly large. Therefore, the associated items are penalized more in the optimization (see Eq. (15)), forcing these damage indexes hi to zero and realizing the sparse damage detection. 5. Comparison with EM algorithm We proposed a new structural damage detection method using the expectation–maximization (EM) technique [12], in which the evidence is maximized through expectation and maximization steps. In short, in the E-step the expectation of the logarithm of the complete-data is solved as

10

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

Fig. 4. Damage identification results of DS3.

8 " #2 9 ( )   Nm Nm = N N
 n n   1 1X 1X þ ln ai  ai E h2i 2p 2 i¼1 2 i¼1

ð28Þ

In the M-step, to maximize the expectation, Eq. (28) is differentiated with respect to a, b, and c, and set to zero, that is

n  o b @E ln p h; njb k; w @ ai n  o b @E ln p h; njb k; w @b n  o b @E ln p h; njb k; w @c

¼

1 1    E h2i ¼ 0 2ai 2

8 " #2 9 Nm b k r  kr ðhÞ = Nm 1 1
ð29Þ

ð30Þ

ð31Þ

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

11

Fig. 5. Variation of hyper-parameters during the iterative process of DS3.

The hyper-parameters are then solved as

1

ai ¼  2  E hi

ð32Þ

12

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

Nm b¼ (

 ) PNm bk r kr ðhÞ 2 E r¼1 bk r

c ¼ nP E

ð33Þ

Np  Nm

Nm b r¼1 k / r

2

 /r ðhÞk2

ð34Þ

o

The damage index h is similarly obtained by minimizing the objective function Eq. (15). Given the expectations in Eqs. (32)–(34) cannot be calculated directly, the numerical sampling technique is required. The posterior and likelihood samplings were employed in Ref. [12]. By contrast, in the present study the solutions of hyper-parameters are expressed in a closed form as shown in Eqs. (24)– (26). The mean and variance of the posterior PDF are obtained based on the Laplace approximation. Consequently, the damage index and hyper-parameters are solved directly without sampling. To compare these two algorithms, the numerical cantilever beam studied in Ref. [12] is utilized here, as shown in Fig. 6. The overall length of the beam is 0.9 m. The cross section is 50.75  6.0 mm2 and the mass density is 7.67  103 kg/m3. The Young’s Modulus is estimated as 2.0  1011N/m2. The beam is divided into 45 equal Euler-Bernoulli beam elements (i.e., n = 45), each with a length of 20 mm. Two cuts are located in element 1 and element 23, i.e., h1 ¼ h23 ¼ 0:5. The first six natural frequencies of the beam in the undamaged and damaged conditions as listed in Table 3 are used for damage identification. ð0Þ

The initializations of the hyper-parameters and damage index are the same as above, that is, ai ¼1=ð10%Þ2 ¼ 100 , bð0Þ ¼ 1=ð0:01Þ2 ¼ 1  104 and hð0Þ ¼ f0; :::; 0gT . Given that only frequencies are taken into calculation, the items associated with mode shapes are dropped. The damage identification results corresponding to each iteration step using the Laplace approximation are shown in Fig. 7. It can be seen that the damaged elements are accurately detected. The identification error and computation time are compared with those using EM technique in Table 4. Both calculations are carried out on a PC with Intel Core i7 3.20 GHz CPU and 4 GB RAM. The damage identification errors of both techniques are very small. The computational time of the proposed Laplace approximation is much lesser than the EM technique as no sampling is required in the proposed method. Consequently, the Laplace approximation is very efficient.

Fig. 6. Geometric configuration of the beam structure (unit: mm).

Table 3 Frequencies of the beam in the undamaged and damaged states. Mode no.

Undamaged

Damaged

Freq. (Hz)

Freq. (Hz)

1 6.02 2 37.75 3 105.73 4 207.25 5 342.70 6 512.07 Average of frequency change (%)

5.75 35.67 102.44 197.69 333.96 492.45

Change ratio (%)

4.56 5.50 3.11 4.61 2.55 3.83 4.03

13

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

Fig. 7. Damage detection results of the beam using the Laplace approximation.

Table 4 Comparison of damage detection results using EM and Laplace approximation. Algorithms

Iteration No.

Identification error

Time consumption

EM (Likelihood sampling) Laplace approximation

4 3

0.41% 0.37%

261.7 min 3.5 min

6. Conclusions The paper proposed a sparse Bayesian model for structural damage detection based on the Laplace approximation. The ARD model with Gaussian distributions is incorporated as the sparsity induced prior in this framework, which equivalently serving as the regularization item to tackle the ill-posed inverse problem. In view of the intractable posterior PDF caused by the nonlinear relationship between the damage index and modal parameters, the Laplace approximation is employed to obtain the asymptotic solution of the evidence. In this manner, the solutions of the damage index and hyper-parameters are expressed in an analytical form and solved iteratively. An experimental frame is utilized to verify the effectiveness of the proposed algorithm. All damages in the three DSs can be correctly located and quantified without false identification. Moreover, compared with the EM technique, the merit of this proposed algorithm lies in the avoidance of sampling, which makes the computation efficient. Acknowledgements The research in this paper was supported by the PolyU Research Grant (Project No. 1-ZVP2) and the National Natural Science Foundation of China (Project No. 51678364).

14

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

Appendix A Hyper-parameters a; b and c are solved as follows.   b a; b; c in Eq. (23) with respect to a as zero, we have Setting the derivative of lnp b k; wj

  b a; b; c @lnp b k; wj @ ai

 1 hi 2 1 @ln jWþbH þ cPj 1 hi 2 1 @ ðWþbH þ cPÞ  ¼    Tr ðWþbH þ cPÞ1 2ai 2 @ ai 2ai @ ai 2 2 2 i 1 hi 2 1 h  ¼  ðWþbH þ cPÞ1 ¼ 0 ii 2ai 2 2

¼

ðA:1Þ

where Tr½ represents the trace of a matrix which is defined as the sum of elements on the leading diagonal. In the above equation, the following derivative equation of a matrix is used:

  @lnjAj @A ¼ Tr A1 @a @a

ðA:2Þ

cPÞ As W is a diagonal matrix with Wii ¼ ai , @ ðWþbHþ is a diagonal zero matrix except a unit at the ith element. From Eq. (A.1), @ ai

we have

ai ¼

h

1

hi þ ðWþbH þ cPÞ1 2

i

ðA:3Þ ii

  b a; b; c with respect to b is obtained below. Similarly, the derivative of lnp b k; wj

  b a; b; c @lnp b k; wj @b

" Nm b kr Nm 1 X  2b 2 r¼1 " Nm b kr Nm 1 X ¼  2b 2 r¼1 " Nm b kr Nm 1 X ¼  2b 2 r¼1

¼

 #2  kr h 1 @ln jWþbH þ cPj  b 2 @b kr  #2 i  kr h 1 h  Tr ðWþbH þ cPÞ1 H b 2 kr  #2 n h i  kr h 1X  ðWþbH þ cPÞ1 H b ii 2 i¼1 kr

ðA:4Þ

Setting Eq. (A.4) as zero, b is obtained as

b¼

Nm

 i PNm bk r kr ðhÞ 2 Pn h þ i¼1 ðWþbH þ cPÞ1 H r¼1 bk r ii

ðA:5Þ

  b a; b; c with respect to c is obtained below. The derivative of lnp b k; wj

  b a; b; c @lnp b k; wj @c

¼

Np h Nm X  i2 1 @ln jWþbH þ cPj Nm Np 1 X b / h /   k;r k;r 2 r¼1 k¼1 2 @c 2c

¼

Np h Nm X i  i2 1 h Nm Np 1 X b / h /   Tr ðWþbH þ cPÞ1 P k;r 2 r¼1 k¼1 k;r 2 2c

¼

Np h Nm X n h i  i2 1 X Nm Np 1 X b / h /   ðWþbH þ cPÞ1 P k;r k;r ii 2 r¼1 k¼1 2 i¼1 2c

ðA:6Þ

Setting Eq. (A.6) as zero, c is obtained as

N N

m p c¼P P h i  i2 Pn h Np Nm b þ i¼1 ðWþbH þ cPÞ1 P r¼1 k¼1 / k;r  /k;r h

ðA:7Þ ii

References [1] S.W. Doebling, C.R. Farrar, M.B. Prime, D.W. Shevitz, Damage identification and health monitoring of structural and mechanical systems from changes in their vibration characteristics: a literature review, Los Alamos National Laboratory Report, 1996.

X. Wang et al. / Mechanical Systems and Signal Processing 140 (2020) 106701

15

[2] H. Sohn, C.R. Farrar, F.M. Hemez, D.D. Shunk, D.W. Stinemates, B.R. Nadler, J.J. Czarnecki, A review of structural health monitoring literature: 19962001, Los Alamos National Laboratory Report, 2003. [3] Y.J. Yan, L. Cheng, Z.Y. Wu, L.H. Yam, Development in vibration-based structural damage detection technique, Mechan. Syst. Signal Processing 21 (2007) 2198–2211. [4] E.P. Carden, P. Fanning, Vibration based condition monitoring: a review, Struct. Health Monitor. 3 (2004) 355–377. [5] F.M. Hemez, Uncertainty quantification and the verification and validation of computational models, Damage Prognosis for Aerospace, Civil and Mechanical Systems, 2004, 201–220. [6] J. Zhang, Y.L. Xu, Y. Xia, J. Li, A new statistical moment-based structural damage detection method, Struct. Eng. Mechan. 30 (2008) 445–466. [7] Y.L. Xu, J. Zhang, J.C. Li, Y. Xia, Experimental investigation on statistical moment-based structural damage detection method, Struct. Health Monitor. 8 (2009) 555–571. [8] Y.Q. Bao, Z.C. Chen, S.Y. Wei, Y. Xu, Z.Y. Tang, H. Li, The state of the art of data science and engineering in structural health monitoring, Engineering 5 (2019) 234–242. [9] O.S. Salawu, Detection of structural damage through changes in frequency: a review, Eng. Struct. 19 (1997) 718–723. [10] Y. Xia, H. Hao, G. Zanardo, A.J. Deeks, Long term vibration monitoring of an RC slab: temperature and humidity effect, Eng. Struct. 28 (2006) 441–452. [11] X.Q. Zhou, Y. Xia, S. Weng, l1 regularization approach to structural damage detection using frequency data, Struct. Health Monitor. 14 (2015) 571–582. [12] R.R. Hou, Y. Xia, X.Q. Zhou, Y. Huang, Sparse Bayesian learning for structural damage detection using expectation-maximization technique, Struct. Control Health Monitor. 26 (2019), e2343. [13] M.E. Tipping, Sparse Bayesian learning and the relevance vector machine, J. Mach. Learn. Res. 1 (2001) 211–244. [14] J.L. Beck, Bayesian system identification based on probability logic, Struct. Control Health Monitor. 17 (2010) 825–847. [15] J.L. Beck, L.S. Katafygiotis, Updating models and their uncertainties. I: Bayesian statistical framework, J. Eng. Mech. 124 (1998) 455–461. [16] L.S. Katafygiotis, J.L. Beck, Updating models and their uncertainties. II: Model identifiability, J. Eng. Mech. 124 (1998) 463–467. [17] Y. Huang, J.L. Beck, H. Li, Bayesian system identification based on hierarchical sparse Bayesian learning and Gibbs sampling with application to structural damage assessment, Comput. Methods Appl. Mech. Eng. 318 (2017) 382–411. [18] Y. Huang, J.L. Beck, Hierarchical sparse Bayesian learning for structural health monitoring with incomplete modal data, Int. J. Uncertain. Quantif. 5 (2015) 139–169. [19] Y. Huang, J.L. Beck, H. Li, Hierarchical sparse Bayesian learning for structural damage detection: theory, computation and application, Struct. Saf. 64 (2017) 37–53. [20] M.J. Beal, Z. Ghahramani, The variational Bayesian EM algorithm for incomplete data: with application to scoring graphical model structures, Bayesian Statist. 7 (2003) 453–463. [21] M.J. Beal, Variational algorithms for approximate Bayesian inference. PhD thesis, University of London, 2003. [22] S.D. Babacan, R. Molina, A.K. Katsaggelos, Bayesian compressive sensing using Laplace priors, IEEE Trans. Image Process. 19 (2010) 53–63. [23] K.V. Yuen, J.L. Beck, L.S. Katafygiotis, Efficient model updating and health monitoring methodology using incomplete modal data without mode matching, Struct. Control Health Monitor. 13 (2006) 91–107. [24] J.N. Rouder, J. Lu, An introduction to Bayesian hierarchical models with an application in the theory of signal detection, Psychon. Bull. Rev. 12 (2005) 573–604. [25] I. Behmanesh, B. Moaveni, G. Lombaert, C. Papadimitriou, Hierarchical Bayesian model updating for structural identification, Mech. Syst. Sig. Process. 64 (2015) 360–376. [26] O. Sedehi, C. Papadimitriou, L.S. Katafygiotis, Probabilistic hierarchical Bayesian framework for time-domain model updating and robust predictions, Mech. Syst. Sig. Process. 123 (2019) 648–673. [27] M. Song, B. Moaveni, C. Papadimitriou, A. Stavridis, Accounting for amplitude of excitation in model updating through a hierarchical Bayesian approach: application to a two-story reinforced concrete building, Mech. Syst. Sig. Process. 123 (2019) 68–83. [28] H.F. Lam, J.H. Yang, S.K. Au, Bayesian model updating of a coupled-slab system using field test data utilizing an enhanced Markov chain Monte Carlo simulation algorithm, Eng. Struct. 102 (2015) 144–155. [29] H.F. Lam, J.H. Yang, S.K. Au, Markov chain Monte Carlo-based Bayesian method for structural model updating and damage detection, Struct. Control Health Monitor. 25 (2018), e2140. [30] Q. Long, M. Scavino, R. Tempone, S.J. Wang, Fast estimation of expected information gains for Bayesian experimental designs based on Laplace approximations, Comput. Methods Appl. Mech. Eng. 259 (2013) 24–39. [31] P.S. Laplace, Memoir on the probability of the causes of events, Statist. Sci. 1 (1986) 364–378. [32] A. Onar, Laplace Approximations in Bayesian lifetime analysis, Encyclop. Statist. Quality Reliab. 2 (2008). [33] L. Tierney, J.B. Kadane, Accurate approximations for posterior moments and marginal densities, J. Am. Stat. Assoc. 81 (1986) 82–86. [34] C.M. Bishop, Pattern recognition and machine learning, Springer, Berlin, 2006. [35] R.B. Nelson, Simplified calculation of eigenvector derivatives, AIAA J. 14 (1976) 1201–1205. [36] S. Weng, Y. Xia, Y.L. Xu, H.P. Zhu, An iterative substructuring approach to the calculation of eigensolution and eigensensitivity, J. Sound Vib. 330 (2011) 3368–3380. [37] S. Weng, H.P. Zhu, Y. Xia, X.Q. Zhou, L. Mao, Substructuring approach to the calculation of higher-order eigensensitivity, Comput. Struct. 117 (2013) 23–33. [38] R.R. Hou, Y. Xia, Q. Xia, X.Q. Zhou, Genetic algorithm based optimal sensor placement for l1 regularized damage detection, Struct. Control Health Monitor. 26 (2019). [39] D. Formenti, M. Richardson, Parameter estimation from frequency response measurements using rational fraction polynomials (twenty years of progress), Proceeding of International Modal Analysis Conference, 4753, 2002, 373–382. [40] J.E. Mottershead, M.I. Friswell, Model updating in structural dynamics: a survey, J. Sound Vib. 167 (1993) 347–375.