- Email: [email protected]
Hierarchical sparse Bayesian learning for structural damage detection: Theory, computation and application
Structural Safety 64 (2017) 37–53
Contents lists available at ScienceDirect
Structural Safety journal homepage: www.elsevier.com/locate/strusafe
Hierarchical sparse Bayesian learning for structural damage detection: Theory, computation and application Yong Huang a,b, James L. Beck b,⇑, Hui Li a a b
Key Lab of Structural Dynamic Behavior and Control of the Ministry of Education, School of Civil Engineering, Harbin Institute of Technology, Harbin, China Division of Engineering and Applied Science, California Institute of Technology, CA, USA
a r t i c l e
i n f o
Article history: Received 10 March 2016 Received in revised form 2 September 2016 Accepted 11 September 2016 Available online 31 October 2016 Keywords: Inverse problems Hierarchical sparse Bayesian learning Model updating Structural health monitoring Damage detection Bayesian Ockham razor IASC-ASCE benchmark problems
a b s t r a c t Structural damage due to excessive loading or environmental degradation typically occurs in localized areas (in the absence of collapse) where it leads to local stiffness reductions. This prior information about the spatial sparseness of structural damage and the associated stiffness loss is exploited here by a hierarchical sparse Bayesian learning (SBL) framework, with the goal to reduce the ill-conditioning in the stiffness loss inversion problem for damage detection. We have previously proposed a SBL approach to establish the probability of localized stiffness reductions caused by damage by using noisy incomplete modal data from before and after possible damage. The excellent performance achieved by introducing sparseness in the damage pattern was demonstrated by using synthetic data where there are only small modeling errors. In this research, a more rigorous formulation along with a corresponding efficient and scale-invariant SBL algorithm are developed. The algorithm is first applied to synthetic data, then to real vibration response from a steel-frame test structure where there are large modeling errors. These data are from the Phase II simulated and experimental benchmark studies that were sponsored by the IASC-ASCE Task Group on Structural Health Monitoring. The results show that, even for the real data, the proposed method can reliably detect, locate and assess damage of the benchmark structure by inferring substructure stiffness losses using the identified modal parameters from the calibration and monitoring stages. The occurrence of missed and false damage alerts is effectively suppressed, and we show that the new algorithm gives a better performance than our previous SBL method in the real data case where there is significant modeling error. Several appealing features of our method are summarized at the end of the paper. Ó 2016 Published by Elsevier Ltd.
1. Introduction The challenge of developing automated sensor-based structural health monitoring (SHM) systems for accurately detecting, locating and assessing damage from severe loading events (such as earthquakes, hurricanes or explosions) or progressive structural deterioration (such as fatigue) has received much attention over the last several decades. Numerous vibration-based global SHM techniques have been developed for estimating both damage location and damage severity [e.g. 1–9], most of which utilize the dependence of the identified structural modal parameters, such as natural frequencies and mode shapes, on physical properties of structures, i.e., stiffness, mass and damping. In practice, most SHM systems record structural vibration data at low-amplitude excitation, e.g., ambient vibration from wind and ground motion, and so a linear ⇑ Corresponding author. E-mail address: [email protected] (J.L. Beck). http://dx.doi.org/10.1016/j.strusafe.2016.09.001 0167-4730/Ó 2016 Published by Elsevier Ltd.
dynamic model with classical normal modes is adequate for damage detection purposes. The typical strategy involves comparing structural models identified from sets of measured modal data (natural frequencies and mode shapes) from a structure before and after possible damage. For applications in structural health monitoring, there are some inherent difficulties in this modelbased approach. One of the main difficulties is that it is impossible to exactly model the full behavior of a structure by using the limited knowledge available. For damage inference problems, there are always modeling uncertainties involved and hence the need to account for uncertainty in model updating arises. These uncertainties may arise from many sources: simplifying approximations to develop the structural models; limited number of sensors installed on a structure; sensor noise; thermally-induced daily variations in structural stiffness, etc. This presence of substantial modeling uncertainty has motivated numerous researchers to tackle the problem of model updating in structural system identification
38
Y. Huang et al. / Structural Safety 64 (2017) 37–53
Nomenclature Nm No Nd Ns Nh M, K h /r ; x2r b ^ r;i ^ 2r;i ; w x
number of extracted modes in the modal identification we use (i = 1, . . ., Nm) number of observed degrees of freedom number of degrees of freedom of the structural model for identification we use (k = 1, . . ., Nd) number of time segments of measured modal data we use (r = 1, . . ., Ns) number of substructures considered (j = 1, . . ., Nh) mass and stiffness matrices of the structural model for identification structural stiffness scaling parameters system mode shape and system natural frequencies of the rth mode equation error precision parameter MAP estimates of modal frequency and mode shape of the ith mode from the rth data segment from modal identification
and damage assessment from a Bayesian perspective [2,3,10–26]. Rather than considering only a point estimate as in the conventional deterministic or frequentist probabilistic methods, Bayesian methods use Bayes’ Theorem to quantify the relative plausibility of all possible values of the structural model parameters via their Bayesian posterior PDF (probability density function). The relative plausibility of each model within a set of candidate model classes chosen to represent the uncertain structural behavior can also be assessed by applying Bayes’ Theorem at the model class level, enabling the consideration of all plausible models and the construction of robust predictions of future structural behavior [21]. We emphasize that probability in the Bayesian sense should be interpreted as the degree of plausibility of a statement on the basis of the specific conditioning information [21,27–29], where the statement may refer to the structural parameters or to the structural model itself. It quantifies uncertainty due to our incomplete information and does not need the postulate of ‘‘inherent randomness”. Another fundamental difficulty for damage detection lies in the fact that the mathematical nature of the underdetermined structural inverse problem based on the constraints imposed by the available data requires a tradeoff between the resolution of the damage locations and the reliability of the probabilisticallyinferred damaged state. In this case, defining a proper threshold to determine whether the damage features shift from their healthy state is very important in order to alleviate false positive (false alarm) and false negative (missed alarm) detections, but it is very challenging to establish a reliable threshold value in a rigorous manner for issuing a timely damage alarm [30,31]. In this paper, we aim to extend the applicability of Bayesian methods to produce reliable damage alarms even for highdimensional model parameter spaces (higher-resolution damage localization), by the exploration of recent developments in sparse Bayesian learning (SBL) [32,33] and Bayesian compressive sensing/sampling [34–36]. This new perspective for the stiffness loss inverse problem exploits the prior knowledge that structural damage typically occurs at a limited number of locations in a structure in the absence of its collapse. A specific hierarchical Bayesian model and the corresponding evidence maximization strategy for Bayesian inference is utilized to promote spatial sparseness in the inferred stiffness reductions in a way that is consistent with the Bayesian Ockham razor [21,37]. We have proposed previously to use a sparse Bayesian perspective to tackle the stiffness loss inverse problem [38,39]. In this
^hu
C T
a g, q m, s, j d
MAP estimate of h determined from the calibration test data. matrix that picks the measured degrees of freedom from each system mode shape matrix relating the vector of Ns sets of identified natural ^ 2 and the system natural frequencies x2 frequencies x variance parameter for the likelihood function of structural stiffness scaling parameters h measurement-error precision parameters for mode shapes and natural frequencies rate parameters controlling the exponential prior distributions of g, q, b0 set of all hyper-parameters x2, q, s, /, g, m, f, b, a0, b0, j
paper, we present an improved version of our SBL method that eliminates two possibly unreliable approximations in our original theoretical formulation. We also propose effective strategies to alleviate the ill-conditioning for stiffness inversions when the modeling error is large, which enhance its applicability and accuracy for damage detection and assessment with real data. In addition, the new algorithm allows more efficient computation and is independent of the unit selections of mass and stiffness matrices, as well as the scaling of the measured mode shapes. We also apply our method to analyze the IASC–ASCE Phase II experimental benchmark data for both brace and joint damage cases, not only the synthetic data that we used previously [38,39]. For both the simulated and experimental data, we show that the performance of the new method is better than our previous method in [39], especially for the joint damage cases. The remainder of the paper is organized as follows. In Section 2, we present our hierarchical Bayesian model, which is similar to that in [39] but which has several modifications. The Bayesian model is then used in an evidence maximization procedure in Section 3, yielding an improved Bayesian inference framework that lends itself to reliable, efficient and scale-invariant computation of the posterior PDF for the structural stiffness parameters. In Section 4, the algorithms for stiffness identification in the calibration and monitoring stages are presented. We then compare in Section 5 our new SBL damage assessment algorithm with that in [39] in terms of the theory and computation. Applications of our Bayesian methods to simulated and real structural data are presented in Section 6, and conclusions and future work are discussed in Section 7. 2. Hierarchical Bayesian model class Suppose that Ns sets of measured vibration time histories are available from a structure and Nm modes of vibration have been identified for each set of time histories so that we have a vector of h iT ^2¼ x ^ 21;1 ; ...; x ^ 21;Nm ; x ^ 22;1 ; ... ; x ^ 2Ns ;Nm identified natural frequencies: x 2 RNs Nm 1 and a vector of identified mode shapes: h iT N N N 1 T T T T s m o ^¼ w ^ ;...;w ^ ^ ^ ^ r;i 2 RNo w 2R , where w 1;1 1;Nm ; w2;1 ; . . . ; wNs ;Nm gives the identified components of the system mode shape of the ith mode (i = 1, . . ., Nm) at the No measured DOF (degrees of freedom) from the rth data segment (r = 1, . . ., Ns). The No measured DOF in a structure are usually a smaller subset of the Nd DOF of an appropriate structural model. It has been
39
Y. Huang et al. / Structural Safety 64 (2017) 37–53
shown to be advantageous [3,10,12–15] to introduce the system T frequencies x2 ¼ x21 ; . . . ; x2Nm 2 RNm 1 and system T T mode shapes / ¼ /1 ; . . . ; /TNm 2 RNd Nm 1 to represent the actual underlying modal parameters of the linear dynamics of the structural system at all DOFs corresponding to those of the structural model. The system modal parameters are therefore distinct from the corresponding modal parameters of the structural model. Based on a defined stochastic model class M (see ^ to update ^ 2 and w [21]), one can use the available modal data x
natural
the system modal parameters x2 and / together with structural model parameters h, without using mode matching [2,3,10,12– 15]. 2.1. Structural model class
As part of the definition of the stochastic model class M, we choose a set of linear structural models with Nd DOFs and with classical normal modes (so a damping matrix need not be explicitly modeled) where each model has the same known mass matrixM 2 RNd Nd and an uncertain stiffness matrix K parameterized as a linear combination of (Nh + 1) substructure stiffness matrices Kj, j = 0, 1, . . ., Nh: Nh X KðhÞ ¼ K0 þ hj Kj
ð1Þ
j¼1
where K0 is any part of the global stiffness matrix that is considered known and that corresponds to parts of the structure that are believed to be not vulnerable to damage, and Kj 2 RNd Nd ; j ¼ 1; . . . ; N h ; is the substructure stiffness matrix representing the nominal contribution of the jth substructure to the global stiffness matrix K from a structural model (e.g., specified by a finite element model of a structure). The corresponding stiffness scaling parameter hj, j = 1, . . ., Nh, is a factor that allows modification of the nominal jth substructure stiffness so as to be more consistent with the real structure behavior and any reduction of hj, j = 1, . . .Nh, is assumed to correspond to damage in the jth substructure. The parameter vectorh ¼ ½h1 ; . . . ; hn 2 RNh ; therefore specifies the structural model in M: 2.2. Prior distribution for structural model parameters and system modal parameters We do not assume that the Nm system natural frequencies x2 and Nm system mode shapes / satisfy the eigenvalue problem corresponding to any structural model specified by h because there will always be modeling errors e ¼ eT1 ; . . . ; eTNm 2 RNm Nd ; where ðKðhÞ x2i MÞ/i ¼ ei ; i ¼ 1; . . . ; N m ; for the ith system mode and the structural model specified by h. The uncertain eigenvalue equation error e therefore provides a bridge between the behaviors of the real system and a deterministic structural model. As described in [39], the following prior PDF conditional on the equation-error precision hyper-parameter b is chosen:
pðx2 ; /; hjbÞ ¼ c0 ð2p=bÞNm Nd =2 ( ) Nm bX 2 2 exp kðKðhÞ xi MÞ/i k 2 i¼1
can deduce the prior PDF for h conditional on system modal parameters [x2, /] as: 1
1
pðhjx2 ;/;bÞ ¼ pðx2 ;/;hjbÞ=pðx2 ;/jbÞ ¼ NðhjðHT HÞ HT b;ðbHT HÞ Þ ¼ ð2p=bÞNh =2 jHT Hj1=2 T 1 1 b exp ðh ðHT HÞ HT bÞ HT Hðh ðHT HÞ HT bÞ 2 ð3Þ where the marginal prior PDF for the system modal parameters [x2, /] and the parameter vector b and matrix H are defined by:
pðx2 ; /jbÞ ¼ c0 ð2p=bÞðNh Nd Nm Þ=2 jHT Hj1=2 1 b T T exp ðb b b HðHT HÞ HT bÞ 2 3 ðx21 M K0 Þ/1 6 7 .. 7 b¼6 . 4 5 ðx2Nm M K0 Þ/Nm N
ð4Þ
2
2
K1 /1 6 . H¼6 4 .. K1 /Nm
ð5Þ m N d 1
3 KNh /1 7 .. 7 . 5 KNh /Nm Nm N .. .
ð6Þ d N h
We model the uncertainty in b by the following Exponential hyper-prior: 1
1
pðbjb0 Þ ¼ b0 expðb0 bÞ
ð7Þ
where b0 is the scale parameter of the distribution. The Exponential prior for b is the maximum entropy distribution with the constraint of E(b|b0) = b0 and is a special case of the widely used Gamma prior. To penalize values of b0 that are too large, we assign an Exponential hyper-prior for the scale parameter b0 as the last stage of the hierarchical model:
pðb0 jjÞ ¼ Expðb0 jjÞ ¼ j expðjb0 Þ
ð8Þ
Remark 2.1. The joint prior in (2) is defined based on the Gaussian probability model for uncertain equation errors e, which is the maximum entropy PDF subject to the first two moments as constraints. Those values of {x2, /, h} that make the equation errors e smaller in norm are considered more plausible a priori; i.e. they have larger values for their prior PDF. However, the joint PDF in (2) has an undesirable peak where all the components of the system mode shapes / are zeros since the corresponding equation errors e are then also zero. For real damage detection when the model is expected to have large modeling errors, the scale parameter b0 is important for scaling the equation-error precision parameter b appropriately to prevent equation error e, and so the components of /, becoming zero. This point is discussed further at the end of the Appendix. 2.3. Likelihood function for structural model parameters and system modal parameters
ð2Þ
where c0 is a normalizing constant and |||| denotes the Euclidean vector norm. Note that the finite value of the equation-error precision parameter b in (2) provides a soft constraint for the eigenequation and it allows for the explicit control of how closely the system and model modal parameters agree. Following [39], we
Since the onset of damage is typically in a small number of locations in a structure in the absence of structural collapse, the poten^u can tial change in the parameter vector from damage Dh ¼ h h be considered as a sparse vector with relatively few non-zero components, where h and hu are the stiffness scaling parameters for the current state (possibly damaged) and the calibration (undamaged) state, and ^ hu is the unique MAP estimate of hu from applying Baye-
40
Y. Huang et al. / Structural Safety 64 (2017) 37–53
sian updating using a large amount of time-domain vibration data from the undamaged structure. Following [39], we choose the MAP value ^ hu from the calibration state as pseudo-data for h to define a likelihood function as:
pð^hu jh; aÞ ¼ Nð^hu jh; AÞ ¼
Nh Y Nð^hu;j jhj ; aj Þ
ð9Þ
j¼1
where A ¼ diagða1 ; . . . ; aNh Þ. The idea here is to learn each aj from hu;j , which is interpreted the modal data and if aj ? 0, then hj ! ^ as the jth substructure being undamaged. The choice of likelihood function in (9) is motivated by the closely-related SBL framework. Although the conventional strategy in SBL is to use an automatic relevance determination (ARD) Gaussian prior PDF [37] to model sparseness, here we incorporate the ARD concept in the likelihood function along with the prior on h in (3). Following the same approach as in [39], the likelihood function for the system modal parameters x2 and / is:
^ x2 ; /; hÞ ¼ pðx ^ ^ ^ 2 ; wj ^ 2 jx2 Þpðwj/Þ ^ 2 jLx2 ; EÞNðwjC/; pðx ¼ Nðx CÞ ð10Þ where the selection matrixC 2 RNo Nm Ns Nm Nd with ‘‘1s” and ‘‘0s” picks ^ the observed DOF of the whole ‘‘measured” mode shape data set w from the system mode shapes /; L ¼ ½INm ; . . . ; INm T 2 RN0 Nm Nm is the matrix between the vector of Ns sets of Nm identified natural fre^ 2 and the Nm system natural frequencies x2; quencies x E ¼ block diagðE1 ; . . . ; ENs Þ is a block diagonal covariance matrix 1 with the diagonal block Er ¼ diag q1 1 ; . . . ; qN m ; r ¼ 1; . . . ; N s ; C ¼ g1 INo Ns Nm ; INm andINo Ns Nm denote the identity matrices of
corresponding size; and q ¼ ½q1 ; . . . ; qNm T and g are prescribed precision parameters for the predictions of the identified natural fre^ from the system modal ^ 2 and mode shapes w quencies x parameters, respectively. In a hierarchical manner, exponential priors are placed on the parameters qi and g:
pðqi jsi Þ ¼ Expðqi jsi Þ ¼ si expðsi qi Þ; i ¼ 1; . . . ; Nm pðgjmÞ ¼ ExpðgjmÞ ¼ m expðmgÞ
ð11Þ
which are the maximum entropy priors with support [0, 1) for and m1 of qi and g, respectively. Then the given mean values s1 i prior PDF for the parameter vector q is given by:
pðqjsÞ ¼
Nm Nm Nm Y Y X pðqi jsi Þ ¼ si exp si qi i¼1
i¼1
!
ð12Þ
i¼1
where s ¼ ½s1 ; . . . ; sNm . T
2.4. Joint PDF for hierarchical Bayesian model By combining all stages of the hierarchical Bayesian model, the joint PDF for all of the stochastic variables is:
^ ^hu Þ ¼ ^ 2 ; w; pðx2 ; q; s; /; g; m; h; a; b; b0 ; j; x 2 2 ^ ^ ^ jx ; qÞpðwj/; gÞpðhu jh; aÞ pðx pðx2 ; /; hjbÞpðqjsÞpðgjmÞpðaÞpðsÞpðmÞpðbjb0 Þpðb0 jjÞpðjÞ
ð13Þ
^ ^ 2 jx2 ; qÞpðwj/; where PDFs pðx gÞ; pð^hu jh; aÞ; pðx2 ; /; hjbÞ; pðqjsÞ; pðgjmÞ; pðbjb0 Þ and p(b0|j) are defined in (10), (9), (2), (12), (11), (7) and (8), respectively. The prior PDFs for aj, si, m and j are taken as uniform over large open intervals that start at zero. The acyclic graph of the hierarchical Bayesian model is shown in Fig. 1, where each arrow denotes the conditional dependencies used in the joint probability model. The bidirectional arrow in the graph of the hierarchical Bayesian model represents the statistical interaction between the system modal parameters x2 and /, which comes from
Fig. 1. Acyclic graph representing the information flow in the hierarchical Bayesian model.
the joint prior p(x2, /|b) in Eq. (4). Compared with non-hierarchical Bayesian models where b, q and g are selected a priori [15], our hierarchical Bayesian model allows all sources of uncertainty and correlation to be learned from the data, and hence potentially produces more accurate damage identification. Remark 2.2. The hierarchical Bayesian model in Fig. 1 is similar to that of [39]; however, the exponential hyper-prior PDF pðajkÞ for a is replaced by a simple uniform distribution since the learning of this hyper-prior using real-data with large modeling errors may produce over-sparse models with too few stiffness reductions. In addition, an exponential hyper-prior is assigned for scale parameter b0 for appropriately scaling the prior distribution of the equation-error precision parameter b when learning from data where the model is expected to have large modeling errors.
3. Bayesian inference To facilitate the goal for a rigorous SBL formulation, we collect all uncertain parameters except stiffness scaling parameter h in the vector d = [(x2)T, qT, sT, /T, g, m, aT, b, b0, j]T and the full posterior uncertainty in h is explicitly incorporated when finding the MAP estimates of all parameters in d using the evidence maximization procedure. 3.1. Approximation of the full posterior PDF using Laplace’s method As widely known, Bayesian inference is based on the full poste^ ^ ^ 2 ; w; rior PDF over all uncertain parameters pðh; djx hu Þ which is given by Bayes’ theorem, but which is almost always analytically intractable because of the high-dimensional normalization integral. Nevertheless, we can turn to a hierarchical Bayesian procedure combined with Laplace’s asymptotic approximation [40] to find an effective approximation of the full posterior PDF ^ ^ ^ 2 ; w; pðh; djx hu Þ: Using the probability product rule, we rewrite the full posterior as:
^ ^hu Þ ¼ pðhjd; x ^ ^hu Þpðdjx ^ ^hu Þ ^ 2 ; w; ^ 2 ; w; ^ 2 ; w; pðh; djx
ð14Þ
The first factor in (14) is the posterior PDF of the stiffness scaling parameter vector conditional on d, which is given by:
41
Y. Huang et al. / Structural Safety 64 (2017) 37–53
^ ^hu Þ ¼ pðhjd; ^hu Þ ¼ pð^hu jh; dÞpðhjdÞ=pð^hu jdÞ ^ 2 ; w; pðhjd; x
ð15Þ
^ ^ 2 and w where the first equality holds because h is independent of x when conditional on x2 and / in d (see Fig. 1) and the second equality is due to Bayes’ Theorem. As a consequence of combining a Gaussian prior (3) and a linear model within a Gaussian likelihood (9), the conditional posterior is also Gaussian:
^ ^hu Þ ¼ Nðhjl; RÞ ^ 2 ; w; pðhjd; x
ð16Þ
stituting the MAP estimates of m into that of g and letting all other ~ g ~ parameters be fixed at their MAP values, the MAP estimates /; ~ and m are derived as (using (A3), (A5) and (A10) in (A1) in the Appendix):
1 ~T F ~ þ bdiag ~¼ b ~F ~ ~ Nm N Þ þ g ~ 1 Þ; . . . ; trðRT ~ CT C / trðRT d ^ b ~ trðRU ~ 1 Þ; . . . ; trðRU ~ Nm N Þ T ~ CT w g d
ð22Þ
with mean and covariance matrix:
N N N 2
l ¼ RðbHT b þ A1 ^hu Þ
ð17aÞ
1 R ¼ bHT H þ A1
g~ ¼ s^ m o~ 2 ¼ ~ m kw C/k
ð17bÞ
~ a ~ in H in (6), ~ is given by (17b) using b, ~ and also / where R
The normalizing constant for the posterior PDF in (15) is called the pseudo-evidence (or pseudo-marginal likelihood) for the model classMðdÞ given by pseudo-data^ hu and it can be analytically evaluated:
pð^hu jdÞ ¼
Z
pð^hu jh; dÞpðhjdÞdh ¼
Z
T
pð^hu jh; aÞpðhjx ; /; bÞdh ð18Þ
1
The second factor in (14) is given by Bayes’ Theorem:
ð19Þ
We assume that the model class is globally identifiable based on ^ h ^u jdÞ in ^ 2 ; w; the modal data, meaning here that the likelihood pðx (19) has a unique global maximum over d and then so does the posterior at ~ d (the MAP value of d). We then treat d in (14) as a ‘nuisance’ parameter vector and integrate it out using Laplace’s asymptotic approximation [40]:
^ ^hu Þpðdjx ^ ^hu Þdd pðhj~d; x ^ ^hu Þ ^ 2 ; w; ^ 2 ; w; ^ 2 ; w; pðhjd; x ð20Þ
^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; where ~ d ¼ arg max pðdjx hu Þ ¼ arg max½pðx hu jdÞ pðdÞ because the denominator in (19) is independent of d. In the next ^ ^ ^ 2 ; w; section, we will focus on the task of maximizing pðdjx hu Þ for ~ finding the MAP value d. ^ ^ 3.2. MAP estimate of d from PDF pðdjx ; w; hu Þ ^2
Finding the MAP value ~ d is equivalent to maximizing log J1(d) with respect to d where:
^ ^hu jdÞpðdÞ ¼ pðx ^ gÞpð^hu jx2 ; /; b; aÞ ^ 2 ; w; ^ 2 jx2 ; qÞpðwj/; J 1 ðdÞ ¼ pðx pðx ; /jbÞpðqjsÞpðgjmÞpðbjb0 Þpðb0 jjÞ 2
3
0 .. .
7 7 5
~Þ x Kðl
~ 21 M
ð24Þ
N m N d Nm Nd
~ a ~ x ~ 2 in b and H in (5) ~ is given by (17a) using b, ~ and also /, where l and (6), and for q = 1, . . ., NmNd:
Tq ¼ PTq Pq
Uq ¼ PTq
NX m Nd
! ~ s Ps /
s–q
^2 ^ ^ ^ ^hu Þ ¼ pðx ; w; hu jdÞpðdÞ pðdjx ; w; ^ h^u Þ ^ 2 ; w; pðx ^2
Z
~Þ x ~ 21 M Kðl .. .. . . 0
where D = A + b (H H) and we have utilized the hierarchical structure exhibited in Fig. 1. This is a Gaussian distribution over the Nh-dimensional pre-damage parameter vector hu evaluated at ^u , and is readily evaluated for arbitrary values of d. its MAP value h
^ ^hu Þ ¼ ^ 2 ; w; pðhjx
6 ~¼6 F 4
ð23Þ
2
1 ¼ Nð^hu jðHT HÞ HT b; DÞ 1
2
1
ð21Þ
where we have used the hierarchical structure exhibited in Fig. 1. Substituting equations (10), (18), (4), (12), (11), (7) and (8) into (21), it is shown in the Appendix that explicit expressions can be obtained for iterative maximization of log J1(d) where we update all of the parameters in d successively and repeat until convergence. By maximizing log J1(d) with respect to the system mode shapes / and the associated hyper-parameters g and m successively, sub-
Pq ¼
@H @/q
ð25Þ
~ requires repeated evaluation of the Determining the MAP value/ ~T F ~þg ~F ~ ~ 1 Þ; . . . ; trðRT ~ Nm N ÞÞÞ for dif~ CT C þ bdiagðtrð inverse of X ¼ ðb RT d ferent trial values for the other parameters. The dimension NmNd of / will often be large in applications (e.g., NmNd > 100). To facilitate the goal of efficient computation of the inverse of the matrix X, we rewrite it in the block diagonal form: X ¼ diagðX1 . . . ; XNm Þ; where the Nd Nd matrices Xi for i = 1, . . ., Nm are given by: ~ ~ ~ ði1ÞN þ1 Þ;. . . ;trðRT ~ iN ÞÞ þ g ~ 2 MÞ2 þ bdiagðtrð ~ N m CT Cs ; Xi ¼ bðKð l~ Þ x RT i
~ is given by: and then the MAP value /
d
d
s
~ ¼ diag X1 ; . . . X1 g ^ b ~ trðRU ~ 1 Þ; . . . ; trðRU ~ Nm N Þ T ~ CT w / 1 Nm d ð26Þ where Cs is the sub-matrix which contains the elements of the first ~ now involves Ns rows and Nd columns of C. The estimation of / inversions of Nm square matrices Xi (i = 1, . . ., Nm) of order Nd, ~ in (22) involves an inversion whereas the original estimation of / of one square matrix of order (NmNd). The latter requires a factor of N2m more computational effort than (26). Similarly, the MAP estimates of the system natural frequencies x2 and their associated hyper-parameters q and s are also derived as (using (A3), (A5) in (A1) in the Appendix):
h
i1 ~ T ~c ~ 1 x ~G ^2 þb LT E
~T G ~ ~ 1 L þ b ~G ~ 2 ¼ LT E x
and for each i = 1, . . ., Nm:
ð27Þ
42
Y. Huang et al. / Structural Safety 64 (2017) 37–53
Ns 2
q~ i ¼ P Ns r¼1
1 2 ¼ ~ si 2
ð28Þ
^ 2r;i x ~i x
6 ~¼6 G 4
~1 M/ .. . 0
3
0 .. .. 7 7 . . 5 ~N M/ m N d N m Nm
ð29Þ
ð30Þ
N d N m 1
Similarly, the MAP estimates of the hyper~ j ; j ¼ 1; . . . ; N h ; are derived as (using (A3) and (A11) parametersa in (A1) in the Appendix): 2 ~ þ ð^hu l ~ Þj a~ j ¼ ðRÞ jj
ð31Þ
~0 and j ~ b ~ are obtained by taking the Finally, the MAP values b, derivatives of log J1(d) with respect to b, b0 and j and substituting the MAP estimates of j into that of b and b0 (see (A13), (A15) and (A16) in the Appendix):
~¼ b
Nd Nm
PNh j¼1
~ ~ 1 1a j Rjj
~ þ 2b ~1 ~l ~ bk kH 0
~0 ¼ b=2 ~ ¼ 1=j ~ b
¼
log½pð^hu jx2 ; /; b; aÞpðhj^hu ; x2 ; /; b; aÞdh log½pðhjx2 ; /; bÞpð^hu jh; aÞ= pðhj^hu ; x2 ; /; b; aÞpðhj^hu ; x2 ; /; b; aÞdh log½pð^hu jh; aÞpðhj^hu ; x2 ; /; b; aÞdh Z log½pðhj^hu ; x2 ; /; b; aÞ=pðhjx2 ; /; bÞ pðhj^hu ; x2 ; /; b; aÞdh
~1 ~ Þ/ Kðl 7 6 . 7 ~c ¼ 6 .. 5 4 ~N ~ Þ/ Kðl m
Z
Z
3
2
Z
¼
where:
2
log½pð^hu jdÞ ¼
2
ð32Þ
ð33Þ
~ is inversely dependent on the A nice feature of (32) is that b square of the scaling of the mode shape and the units of the stiffness and mass matrices, which leads to the desirable property that the inferred stiffness reduction results are scale-invariant, since the terms bHTH and bHTb in (17a) and (17b) are independent of these scale selections. The MAP estimate ~ d can be calculated by going through the sequence of equations (26), (23), (27), (28), (31), (32) and (33) and then updating the Gaussian posterior in (20) for h by substituting ~ d into (17a,b) before repeating this process. For damage detection with real data where there are large modeling errors, it was sometimes observed during iterative solution for the MAP estimates ~ d that the components of the estimated sys~ kept decreasing towards zero while the tem mode shape / ~ tended to infinity. To avoid equation-error precision parameterb ~0 for the iterative procedure this occurring, a good initial choice of b
~0 is is required, which is discussed at the end of the Appendix, and b updated separately in an outer loop until convergence of the optimization occurs for the other parameters in d. 3.3. Information-theoretic interpretation of producing sparseness in the inferred stiffness reductions
For maximization of the pseudo-evidencepð^ hu jdÞ with respect to d, there is an information-theoretical interpretation of the trade-off between data fitting and model complexity [21,36] that can be demonstrated as follows, where we use Bayes’ Theorem and the hierarchical model structure in Fig. 1:
ð34Þ
Eq. (34) shows that the log pseudo-evidence, which is to be maximized, is the difference between the posterior mean of the log likelihood function for h (the first term) and the relative entropy (or Kullback-Leibler information) of the posterior with respect to the prior distribution for h (the second term). The first term quantifies the ability of the hierarchical model to match the MAP vector ^ hu which is obtained from the calibration state; it is
hu;j in the posterior maximized if all hj are tightly clustered around ^ for h (i.e. all aj ? 0). The second term reflects the amount of information extracted from the pseudo-data ^ hu ; it is minimized if h is
independent of ^ hu ; conditional on (x2, /, b) (i.e. all aj ? 1). This term penalizes models that have fewer parameter components hj ^u ; and therefore forces the model updating differing from those in h
^u and relatively more from the to extract less information fromh system modal parameters x2 and /, and so from the ‘‘measured” ^ as implied by the observations of the hierar^ 2 and w, modal data x chical model structure exhibited in Fig. 1. However, overextraction of information from the system modal parameters x2 and / will produce a structural model vector h with too large of a difference from ^ hu that is overly sensitive to the details of the
information in the specified system modal parameters x2 and /, ^ In ^ 2 and w: and therefore in the ‘‘measured” modal parameters x other words, the sensor noise and other environmental effects may not be ‘‘smoothed out” so they may have a large effect on the damage detection performance. To summarise, (34) demonstrates that the evidence maximization procedure over hyper-parameter vector d automatically involves a trade-off between the average data-fit of the model class defined by d and the information it extracts from the data. This is the Principle of Model Parsimony or the Bayesian Ockham razor [21,37] at work. The Bayesian procedure is effectively implementing Ockham’s Razor by assigning lower probabilities to a structural model whose parameter vector h has too large or too small differ^u obtained at the calibration state (too few or too ences from h
many aj ? 0), therefore suppressing the occurrence of false and missed damage detection alarms.
4. Proposed algorithms for identification of stiffness loss The above formulation leads to the following algorithms for probabilistic inference of structural stiffness loss. There are two variants of the algorithm. For the calibration stage, model sparseness is not expected and hence Algorithm 1 is used without optimization of the hyper-parameter vector a; instead, we fix all components aj with some large values. For the monitoring stage, Algorithm 2 is used where all hyper-parameters in d are estimated to ensure sparse stiffness loss inference.
43
Y. Huang et al. / Structural Safety 64 (2017) 37–53
4.1. Algorithm 1: Bayesian inference for the structural model in the calibration stage
structural stiffness reductions in the monitoring stage, which is summarized in Algorithm 2.
Algorithm 1
Algorithm 2
INPUTS:
INPUTS:
^ ¼w ^ u for Nm ^2 ¼ x ^ 2u and w . Identified Ns sets of modal data x modes from the calibration stage (N s P 3) . A chosen nominal value h0 of hu (e.g. from a finite-element model) and associated prior covariance matrix R0 ¼ r20 INh with large variance
r20 .
^ 2; ~ ¼ 10ðN s N m N o 2Þ=kwk . Initial parameter values for g PN s 4 ^ r;i ; i ¼ 1; . . . ; N m , q~ i ¼ 10ðNs 2Þ= r¼1 x PN s ~ ~ , where w ~ ¼ 20b ^ r;1 k2 and b ^ 21 MÞ b0 ¼ 2N s N o = r kðKðh0 Þ x 0 K and M are the reduced stiffness and mass matrices that correspond to the measured DOFs ALGORITHM: P s 2 ^ r =N s ~ 2 ¼ Nr¼1 1. Initialize l = h0, R = R0 and x x ~ ~ from input values ~ ; b0 and b ~; q 2. Initialize g ~ with large values (e.g. a ~ j ¼ 109 ) 3. Fix all components in a 4. While convergence criterion on ^ hu is not met (Outer loop) ~ is not met (Inner 5. While convergence criterion on b loop) ~ using (26), and then update g ~ using (23) 6. Update MAP / ~ 2 using (27), and then update q ~ using 7. Update MAP x (28) ~ using (32) 8. Update MAP b 9. Calculate the conditional posterior mean ^ hu ¼ l and covariance matrix Ru = R for h = hu using (17a) and (17b) ~ between the current 10. End while (the ratio of change of b iteration and the previous iteration of the inner loop is sufficiently small (e.g. smaller than 0.001)) ~0 using (33) 11. Updateb 12. End while (the ratio of change in all components of ^ hu between the current iteration and the previous iteration of the outer loop is sufficiently small (e.g. smaller than 0. 01)) OUTPUTS: . MAP estimates of all uncertain parameters in d . Posterior mean ^ hu and covariance matrix Ru of the stiffness scaling parameters hu at the calibration stage
^ ¼w ^ d for Nm ^2 ¼ x ^ 2d and w . Identified Ns sets of modal data x modes from the monitoring stage (N s P 3) . MAP estimate ^ hu and the associated posterior covariance matrix Ru from the calibration stage . Initial parameters values for ~ ¼ 2N s N o =PNs kðKð ~ , and ^ w ~ ¼ 20b ^ r;1 k2 and b ^ 21 MÞ b hu Þ x 0 0 r ~ i are the same as in Algorithm 1 ~ and q initial values for g ~ u of the system natural frequencies ~ 2 and / . MAP estimates x u
and mode shapes from the calibration stage ALGORITHM: ~ j ¼ N 2h , j = 1, . . ., Nh 1. Initialize a ~ as their MAP estimates from the cal~ 2 and / 2. Initialize x ibration stage ~0 and b ~ from input values ~; b ~; q 3. Initialize g 4. While convergence criterion on ~ hd is not met (Outer loop) ~ is not met (inner 5. While convergence criterion on a loop) ~ j ¼ 0 if aj (j = 1, . . ., Nh) becomes smaller than amin 6. Set a (chosen as 109 in the example later) ~ j ¼ 0 if ~ 7. Set a hj > ^ hu;j ðj ¼ 1; . . . ; N h Þ when the iteration number l > 2 8. Calculate the conditional posterior MAP value ~ hd ¼ l and covariance matrix Rd = R using (17a) and (17b) 9. Do steps 6–8 in Algorithm 1 ~ using (31) 10. Update a ~ 11. End while (the ratio of change in all components of a between the current iteration and the previous iteration of the inner loop is sufficiently small (e.g. smaller than 0.001)) ~0 using (33) 12. Update b 13. End while (the ratio of change in all components of ~ hd between the current iteration and the previous iteration of the outer loop is sufficiently small (e.g. smaller than 0.01)) OUTPUTS: . MAP estimates of all uncertain parameters in d . Posterior mean ~ hd and covariance matrix Rd of the stiffness scaling parameters hd at the monitoring stage
4.2. Algorithm 2: Sparse Bayesian learning for inference of structural stiffness reductions in the monitoring stage When real data from a structure is used, the modeling uncertainty in the structural stiffness may be too large to have the stiffness reduction reliably detected by most existing methods. To regularize a possibly severe ill-posed inverse problem, Ching and Beck [12,14] imposed a constraint in their Expectation–Maximiza tion algorithm that the values of the stiffness parameters for a damaged state cannot exceed those for the calibration state. However, we have not found an analytically tractable solution for SBL when imposing this additional prior information in the formulation. Therefore, we present an alternative strategy: a0j s for those substructures that have their MAP estimate ~ hj > ^ hu;j for the current iteration will be set to zero directly (the first two iterations will be discarded to make sure that the parameters in d are not directly influenced by the initial setting). Then the prior information that any stiffness increase from the calibration model is unacceptable is effectively utilized in our SBL procedure for inference of
Remark 4.1. Algorithm 2 is performed by iterating between the two groups of parameters, hd and d. It starts by considering all substructures as possibly damaged (the reason for the choice a~ j ¼ N2h ; j = 1, . . ., Nh, at the first iteration) and then causes the ~ j < 109 ; to be exactly ‘‘inactive” components hj of hd, which have a ^u;j from the calibration stage when optimizing over the equal to h hyper-parameters aj. After convergence, there are only a few ‘‘active”h0j s that are changed from their calibration values ^ hu;j and their corresponding substructures are considered to be damaged.
Remark 4.2. After obtaining the posterior PDFs of the stiffness scaling parameters hd and hu for both the current monitoring stage (where damage is possible) and the undamaged calibration stage, respectively, the posterior mean and standard deviations for hd and hu can be used to construct an approximation for the integrals involved (see (54) of [39]), allowing computation of the probability
44
Y. Huang et al. / Structural Safety 64 (2017) 37–53
that each substructure stiffness parameter has been reduced by more than a specified fraction of its value in the initial calibration stage [2,3,12,15,39].
explicitly addressed in the new SBL algorithm. We discuss these issues in this section. 5.1. Theory
Remark 4.3. The chosen Gaussian likelihood function for the stiffness scaling parameter h in (9) allows the stiffness to increase and ^u with equal chance. For many decrease from the calibration state h structures of interest, a prior judgment can be made that the structure is expected to lose stiffness where it is damaged. For example, structural damage in steel members tends to be at connections with local buckling or weld fracture, while in R/C members, damage usually occurs because of spalling of concrete and breaking of the bond between the reinforcing steel and the surrounding concrete. In these cases, laboratory tests of structural members have shown that there is a reduction in stiffness. For such cases, we impose the constraint in Step 7 of Algorithm 2 to alleviate the ill-conditioning in real damage inversions, that is, the solution for the MAP parameter~ d in Algorithm 2 is obtained by solving the following constrained optimization:
^ ^hu Þ; subject to ^ 2 ; w; max pðdjx
lj 6 ^hu;j ðj ¼ 1; . . . ; Nh Þ
where the posterior mean l of hd is a function of d. Even though ~ j ¼ 0 for any ~ it is a ‘‘greedy” optimization if we set a hj > ^ hu;j ðj ¼ 1; . . . ; N h Þ during any iteration (except the first two, once an aj = 0 in a certain iteration, it stays at zero), the strategy suppresses those local maxima that correspond to models having some substructures with stiffness increases from the calibration state; this significantly increases the probability of finding the global maximum with respect to d. In addition, under the condition that the modeling errors are zero-mean, there is a possible compensation between the inferred stiffness increases and decreases for different substructures that result from modeling errors alone. Therefore, by suppressing stiffness increases, the strategy helps to reduce the occurrence of inaccurate stiffness reductions and also induces a higher degree of sparseness in the identified stiffness reduction results. Remark 4.4. In our new formulation, the Gaussian and Exponential PDFs have a scale invariant property and their associated hyper-parameters are also properly learned. This leads to the desirable property that the algorithms are automatically scaleinvariant: the identification of the stiffness scaling parameters is independent of any scaling of the measured mode shapes and of the units chosen for the mass and stiffness matrices. Our algorithm in [39] does not have this scale-invariance property and so care must be taken in the selection of the scale parameter b0 for the equation-error precision b.
5. Comparison of new SBL algorithm with that in [39] In [39], the hyper-parameters a for the likelihood function in (9) are learned by maximizing the evidence function R ^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; pðx hu jaÞ ¼ pðx hu jh; n; aÞpðh; njaÞdhdn along with the stiffness scaling parameter h and other parameters n = {(x2)T, qT, sT, /T, g, b} to produce sparseness in the inferred stiffness reductions. For tractable estimation of the evidence function, two ^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; approximations pðx hu jaÞ pðx hu j~ n; aÞ and ~ n^ n are uti2 ~ ¼ argmax pðnjx ^ h ^ is obtained by maxi^u Þ and n ^ ; w; lized, where n ^ h ^u Þ. The excellent performance of the ^ 2 ; w; mizing the PDF pðn; hjx proposed SBL method in [39] has been verified with synthetic data. However, there are several important issues that needed further exploration for reliable application to real data and they have been
The
approximation
of
the
evidence
function
^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; pðx hu jaÞ pðx hu j~ n; aÞ in [39] is valid only if ^ h ^u Þ has a unique maximum at ~ ^ 2 ; w; n with a sharp peak at its pðnjx mode. Furthermore, for real structural damage detection where there are large modeling errors, it is challenging to confirm the accuracy of the approximation ~ n^ n employed in the formulation in [39]. Therefore, these approximations to estimate the evidence function may possibly be unreliable. The new formulation presented in Sections 3.1 and 3.2 establishes the probability of localized stiffness reductions without these possibly unreliable approximations and is therefore more rigorous. Another benefit of the new formulation is that all parameters in d are learned by an evidence maximization procedure and therefore the inferred sub-structural stiffness reductions are expected to be sparser than in [39]. Note that due to the hierarchical nature ~ 2 andq ~; x ~ i are estimated of the model in Fig. 1, the MAP values of g in the same way as those in [39]. However, the new terms ~ RU ~ diagðtrðRT ~ Nm N ÞÞ; ~ 1 Þ; . . . ; trðRU ~ Nm N Þ ~ 1 Þ; . . . ; trðRT b½trð and b d d PNh 1 ~ ~ ~ ~ R in (22) and (32) for the MAP values / and b are 1 a jj j¼1 j important for the algorithm to produce more sparse substructural stiffness reductions that help to suppress the possible occurrences of false damage detections. In addition, the information-theoretic explanation for the sparseness in the inferred stiffness reductions in Section 3.3 is important to interpret how the proposed algorithm suppresses the occurrences of false and missed damage detections, which was not discussed in [39]. 5.2. Computation for real applications Although the formulation in [39] does not involve solving an eigenvalue problem during each iteration, the dimension (NmNd) of the system mode shape vector / is a concern for computational efficiency since its MAP estimation has a computational effort of OðN 3m N 3d Þ, which is computationally demanding when the dimension NmNd is not small. In the new algorithm, we enhance the computational efficiency by using (26) to reduce the computational cost by a factor of N 2m . The second computational improvement is that our new algorithm, unlike our algorithm in [39], has the desirable property of scale-invariance, i.e., the identified stiffness reductions are independent of the choice of units for both the stiffness and mass matrices, as well as the scaling of the measured mode shapes. This property is a consequence of choosing the value of b0 adaptively to scale the prior distribution of the equation-error precision b. In [39], we suggested an a priori selection of b0 based on an initial empirical study but, unfortunately, this can lead to unreliable results for both the undamaged and possibly damage states because they are somewhat sensitive to b0. In the new algorithm, b0 is estimated by its MAP value, which is given by (33), and it is updated separately in an outer loop initialized from a reasonable value as specified in (A18). The numerical and experimental results given later support the effectiveness of this choice of b0. For application to real structural response data, the modeling error sometimes dominates over the useful data information for constraining the parameter updating and then it becomes challenging to implement the damage detection algorithm in [39]. For example, an increase in the identified substructure stiffness can sometimes occur. In the new algorithm, a more effective regu-
Y. Huang et al. / Structural Safety 64 (2017) 37–53
larization strategy is proposed based on our prior expectations in order to further alleviate the ill-conditioned stiffness inversion by suppressing mean stiffness increases during the optimization to find the MAP estimates of the parameters. 6. Illustrative applications The applicability of the proposed methodology to identify substructure stiffness reductions is illustrated with the IASC-ASCE Phase II SHM benchmark problems, first using the simulated data [41], and second using the experimental data [42]. The first application serves to verify the proposed methodology while the second one illustrates its application to real data from a damaged structure. The benchmark structure is a four-story, two-bay by twobay steel braced-frame structure and a diagram for the analytical model is depicted in Fig. 2 along with its dimensions, in which the x-direction is the strong direction of the columns. A detailed description of the benchmark problem, including detailed nominal properties of the structural elements in the analytical model, can be found in [12]. 6.1. Simulated Phase II benchmark problem Simulated Phase II benchmark problem is defined in Bernal et al. [41] and consists of both brace damage and unbraced beam-column joint damage cases (see [12] for the various damage pattern cases considered here). For instrumentation scenarios, only the results of the challenging partial-sensor scenario are presented, where only the measurements at the third floor (mid-height) and the roof are available. For generation of simulated time-domain data and the extraction of the modal properties, the reader can refer to [12] or [14] for detailed information. To locate the faces sustaining brace damage, the identification model class is defined with the stiffness matrix K parameterized as
KðhÞ ¼ K0 þ
XX s
v
hsv Ksv
ð35Þ
where s = 1, . . ., 4 refers to the story number and v = +x, x, +y, y indicates the direction of the outward normal of each face at each floor, yielding a stiffness scaling parameter vector h with sixteen components, corresponding to sixteen substructures (four faces of four stories). The K0sv s are the ‘nominal’ stiffness matrices to make the nominal value of each stiffness scaling parameter to be 1.0.
Fig. 2. The diagram of the benchmark structure [12,14].
45
For the braced structure, they are computed based on shearbuilding assumptions by projecting the full stiffness matrix for the 120-DOF benchmark structural model [41] onto 12 DOFs for a 3-D shear-building model with rigid floors (three DOFs per floor: translations parallel to the x- and y-axes and rotation about the zaxis). K0 is the nominal stiffness matrix contribution from the columns and beams. For the unbraced structure, we project the 120-DOF benchmark structural model [41] onto 20-DOFs for a 3-D shear-building model that assumes rigid floors in the x–y plane and allows rotation along the x- and y-axes (five DOFs per floor: translations parallel to the xand y-axes and rotation about the x-, y- and z-axes). It is assumed that the rotational stiffnesses of all beam-column connections in the same story along the x-axis (or along the y-axis) are identical. Two parameters are used for the rotational stiffness in each story to give eight stiffness parameters in total:
KðhÞ ¼ K0 þ
XX hsu Ksu s
ð36Þ
u
where s = 1, . . ., 4 refers to the story number and u = x, y indicates the axis along which the rotational stiffness is active. The Ksu are the nominal rotational stiffness matrices computed based on the model assumptions for the original undamaged structure; K0 is the nominal stiffness matrix contribution from the columns and beams. Ideally, we would like to treat each structural member (e.g., brace, beam or column) as a substructure so that we can infer from the dynamic data which, if any, members have been damaged by the event. However, the information available from the structure’s local network of sensors will generally be insufficient to support a member-level resolution of stiffness loss from damage, so larger substructures consisting of assemblages of structural members may be necessary in order to reduce the number of model parameters in h. A tradeoff is therefore required between the number of substructures (and hence the resolution of the damage locations) and the reliability of the probabilistically-inferred damage state. In this simulated data example, we only investigate the partialsensor scenario. During the calibration stage, we utilize Algorithm 1 based on 100 independent sets of identified modal parameters from the undamaged structure to find the MAP value of the structural stiffness scaling parameters^ hu and the corresponding c.o.v. (coefficients of variation), which are shown in Tables 1 and 2, for the undamaged braced structure (column RB.ps) and unbraced structure (column RU.ps), respectively. It is evident from the tables that the MAP estimates of all stiffness scaling parameters hj are accurate, that is, the errors are smaller than 3.7% and 2.5% for the braced and unbraced structure, respectively. ^u from During the monitoring stage, we choose the MAP value h the calibration stage as pseudo-data for h and perform Algorithm 2 based on ten sets of identified modal parameters as the primary data. The identified stiffness ratios of the MAP stiffness parameters with respect to those for the undamaged case (RB) are tabulated in Tables 1 and 2, for the four brace and three joint damage patterns, respectively, along with their corresponding c.o.v. The actual damaged locations are made bold for comparison. It is found that the damage patterns are reliably detected in both qualitative and quantitative ways. The identified stiffness ratios for the actual damaged substructures are close to their actual values (0.943 and 0.887 for the brace damage cases and 0.500, 0.667, 0.750 and 0.833 for the joint damage cases, see [12]) and all components with non-bold font have stiffness ratios that are exactly equal to one. In addition, the corresponding c.o.v. for the exactly identified non-damaged components is equal to zero, which means these substructures have no stiffness reduction with full confidence compared with that of the calibration stage. This is a benefit of
46
Y. Huang et al. / Structural Safety 64 (2017) 37–53
Table 1 Brace damage cases in Simulated Phase II benchmark: Identification results from Algorithms 1 and 2 for the partial-sensor scenario. Parameter
h1,+x h2,+x h3,+x h4,+x h1,+y h2,+y h3,+y h4,+y h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
RB.ps
DP1B.ps
DP2B.ps
DP3B.ps
MAP value ð^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
DP3Bu.ps c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
0.984 1.010 1.026 0.963 0.981 0.992 1.005 1.012 0.997 0.985 0.976 1.021 0.998 0.995 0.988 1.006
0.192 0.094 0.094 0.075 0.235 0.122 0.100 0.098 0.325 0.162 0.171 0.125 0.221 0.117 0.093 0.097
1.000 1.000 1.000 1.000 0.865 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.896 1.000 1.000 1.000
0.000 0.000 0.000 0.000 0.203 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.196 0.000 0.000 0.000
1.000 1.000 1.000 1.000 0.941 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.936 1.000 1.000 1.000
0.000 0.000 0.000 0.000 0.183 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.179 0.000 0.000 0.000
1.000 1.000 1.000 1.000 0.858 1.000 0.952 1.000 1.000 1.000 1.000 1.000 0.886 1.000 0.954 1.000
0.000 0.000 0.000 0.000 0.208 0.000 0.087 0.000 0.000 0.000 0.000 0.000 0.201 0.000 0.081 0.000
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.867 1.000 0.957 1.000
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.183 0.000 0.069 0.000
⁄.ps denotes the partial-sensor scenario. ⁄The true ratio value for h1,+y for DP1B and DP3B and h1,y for DP1B, DP3B and DP3Bu are 88.7% of the undamaged values. ⁄The true ratio value for h1,+y for DP2B, h3,+y for DP3B, h1,y for DP2B and h3,y for DP3B and DP3Bu are 94.3% of the undamaged values.
Table 2 Joint damage cases in Simulated Phase II benchmark: Identification results from Algorithms 1 and 2 for the partial-sensor scenario. Parameter
RU.ps
h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
DP1U.ps
DP2U.ps
DP1Uu.ps
MAP value ð^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
0.984 1.010 1.026 0.963 0.981 0.992 1.005 1.012
0.192 0.094 0.094 0.075 0.235 0.122 0.100 0.098
1.000 1.000 1.000 1.000 0.528 0.615 1.000 1.000
0.000 0.000 0.000 0.000 0.902 1.346 0.000 0.000
1.000 1.000 1.000 1.000 0.737 1.000 1.000 1.000
0.000 0.000 0.000 0.000 0.633 0.000 0.000 0.000
1.000 1.000 1.000 1.000 0.797 0.799 1.000 1.000
0.000 0.000 0.000 0.000 0.595 1.085 0.000 0.000
⁄The true ratio values for h1,y and h2,y for DP1U.ps are 50% and 66.7% of the undamaged values, respectively. ⁄The true ratio values for h1,y for DP2U.ps are 66.7% of the undamaged values. ⁄The true ratio values for h1,y and h2,y for DP1U.ps are 75% and 83.3% of the undamaged values, respectively.
the proposed sparse Bayesian formulation, which reduces the uncertainty of the unchanged components. It produces sparse ~ sv ! 0 models by learning the hyper-parameter a, where a ~ su ! 0Þ implies that R(sv)(sv) ? 0 (or R(su)(su) ? 0) and ðor a ~ hu Þ ðor ~ hsu ! ð^ hu Þ Þ. hsv ! ð^ sv
su
In order to further portray the identification performance of the proposed algorithms, the identified stiffness reduction ratios of Algorithm 2 are compared with those of our previous Bayesian inversion algorithm in [39] for various damage patterns, and the results are presented in Fig. 3, where the stiffness reduction ratio is computed as the difference between each MAP value ^u Þ ðor ðh ^u Þ Þ of the stiffness scaling parameters from the caliðh sv
su
~su ) from the current monbration stage and the MAP value ~ hsv (or h ^u Þ (or ð^ h itoring stage, normalized by ðh u Þsu ). A variant of Algorithm sv 2 with the stiffness increase constraint removed (Step 7 is omitted) is also studied in order to investigate the effect of this constraint on the identification performance. The results for the modified Algorithm 2 without constraints and for Algorithm 2 are presented in Fig. 3(a),(b) in the columns (iii) and (iv), respectively, while for the method in [39], the results with the Laplace-prior parameter k = 0 and k optimized are given in Fig. 3(a),(b) in the columns (i) and (ii), respectively. Note that a uniform hyper-prior over a is utilized for Algorithm 2, which is equal to the Laplace prior case with k ¼ 0: In Fig. 3, the 16 (or 8) components hsv (or hsu) of h are presented in the order h1, . . ., h16 (or h1, . . ., h8) that corresponds to the order listed in the first column of Table 1 (or Table 2).
The results for Algorithm 2 without constraints in Fig. 3(b)(iii) are more sparse in comparison with the results of the method in [39] with k ¼ 0 (Fig. 3(b)(i)), which shows the advantage of the new formulation in Section 3 as discussed in Sections 5.1 and 5.2. The results for Algorithm 2 are even better with no false damage detections observed, and hence an appropriate level of sparseness of stiffness reductions is inferred by imposing the constraint in Step 7 to suppress stiffness increases. The algorithm in [39] with the optimal Laplace prior parameterk has the same perfect damage assessment performance as Algorithm 2 for the brace damage patterns (compare Fig. 3(a)(ii),(iv)) but, as shown in Fig. 3(b)(ii),(iv), this is no longer true for the joint damage patterns; in fact, Fig. 3 (b)(ii) shows that totally over-sparse stiffness reduction models with no nonzero components are produced for the algorithm with optimalk in the DP2U.ps and DP1Uu.ps cases. This is because when finding the MAP value of a by maximiz^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; ingpðajx hu ; kÞ / pðx hu jkÞpðajkÞ, the Laplace hyperpriorpðajkÞ in the method of [39] will dominate over the evi^ h ^u jaÞ in these cases. ^ 2 ; w; dencepðx A more complete picture of substructure damage is given by the dam probabilities of damage, Pdam sv ðf Þ (or P su ðf Þ), with different fractional stiffness losses f e [0, 1], calculated using the theory presented in [39]. In Fig. 4, the probability of damage is shown for the sixteen h0sv s (or eight h0su s) for various damage patterns from applying Algorithm 2. It is clear in all cases that the substructures with real stiffness losses have damage with a probability of unity.
47
Y. Huang et al. / Structural Safety 64 (2017) 37–53
DP1B.ps
0.2
0.1
0.1
0.1
0
0
0
0.2 0.1
DP2B.ps 0.2 DP2B.ps 0.1
0 0.2
0.2
DP3B.ps
0 4 8 12 16
DP2B.ps
0.2
DP3B.ps
DP3B.ps
0.1
0.1
0.1
0.1
0
0
0
0
0 4 8 12 16
0 4 8 12 16
0 4 8 12 16
DP3B.ps
0 4 8 12 16
0.2
0.2 DP3Bu.ps 0.2 0.2 DP3Bu.ps DP3Bu.ps DP3Bu.ps 0.1 0.1 0.1 0.1 0
0 0 4 8 12 16 (i)
0
0 0 4 8 12 16 (ii)
0 4 8 12 16 (iii)
0 4 8 12 16 (iv)
(a)
DP 1U.ps
0.8
DP1U.ps
DP 1U.ps
0
0 0 2 4 6 8
0.8 0.4
0.4
0 0 2 4 6 8
0.6 0.4 0.2 0
0.8 0.4
0
0 4 8 12 16 0.2
DP1U.ps
0.4
0 0 4 8 12 16
0.2
0.8
DP2B.ps
0.1
0 0 4 8 12 16
0.2
0.2 DP1B.ps 0.1 0
0.1
0 0 4 8 12 16
DP1B.ps
0 4 8 12 16
0 4 8 12 16
0 4 8 12 16 Stiffness reduction ratios
0.2
DP1B.ps
Stiffness reduction ratios
0.2
0 2 4 6 8
0 2 4 6 8
0.6 0.6 0.6 DP2U.ps DP2U.ps 0.4 DP2U.ps 0.4 DP2U.ps 0.4 0.2 0.2 0.2 0 0 0 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8
0.6 0.6 0.6 0.6 DP1Uu.ps DP1Uu.ps DP1Uu.ps DP1Uu.ps 0.4 0.4 0.4 0.4 0.2 0.2 0.2 0.2 0 0 0 0 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8
(i)
(iii)
(ii)
(iv)
(b)
Fig. 3. Simulated Phase II benchmark for (a) Brace damage cases and (b) joint damage cases: Comparison of stiffness reduction ratios from the two Bayesian inversion methods for the substructures for the investigated damage patterns using the data of the partial-sensor scenario. Column (i) Algorithm in [39] with k ¼ 0; (ii) Algorithm in [39] with k optimized in each iteration; (iii) Algorithm 2 with no constraints; (iv) Algorithm 2. (The bars colored red are ratios that agree with the actual damage locations. Also, the negative bars in (a), (b) (i) and (iii) correspond to small stiffness gains). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
For the substructures that are undamaged, most of the median fractional stiffness losses are close to zero, and the probability is very small that the damage extent exceeds 2%, which is within the tolerable uncertainty level for applications.
6.2. Experimental Phase II benchmark problem In the second example, the proposed methodology is applied to both the brace and joint damage patterns in the IASC-ASCE experimental Phase II benchmark problem [42]. Nine configurations are considered, in which Configs. 1–6 are five brace damage patterns with Config. 1 being the calibration (undamaged) case, and Configs. 7–9 are two joint damage patterns for the unbraced structure with Config. 7 being the calibration case. The five brace damage patterns considered are as follows: 1) Config. 2: removal of all braces (100% stiffness loss) on the –y face; 2) Config. 3: removal of the left-handside brace in each story on the –y face; 3) Config. 4: removal of the left-hand-side braces in the first and fourth stories on the –y face; 4) Config. 5: removal of the left hand-side brace in the first story on the –y face; 5) Config. 6: removal of two braces in the second story on the + x face. For the two joint damage patterns, we have: 1) Config. 8: loosen both ends of the right-hand-side beam at each floor on the y face in Fig. 2; 2) Config. 9: loosen both ends of the right hand-side beam at the first and second floors on the y face. For each configuration, experimental acceleration data from the benchmark structure measured after impact of a sledgehammer are utilized. Sensor data were collected from five Kinemetrics EPI sensors placed near the base and floor centers (Nodes 5, 14, 23, 32, and 41 in Fig. 2) that sensed the accelerations in the +y direction, and from ten force–balance accelerometers mounted on the +y and y faces of the base and all floors (Nodes 2, 8, 11, 17, 20, 26, 29, 35, 38, and 44 in Fig. 2) that sensed accelerations in the +x direction. For each configuration, the time-domain data from each sensor is divided into three segments and three independent sets of
experimental modal parameters are extracted by the modal identification procedure MODE-ID [43,44]. For each time segment, five modes (Nm = 5), consisting of the first and second translation modes in the x and y directions and the first torsion mode are identified for the brace damage cases, and eight modes (Nm = 8), consisting of the first, second and third translation modes in the x and y directions and the first and second torsion modes are identified for the joint damage cases. These modal identification results can be found in [12]. The same 3D 12-DOF and 20-DOF structural models adopted in the simulated data study are used for the brace and joint damage cases, respectively. We first run Algorithm 1 to get the MAP estimates of the stiffness scaling parameter ^ hu in the calibration stage. The available information of three independent sets of identified ^ u is insufficient to support an accurate ^ 2u and w modal parameter x inference of the stiffness scaling parameters for both cases. We therefore assume that the stiffness scaling parameters of +x and x faces for all floors, and +y and y faces for all floors, are identical in order to reduce the number of stiffness model parameters to two for the braced cases, while for the unbraced cases, we assume all stiffness scaling parameters are identical. During the monitor^u from the calibration stage is used as ing stage, the MAP value h pseudo-data for h and Algorithm 2 is implemented based on three ^ d . As in Section 6.1, ^ 2 and w sets of identified modal parameters x d
for brace and joint damage cases, all 16 and 8 components hj of h, respectively, are updated according to the order listed in the first columns of Tables 1 and 2. Based on the results from Algorithm 2, the stiffness ratios and the associated c.o.v. are tabulated in Tables 3 and 4 for the brace and joint damage cases, respectively. For the brace damage shown in Table 3, it is seen that all damage is effectively detected although there is a false detection for h16 = h4,y of Config. 5 with a 3.6% stiffness reduction; however, in this case, the difference in value for no damage (1.000) is less than one standard deviation (5.2%) from the MAP value of 0.964 for h4,y and so the probability of no damage is
48
Y. Huang et al. / Structural Safety 64 (2017) 37–53
1
0.8
1,+y
0.6
1,-y
Probability of damage
Probability of damage
1
0.4 0.2 0 -0.1
-0.05
0
0.05 0.1 Damage extent, f
0.15
0.2
0.8
1,+y
0.6
1,-y
0.4 0.2 0 -0.1
0.25
-0.05
0
0.05 0.1 Damage extent, f
(a)
0.6
1,-y
0.4
3,+y
0.2
3,-y
-0.05
0
0.05 0.1 Damage extent, f
0.15
0.2
Probability of damage
Probability of damage
1,+y
0.8 1,-y
0.6 0.4 0.2
3,-y
0 -0.1
0.25
-0.05
0
0.05 0.1 Damage extent, f
(c)
0.15
0.2
0.25
(d)
1
1 1,y
0.8
Probability of damage
Probability of damage
0.25
1
0.8
0.6 0.4 0.2 0 -0.4
0.2
(b)
1
0 -0.1
0.15
-0.2
0
0.2 0.4 0.6 Dama ge extent, f
0.8
1
0.8
1,y
0.6
2,y
0.4 0.2 0 -0.4
-0.2
0
0.2 0.4 0.6 Dama ge extent, f
(e)
0.8
1
(f) Probability of damage
1 0.8
1,y
0.6
2,y
0.4 0.2 0 -0.4
-0.2
0
0.2 0.4 Damage extent, f
0.6
0.8
1
(g) Fig. 4. Simulated Phase II benchmark: Probability of substructure damage exceeding f calculated using Algorithm 2 for the partial-sensor scenario: (a) DP1B.ps; (b) DP2B.ps; (c) DP3B.ps; (d) DP3Bu.ps;(e) DP1U.ps; (f) DP2U.ps and (g) DP1Uu.ps.
quite high. The MAP estimate of stiffness ratios for a particular face with removal of one and two braces are underestimated as approximately 50% and 30%, respectively, compared to the expected stiffness ratios of 77.4% and 54.9%, respectively [12]. For the quantification of the posterior uncertainty, the c.o.v. is generally much higher compared with those values in Table 1 for the simulated data benchmark due to larger modeling errors in this real
case, especially for those components corresponding to real damage locations. For the joint damage cases shown in Table 4, although most damage except for h3,y is detected, four false detections are found, and some of them indicate significant stiffness reduction (e.g., h4,x and h4,y of Config. 9). It is believed that the relative poor performance compared with the brace damage cases is due to the fact
49
Y. Huang et al. / Structural Safety 64 (2017) 37–53 Table 3 Brace damage cases in Experimental Phase II Benchmark: Identification results from Algorithms 1 and 2. Parameter
h1,+x h2,+x h3,+x h4,+x h1,+y h2,+y h3,+y h4,+y h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
Config. 1
Config. 2
Config. 3
Config. 4
Config. 5
Config. 6
MAP value ^u Þ ðh
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.302 0.370 0.268 0.160
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 36.148 58.243 32.306 45.287
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.443 0.646 0.386 0.293
0.000 0.000 0.000 0.000 1.593 0.000 0.000 0.000 0.000 0.000 0.000 0.000 28.930 34.911 25.493 33.986
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.767 1.000 1.000 0.608
0.000 0.000 0.000 0.000 0.000 1.981 0.000 0.000 0.000 0.000 0.000 0.000 22.028 0.000 1.359 22.014
1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 1.000 1.000 0.677 1.000 1.000 0.964
0.000 0.000 0.000 0.000 0.000 1.345 0.000 0.000 0.000 0.000 0.000 0.000 13.541 0.000 0.000 5.245
1.000 0.835 1.000 1.000 0.994 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
0.000 15.719 0.000 0.000 3.458 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
⁄The true ratio values for h1,y, h2,y, h3,y and h4,y for Config. 2 and h1,y for Config. 6 are 54.9% of the undamaged values. ⁄The true ratio values for h1,y, h2,y, h3,y and h4,y for Config. 3, h1,y and h4,y for Config. 4, h1,y for Config. 5 are 77.4% of the undamaged values.
Table 4 Joint damage cases in Experimental Phase II Benchmark: Identification results from Algorithms 1 and 2. Parameter
h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
Config. 7
Config. 8
Config. 9
MAP value ð^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
8.281 8.281 8.281 8.281 8.281 8.281 8.281 8.281
0.925 1.000 1.000 1.000 0.715 0.766 1.000 0.747
0.832 0.000 0.000 2.100 0.599 0.931 0.000 2.163
0.957 1.000 1.000 0.724 0.736 0.870 1.000 0.862
0.934 0.000 0.000 1.916 0.676 1.157 0.000 2.141
⁄The true ratio values for h1,y, h2,y, h3,y and h4,y for Config. 8 are 83.3% of the undamaged values. ⁄The true ratio values for h1,y and h2,y for Config. 9 are 83.3% of the undamaged values.
that the stiffness matrix of the structure is dominated by the columns and the installed sensors only measure the accelerations in the x and y directions, rendering the identification of rotational stiffness difficult. However, it is seen that the damage assessment performance here is better than in [14], since less false damage detections are produced and the identified stiffness ratios for the damaged substructures are closer to their actual values (0.833). In Fig. 5, the modified Algorithm 2 without constraints (Fig. 5 (c)) and Algorithm 2 (Fig. 5(d)) are compared with the method in [39] by presenting the stiffness reduction ratios for various damage cases, where hyper-parameterk ¼ 0 in Fig. 5(a) andk is optimized for the assigned Laplace prior in Fig. 5(b). It is interesting to see totally over-sparse stiffness reduction models with no nonzero components are produced when k is optimized. In comparison with the method in [39] with k ¼ 0; Algorithm 2 with no constraints has better performance (Fig. 5(c)), although false detections are still observed in h12 = h4,x for Config. 4, h16 = h4,y for Config. 6, and h4 = h4,x, h8 = h4,y for Config. 9 that correspond to significant identified stiffness reductions. In addition, some undamaged substructures show a significant increase in stiffness that is unrealistic. However, incorporating the additional constraint to suppress stiffness increases in Algorithm 2 results in exact localization of all stiffness reductions (Fig. 5(d)) for the brace damage patterns, although the modeling error seems to be too large to have the small effect of rotational stiffness damage reliably detected for
the joint damage patterns (especially Config. 9), which is an inherent difficulty of this particular damage identification case. Another reason for the better performance of Algorithm 2 compared with the method in [39] for the experimental data is that it automatically chooses the optimal scale parameter b0 whereas in [39] we select b0 by trial-and-error but the damage detection results can be sensitive to the value of b0. The damage probability curves for the stiffness scaling parameters h are plotted in Fig. 6 for Configs. 2–6 for brace damage patterns and Configs. 8–9 for joint damage patterns. The non-zero probabilities for negative fractional stiffness losses f are a result of the Gaussian approximations used in the calculation of the curves. For brace damage, all the actual damaged substructures are clearly shown to have a large damage probability when their fractional stiffness loss is close to the actual fraction. For the undamaged substructures, the posterior medians of the estimated damage corresponding to a damage probability of 0.5 are around zero, although in some cases small probabilities of about 0.2 are estimated for damage fractions exceeding 10%, which is due to the larger posterior uncertainty for the stiffness scaling parameters in the calibration and monitoring stages (Table 3). For the joint damage cases, h1,y, h2,y and h4,y for Config. 8 and h1,y, h2,y, h4,y and h4,x for Config. 9 are clearly shown to have a large damage probability, where false detection occurs in h4,y and h4,x for Config. 9 and missed detections occurs in h3,y for Config. 8 (see also Table 4).
50
Y. Huang et al. / Structural Safety 64 (2017) 37–53
1 0.5 0 -0.5 -1
Stiffness reduction ratios
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
Config2
0 4 8 12 16 Config3
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1
0 4 8 12 16 Config4
1 0.5 0 -0.5 -1
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
0
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
Config3
0 4 8 12 16 Config4
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
Config9
Config9
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
0 4 8 12 16 Config4
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
0 4 8 12 16 Config3
0 4 8 12 16 Config4
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
-0.5 0 2 4 6 8
0.5
Config2
0
Config9
0
-0.5
-0.5
1 0.5 0 -0.5 -1
Config3
-0.5
0
0
1 0.5 0 -0.5 -1
0 4 8 12 16 1 0.5 0 -0.5 -1
0 2 4 6 8 0.5
Config2
0
-0.5 0 2 4 6 8
1 0.5 0 -0.5 -1
0 4 8 12 16
0
-0.5 0.5
Config2
0 2 4 6 8 0.5
Config9
0 -0.5
-0.5
0 2 4 6 8
0 2 4 6 8
0 2 4 6 8
0 2 4 6 8
(a)
(b)
(c)
(d)
Fig. 5. Experimental Phase II benchmark: Comparison of stiffness reduction ratios from the two Bayesian inversion methods for the all substructures for five brace damage and two joint damage scenarios: (a) The algorithm in [39] with k ¼ 0; (b) The algorithm in [39] with k optimized in each iteration; (c) Algorithm 2 with no constraints; (d) Algorithm 2. (Note that a negative value corresponds to a stiffness gain).
7. Concluding remarks In the paper, we examine the problem of inferring damageinduced stiffness reductions in a structure from structural vibration response by taking a hierarchical Bayesian model perspective. For each chosen substructure, not only the most probable estimate of the spatially-sparse stiffness changes is inferred from the identified modal parameters but also the associated posterior uncertainties of the stiffness parameters are quantified, including the probability of substructure damage of various amounts. Our approach leads to an improved sparse Bayesian learning (SBL) algorithm for probabilistic damage detection compared to our previously proposed SBL algorithm in [39]. Apart from the three appealing features in the original method listed in [39], that there is no need to solve the nonlinear eigenvalue problem nor to match model and experimental modes explicitly, as well as the automated selection of the algorithmic parameters, the new method has the following four additional advantages. First, two possibly unreliable approximations in the original method in [39] have been eliminated. Second, the computation is more efficient in finding the MAP values of the system mode shapes. Third, the identified stiffness reductions are independent of the units of the stiffness and mass matrices, as well as the scaling of the identified
mode shapes. Fourth, a regularization strategy during hyperparameter estimation has been developed to further alleviate the ill-conditioning of the damage inversion for the case of large modeling errors when applied to real data. This strategy is based on the expectation that damaged structures actually lose stiffness at the damaged locations, so we incorporate an additional constraint during optimization for the MAP (maximum a posteriori) values of the parameters other than the stiffness scaling parameters, that the latter cannot have mean values in a damaged state that are greater than their mean values in the calibration state, which reduces the occurrences of incorrect identifications of stiffness changes. The improved effectiveness and robustness of the proposed method compared with our previously-published Bayesian method [39] was demonstrated through the analysis of the IASC–ASCE Phase II simulated and experimental benchmark problem to identify brace damage and damage to beam-column connections in the steel-frame benchmark structure. The results show that our new damage assessment method is reliable, despite significant error in the structural model in the real experiments. For future studies, it would be useful to explore the uncertainties in the identified system modal parameters and stiffness parameters from the calibration stage, since these are effectively ignored when only their MAP values are utilized.
51
Y. Huang et al. / Structural Safety 64 (2017) 37–53
1 4,-y
0.6
3,-y
0.4
2,-y
0.2
1,-y
0 -0.3
0
0.3 0.6 Damage extent, f
0.9
1.2
1
0.8
4,-y
0.6
3,-y
0.4
2,-y
0.2
1,-y
0 -0.3
0
0.3 0.6 Damage extent, f
(a)
1.2
0.6 0.4 1,-y
0.2 0 -0.3
0
0.3 0.6 Damage extent, f
1.2
(c)
1,-y Probability of damage
0.8 0.6 0.4 0.2 0
0.3 0.6 Damage extent, f
0.9
1.2
0.8 0.6 0.4 2,+x 0.2 0 -0.3
0
0.3 0.6 Damage extent, f
(d)
0.9
1.2
(e) 1
0.8
1,y
0.6
2,y
0.4
Probability of damage
1
4,y
0.2 0 -0.3
0.9
1
0 -0.3
Pro babili ty o f damage
4,-y
0.8
(b) 1
Pro babili ty o f damage
0.9
Probability of damage
0.8
P ro ba bili ty o f damage
Probability of damage
1
0
0.3 0.6 Damage extent, f
0.9
1.2
0.8
1,y
0.6
4,x
0.4
2,y 4,y
0.2 0 -0.3
0
0.3 0.6 Damage extent, f
(f)
0.9
1.2
(g)
Fig. 6. Experimental Phase II benchmark: Probability of substructure damage exceeding f using Algorithm 2 for the brace damage cases: (a) Config. 2; (b) Config. 3; (c) Config. 4; (d) Config. 5 and (e) Config. 6, and joint damage cases: (f) Config. 8 and (g) Config. 9.
Acknowledgements The first author is partially supported by the George W. Housner Earthquake Engineering Research Fund of the California Institute of Technology. The authors gratefully acknowledge this support. The first author is also partially supported by a grant from the National Natural Science Foundation of China (NSFC grant no. 51308161) and the International Postdoctoral Exchange Fellowship Program 2014 of the Office of China Postdoctoral Council (No. 20140018), and this support is also gratefully acknowledged.
Appendix A. MAP estimation of d Iterative equations for calculating the MAP estimate of d = [(x2)T, qT, sT, /T, g, m, aT, b, b0, j]T are obtained by differentiation of the logarithm of J1(d) defined in (21) with respect to d and then setting these derivatives to zero to find the stationary points of J(d) = log J1(d). Ignoring terms in the logarithm that are independent of d, the objective function is defined as:
T g ^ T ^ C/Þ 1 ðx ^ 2 Lx2 Þ E1 ðx ^ 2 L x2 Þ JðdÞ ¼ ðw C/Þ ðw 2 2 1 1 þ ðN d Nm Nh Þ log b log jHT Hj 2 2 1 T b T 1 T T ðb b b HðH HÞ H bÞ þ No Ns N m log g 2 2 Nm 1 1 X 1 þ Ns ðlog qi Þ log jA þ b1 ðHT HÞ j 2 i¼1 2 T 1 1 1 1 ðh^u ðHT HÞ HT bÞ ðA þ b1 ðHT HÞ Þ 2 Nm X 1 ð^hu ðHT HÞ HT bÞ þ log m mg þ ½log si si qi i¼1
1 b0 b
log b0 þ log j jb0
ðA1Þ
Following [39], we can simplify the following two terms in (A1): 1
log jA þ ðbHT HÞ j ¼ Nh logb log jHT Hj þ logjAj þ logjR1 j ðA2aÞ
52
Y. Huang et al. / Structural Safety 64 (2017) 37–53 1
T
1 1 1 ð^hu ðHT HÞ HT bÞ ðA þ ðbHT HÞ Þ ðh^u ðHT HÞ HT bÞ T
1
Then: NX NX m Nd m Nd @ðHT HÞ ¼ PTq H þ HT Pq ¼ PTq /s Ps þ /s PTs Pq @/q s¼1 s¼1 ! NX m Nd ¼ 2/q PTq Pq þ PTq /s Ps
1
¼ ð^hu ðHT HÞ HT bÞ ðbHT H bHT HðA1 þ bHT HÞ bHT HÞ 1
ð^hu ðHT HÞ HT bÞ T
1
1
¼ bðl ðHT HÞ HT bÞ HT Hðl ðHT HÞ HT bÞ T þð^hu lÞ A1 ð^hu lÞ
þ
1 1
where R = (bH H + A ) is the covariance matrix of h defined in (17b) and it is readily evaluated for arbitrary values of parameters /, a and b. Notice that (A2a) shows that the term with log |HTH| in (A1) cancels out. Now consider the sum of two of the terms in (A1): 1
T
@ðlogjR1 jÞ @ðHT HÞ ¼btr R @/q @/q
1
1
ðA3Þ
T
HT Hðl ðHT HÞ HT bÞ þ ð^hu lÞ A1 ð^hu lÞ T 2 ¼ bkHl bk þ ð^hu lÞ A1 ð^hu lÞ
Note that when differentiating (A3) with respect to uncertain parameters /, x2, b and a, the terms involving derivatives of l will drop out since:
bðdlÞ H ðHl bÞ þ bðHl bÞ Hdl þ ðdlÞ A ð^hu lÞ T
T
T
T
1
T
þð^hu lÞ A1 ðdlÞ ¼ ðdlÞ ðbHT H þ A1 Þl þ lT ðbHT H þ A1 Þdl bðdlÞ HT b T
T
ðA4Þ
bb Hdl dlT A1 ^hu ^hTu A1 dl T
T
T
¼ ðdlÞ R1 l þ lT R1 dl ðdlÞ R1 l lT R1 dl ¼ 0 Then for the term b||Hl b||2 in (A3), the derivatives with respect to /, x2 and b are expressed as:
@ðbkHl bk2 Þ=@/ ¼ @ðbkF/k2 Þ=@/ ¼ 2bFT F/ @ðbkHl bk2 Þ=@ x2 ¼ @ðbkGx2 ck2 Þ=@ x2 ¼ 2bðGT Gx2 GT cÞ
ðA5Þ where F is defined in (24) while G and c are defined in (29) and (30), respectively. The derivatives of log|R1| in (A2a) with respect to /, a and b are now derived. To differentiate log|R1| with respect to /, we first rewrite H as the following form:
6 H¼6 4
K1 /1 .. . K1 /Nm
3
KNh /1 7 .. .. 7 . . 5 KNh /Nm Nm N
¼
NX m Nd
/q Pq
ðA6Þ
q¼1 d N h
where /q are the NmNd components of / ¼ /T1 ; . . . ; /TNm and
@H Pq ¼ ; q ¼ 1; . . . ; Nm Nd @/q
as defined in (25).
!
ðA9Þ
Consequently, the derivative of log|R1| with respect to vector / is given by:
@ðlogjR1 jÞ ¼ 2b diag trðRT1 Þ; . . . ; trðRTNm Nd Þ / @/ T þ 2b trðRU1 Þ; . . . ; trðRUNm Nd Þ
ðA10Þ
Using the determinant derivative formula again, the derivatives of log|R1| with respect to aj and b are given by:
@ðlogjR1 jÞ @ðA1 Þ ¼ tr R @ aj @ aj
!
¼ a2 j Rjj
ðA11Þ
Nh X @ðlogjR1 jÞ ¼ trðRHT HÞ ¼ trðb1 RR1 b1 RA1 Þ ¼ b1 ð1 a1 j Rjj Þ @b j¼1
ðA12Þ where Rjj is the jth diagonal element of the posterior covariance matrix of R in (17b). Finally, setting the derivatives of J(d) in (A1) to zero leads to the update formulae for the parameters in d that are given in Section 3.2. We give some more details regarding (32) for the MAP estimation of equation error precision b, which is sensitive to the sparseness level of the identified stiffness reductions. Setting the ~ is estiderivative of J(d) with respect to b to zero, the MAP value b mated as:
@ðbkHl bk2 Þ=@b ¼ kHl bk2
2
¼ 2b/q trðRT q Þ þ 2btrðRU q Þ
¼ bðb b b HðHT HÞ HT bÞ þ bðl ðHT HÞ HT bÞT 1
Nm Nd s–q /s Ps
¼ 2b/q trðRTq Þ þ btrðRUq Þ þ btrðRUTq Þ
1
T
P
Therefore, using the Jacobi’s determinant derivative formula:
1
1 1 ðA þ b1 ðHT HÞ Þ ðh^u ðHT HÞ HT bÞ T
Pq ¼ 2/q Tq þ Uq þ UTq
where Tq ¼ PTq Pq and Uq ¼ PTq
bðb b b HðHT HÞ HT bÞ þ ð^hu ðHT HÞ HT bÞT T
! /s PTs
s–q
ðA2bÞ T
NX m Nd
ðA8Þ
s–q
T
2 RNm Nd 1
ðA7Þ
~¼ b
P h P h ~ ~ ~ 1 ~ 1 Nd N m Nj¼1 ð1 a N d Nm Nj¼1 ð1 a j Rjj Þ j Rjj Þ ¼ P 2 N 2 m ~ þ 2b ~1 ~1 ~l ~ i k þ 2b ~ 2 MÞ/ ~Þ x ~ bk kðKðl kH 0
i¼1
i
ðA13Þ
0
~0 and j ~ from the We continue by estimating the MAP values b derivatives of J(d) as follows:
~0 ¼ b
~ 2b qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ ~b 1 þ 1 þ 4j
j~ ¼ 1=b~0
ðA14Þ
ðA15Þ
~0 into (A 14), we get the MAP estimate ~ ¼ 1=b If we substitute j of b0 given in (33):
~0 ¼ b=2 ~ b
ðA16Þ
As a final point, we note that for real damage detection when the model has large modeling errors, it was found initially that during optimization for the MAP parameters, the components of
Y. Huang et al. / Structural Safety 64 (2017) 37–53
system mode shape / sometimes tended to zero along with the corresponding equation errors e, as mentioned in Remark 2.1. We found that to avoid this problem, we needed to carefully ~0 ~ and its associated scale parameterb choose the initial values of b in the iterative procedure in Algorithms 1 and 2. First, based on (A13) and (A16) but simplifying some terms in (A13):
~¼P b Ns r
4Ns No w ^ r;1 k2 ~ 21 MÞ kðK x
~0 ¼ P b Ns r
2Ns No w x ^ r;1 k2 ~ 21 MÞ kðK
ðA17Þ
ðA18Þ
^ r;1 2 RNo gives the identified mode shape of the 1st mode where w from the rth data segment. (A17) and (A18) are established by P m ~ i k2 =N m ~ 1 k2 ~ 2i MÞ/ ~ 21 MÞ/ replacing Ni¼1 kðK x withkðK x and P 2 N 2 m 2 ^ x ~ 21 MÞ/1 k =Nd with ~ kðK x kð K k =4N N . Here, K MÞ w r;1 s o i¼1 i are the reduced No No stiffness and mass matrices, respecand M tively, that correspond to the measured DOFs. Then, to achieve a good performance in the optimization for the MAP values in Algo~0 initially. Presumably, ~ ¼ 20b rithms 1 and 2 in Section 4, we take b
~ works because it is better if we make the larger initial value for b the model updating extract a larger amount of information from the system modal parameters x2 and / in the beginning, and there^ before pro^ 2 and w, fore from the ‘‘measured” modal parametersx ducing sparseness of the inferred stiffness reductions. References [1] Doebling SW, Farrar CR, Prime MB. A summary review of vibration-based damage identification methods. Shock Vib Digest 1998;30(2):91–105. [2] Vanik MW, Beck JL, Au SK. Bayesian probabilistic approach to structural health monitoring. J Eng Mech 2000;126(7):738–45. [3] Beck JL, Au SK, Vanik MW. Monitoring structural health using a probabilistic measure. Comput-Aid Civ Infrastruct Eng 2001;16(1):1–11. [4] Ou JP, Li H. Structural health monitoring in mainland China: review and future trends. Struct Health Monitoring 2010;9(3):219–31. [5] Farrar CR, Worden K. Structural health monitoring: A machine learning perspective. John Wiley & Sons; 2013. [6] Lynch JP. An overview of wireless structural health monitoring for civil structures. Phil Trans R Soc A 2007;365:345–72. [7] Li H, Huang Y, Ou J, Bao Y. Fractal dimension-based damage detection method for beams with a uniform cross-section. Comput-Aid Civ Infrastruct Eng 2011;26(3):190–206. [8] Moaveni B, Stavridis A, Lombaert G, Conte JP, Shing PB. Finite-element model updating for assessment of progressive damage in a 3-story infilled RC frame. J Struct Eng 2013;139:1665–74. [9] Bernal D. Load vectors for damage localization. J Eng Mech 2002;128(1):7–14. [10] Vanik MW. A Bayesian probabilistic approach to structural health monitoring, Technical Report EERL 97-07. Caltech, Pasadena, CA: Earthquake Engineering Research Laboratory; 1997. [11] Sohn H, Law KH. A Bayesian probabilistic approach for structure damage detection. Earthquake Eng Struct Dyn 1997;26(12):1259–81. [12] Ching J, Beck JL. Two-step Bayesian structure health monitoring approach for IASC-ASCE phase II simulated and experimental benchmark studies, Technical Report EERL 2003-02. California Institute of Technology, Pasadena, CA: Earthquake Engineering Research Laboratory; 2003. [13] Yuen KV, Au SK, Beck JL. Two-stage structural health monitoring approach for phase I benchmark studies. J Eng Mech 2004;130(2):16–33. [14] Ching J, Beck JL. Bayesian analysis of the phase II IASC–ASCE structural health monitoring experimental benchmark data. J Eng Mech 2004;130(10):1233–44. [15] Yuen KV, Beck JL, Katafygiotis LS. Efficient model updating and health monitoring methodology using incomplete modal data without mode matching. Struct Control Health Monitoring 2006;13(1):91–107. [16] Yuen KV. Bayesian methods for structural dynamics and civil engineering. Wiley; 2010.
53
[17] Yin T, Lam HF, Chow HM. A Bayesian probabilistic approach for crack characterization in plate structures. Comput-Aid Civ Infrastruct Eng 2010;25 (5):375–86. [18] Sankararaman S, Mahadevan S. Bayesian methodology for diagnosis uncertainty quantification and health monitoring. Struct Control Health Monitoring 2013;19:88–106. [19] Figueiredo E, Radub L, Worden K, Farrar CR. A Bayesian approach based on a Markov-chain Monte Carlo method for damage detection under unknown sources of variability. Eng Struct 2014;80(1):1–10. [20] Arangio S, Bontempi F. Structural health monitoring of a cable-stayed bridge with Bayesian neural networks. Struct Infrastruct Eng 2015;11(4):575–87. [21] Beck JL. Bayesian system identification based on probability logic. Struct Control Health Monitoring 2010;17(7):825–47. [22] Girolami M, Calderhead B, Vyshemirsky V. System identification and model ranking: The Bayesian perspective. In: Learning and Inference in Computational Systems Biology. MIT Press; 2010. [23] Green PL, Cross EJ, Worden K. Bayesian system identification of dynamical systems using highly informative training data. Mech Syst Sig Process 2015;56–57:109–22. [24] Worden K, Hensman JJ. Parameter estimation and model selection for a class of hysteretic systems using Bayesian inference. Mech Syst Sig Process 2012;32:153–69. [25] Au SK, Zhang FL. Fundamental two-stage formulation for Bayesian system identification. Part I: General theory. Mech Syst Sig Process 2016;66– 67:31–42. [26] Zhang FL, Au SK. Fundamental two-stage formulation for Bayesian system identification. Part II: Application to ambient vibration data. Mech Syst Sig Process 2016;66–67:43–61. [27] Cox RT. Probability, frequency and reasonable expectation. Am J Phys 1946;14 (1):1–13. [28] Jaynes ET. Information theory and statistical mechanics. Phys Rev 1957;106:620–30. [29] Jaynes ET. Probability theory: the logic of science. Cambridge, U.K.: Cambridge University Press; 2003. [30] Zhou H, Ni Y, Ko J. Eliminating temperature effect in vibration-based structural damage detection. J Eng Mech 2011;137(12):785–96. [31] Sohn H, Allen DW, Worden K, Farrar CR. Structural damage classification using extreme value statistics. J Dyn Sys Meas Control 2005;127(1):125–32. [32] Tipping ME. Sparse Bayesian learning and the relevance vector machine. J Mach Learning Res 2001;1:211–44. [33] Wipf D, Palmer J, Rao B. Perspectives on sparse Bayesian learning. Adv Neural Inf Process Syst 2004;16. [34] Ji S, Xue Y, Carin L. Bayesian compressive sensing. IEEE Trans Sig Process 2008;56(6):2346–56. [35] Huang Y, Beck JL, Wu S, Li H. Bayesian compressive sensing for approximately sparse signals and application to structural health monitoring signals for data loss recovery. Probabilistic Eng Mec. 2016. doi: http://dx.doi.org/10.1016/j. probengmech.2016.08.001. [36] Huang Y, Beck JL, Wu S, Li H. Robust Bayesian compressive sensing for signals in structural health monitoring. Comput-Aid Civ Infrastruct Eng 2014;29 (3):160–79. [37] Mackay DJC. Bayesian methods for adaptive models (Ph.D. Thesis in Computation and Neural Systems). Pasadena, California: California Institute of Technology; 1992. [38] Huang Y, Beck JL. Novel sparse Bayesian learning for structural health monitoring using incomplete modal data. In: Proceedings of 2013 ASCE International Workshop on Computing in Civil Engineering, Los Angeles, CA; June 2013. [39] Huang Y, Beck JL. Hierarchical sparse Bayesian learning for structural health monitoring with incomplete modal data. Int J Uncertain Quantif 2015;5 (2):139–69. [40] Beck JL, Katafygiotis LS. Updating models and their uncertainties. I: Bayesian statistical framework. J Eng Mech 1998;124(4):455–61. [41] Bernal D, Dyke SJ, Lam HF, Beck JL. Phase II of the ASCE benchmark study on SHM. In: Proceedings of 15th eng. mechanics conf., ASCE; 2002. p. 1048–55. [42] Dyke SJ, Bernal D, Beck JL, Ventura C. Experimental Phase II of the Structural Health Monitoring Benchmark Problem. In: Proceedings of 16th Eng. Mechanics Conf., ASCE, Reston, VA; 2003. [43] Beck JL, Jennings PC. Structural identification using linear models and earthquake records. Earthquake Eng Struct Dyn 1980;8(2):145–60. [44] Beck JL, May BS, Polidori DC. Determination of modal parameters from ambient vibration data for structural health monitoring. In: Proceedings of first world conference on structural control, Pasadena, California; 1994. p. TA3:3–TA3:12.
Contents lists available at ScienceDirect
Structural Safety journal homepage: www.elsevier.com/locate/strusafe
Hierarchical sparse Bayesian learning for structural damage detection: Theory, computation and application Yong Huang a,b, James L. Beck b,⇑, Hui Li a a b
Key Lab of Structural Dynamic Behavior and Control of the Ministry of Education, School of Civil Engineering, Harbin Institute of Technology, Harbin, China Division of Engineering and Applied Science, California Institute of Technology, CA, USA
a r t i c l e
i n f o
Article history: Received 10 March 2016 Received in revised form 2 September 2016 Accepted 11 September 2016 Available online 31 October 2016 Keywords: Inverse problems Hierarchical sparse Bayesian learning Model updating Structural health monitoring Damage detection Bayesian Ockham razor IASC-ASCE benchmark problems
a b s t r a c t Structural damage due to excessive loading or environmental degradation typically occurs in localized areas (in the absence of collapse) where it leads to local stiffness reductions. This prior information about the spatial sparseness of structural damage and the associated stiffness loss is exploited here by a hierarchical sparse Bayesian learning (SBL) framework, with the goal to reduce the ill-conditioning in the stiffness loss inversion problem for damage detection. We have previously proposed a SBL approach to establish the probability of localized stiffness reductions caused by damage by using noisy incomplete modal data from before and after possible damage. The excellent performance achieved by introducing sparseness in the damage pattern was demonstrated by using synthetic data where there are only small modeling errors. In this research, a more rigorous formulation along with a corresponding efficient and scale-invariant SBL algorithm are developed. The algorithm is first applied to synthetic data, then to real vibration response from a steel-frame test structure where there are large modeling errors. These data are from the Phase II simulated and experimental benchmark studies that were sponsored by the IASC-ASCE Task Group on Structural Health Monitoring. The results show that, even for the real data, the proposed method can reliably detect, locate and assess damage of the benchmark structure by inferring substructure stiffness losses using the identified modal parameters from the calibration and monitoring stages. The occurrence of missed and false damage alerts is effectively suppressed, and we show that the new algorithm gives a better performance than our previous SBL method in the real data case where there is significant modeling error. Several appealing features of our method are summarized at the end of the paper. Ó 2016 Published by Elsevier Ltd.
1. Introduction The challenge of developing automated sensor-based structural health monitoring (SHM) systems for accurately detecting, locating and assessing damage from severe loading events (such as earthquakes, hurricanes or explosions) or progressive structural deterioration (such as fatigue) has received much attention over the last several decades. Numerous vibration-based global SHM techniques have been developed for estimating both damage location and damage severity [e.g. 1–9], most of which utilize the dependence of the identified structural modal parameters, such as natural frequencies and mode shapes, on physical properties of structures, i.e., stiffness, mass and damping. In practice, most SHM systems record structural vibration data at low-amplitude excitation, e.g., ambient vibration from wind and ground motion, and so a linear ⇑ Corresponding author. E-mail address: [email protected] (J.L. Beck). http://dx.doi.org/10.1016/j.strusafe.2016.09.001 0167-4730/Ó 2016 Published by Elsevier Ltd.
dynamic model with classical normal modes is adequate for damage detection purposes. The typical strategy involves comparing structural models identified from sets of measured modal data (natural frequencies and mode shapes) from a structure before and after possible damage. For applications in structural health monitoring, there are some inherent difficulties in this modelbased approach. One of the main difficulties is that it is impossible to exactly model the full behavior of a structure by using the limited knowledge available. For damage inference problems, there are always modeling uncertainties involved and hence the need to account for uncertainty in model updating arises. These uncertainties may arise from many sources: simplifying approximations to develop the structural models; limited number of sensors installed on a structure; sensor noise; thermally-induced daily variations in structural stiffness, etc. This presence of substantial modeling uncertainty has motivated numerous researchers to tackle the problem of model updating in structural system identification
38
Y. Huang et al. / Structural Safety 64 (2017) 37–53
Nomenclature Nm No Nd Ns Nh M, K h /r ; x2r b ^ r;i ^ 2r;i ; w x
number of extracted modes in the modal identification we use (i = 1, . . ., Nm) number of observed degrees of freedom number of degrees of freedom of the structural model for identification we use (k = 1, . . ., Nd) number of time segments of measured modal data we use (r = 1, . . ., Ns) number of substructures considered (j = 1, . . ., Nh) mass and stiffness matrices of the structural model for identification structural stiffness scaling parameters system mode shape and system natural frequencies of the rth mode equation error precision parameter MAP estimates of modal frequency and mode shape of the ith mode from the rth data segment from modal identification
and damage assessment from a Bayesian perspective [2,3,10–26]. Rather than considering only a point estimate as in the conventional deterministic or frequentist probabilistic methods, Bayesian methods use Bayes’ Theorem to quantify the relative plausibility of all possible values of the structural model parameters via their Bayesian posterior PDF (probability density function). The relative plausibility of each model within a set of candidate model classes chosen to represent the uncertain structural behavior can also be assessed by applying Bayes’ Theorem at the model class level, enabling the consideration of all plausible models and the construction of robust predictions of future structural behavior [21]. We emphasize that probability in the Bayesian sense should be interpreted as the degree of plausibility of a statement on the basis of the specific conditioning information [21,27–29], where the statement may refer to the structural parameters or to the structural model itself. It quantifies uncertainty due to our incomplete information and does not need the postulate of ‘‘inherent randomness”. Another fundamental difficulty for damage detection lies in the fact that the mathematical nature of the underdetermined structural inverse problem based on the constraints imposed by the available data requires a tradeoff between the resolution of the damage locations and the reliability of the probabilisticallyinferred damaged state. In this case, defining a proper threshold to determine whether the damage features shift from their healthy state is very important in order to alleviate false positive (false alarm) and false negative (missed alarm) detections, but it is very challenging to establish a reliable threshold value in a rigorous manner for issuing a timely damage alarm [30,31]. In this paper, we aim to extend the applicability of Bayesian methods to produce reliable damage alarms even for highdimensional model parameter spaces (higher-resolution damage localization), by the exploration of recent developments in sparse Bayesian learning (SBL) [32,33] and Bayesian compressive sensing/sampling [34–36]. This new perspective for the stiffness loss inverse problem exploits the prior knowledge that structural damage typically occurs at a limited number of locations in a structure in the absence of its collapse. A specific hierarchical Bayesian model and the corresponding evidence maximization strategy for Bayesian inference is utilized to promote spatial sparseness in the inferred stiffness reductions in a way that is consistent with the Bayesian Ockham razor [21,37]. We have proposed previously to use a sparse Bayesian perspective to tackle the stiffness loss inverse problem [38,39]. In this
^hu
C T
a g, q m, s, j d
MAP estimate of h determined from the calibration test data. matrix that picks the measured degrees of freedom from each system mode shape matrix relating the vector of Ns sets of identified natural ^ 2 and the system natural frequencies x2 frequencies x variance parameter for the likelihood function of structural stiffness scaling parameters h measurement-error precision parameters for mode shapes and natural frequencies rate parameters controlling the exponential prior distributions of g, q, b0 set of all hyper-parameters x2, q, s, /, g, m, f, b, a0, b0, j
paper, we present an improved version of our SBL method that eliminates two possibly unreliable approximations in our original theoretical formulation. We also propose effective strategies to alleviate the ill-conditioning for stiffness inversions when the modeling error is large, which enhance its applicability and accuracy for damage detection and assessment with real data. In addition, the new algorithm allows more efficient computation and is independent of the unit selections of mass and stiffness matrices, as well as the scaling of the measured mode shapes. We also apply our method to analyze the IASC–ASCE Phase II experimental benchmark data for both brace and joint damage cases, not only the synthetic data that we used previously [38,39]. For both the simulated and experimental data, we show that the performance of the new method is better than our previous method in [39], especially for the joint damage cases. The remainder of the paper is organized as follows. In Section 2, we present our hierarchical Bayesian model, which is similar to that in [39] but which has several modifications. The Bayesian model is then used in an evidence maximization procedure in Section 3, yielding an improved Bayesian inference framework that lends itself to reliable, efficient and scale-invariant computation of the posterior PDF for the structural stiffness parameters. In Section 4, the algorithms for stiffness identification in the calibration and monitoring stages are presented. We then compare in Section 5 our new SBL damage assessment algorithm with that in [39] in terms of the theory and computation. Applications of our Bayesian methods to simulated and real structural data are presented in Section 6, and conclusions and future work are discussed in Section 7. 2. Hierarchical Bayesian model class Suppose that Ns sets of measured vibration time histories are available from a structure and Nm modes of vibration have been identified for each set of time histories so that we have a vector of h iT ^2¼ x ^ 21;1 ; ...; x ^ 21;Nm ; x ^ 22;1 ; ... ; x ^ 2Ns ;Nm identified natural frequencies: x 2 RNs Nm 1 and a vector of identified mode shapes: h iT N N N 1 T T T T s m o ^¼ w ^ ;...;w ^ ^ ^ ^ r;i 2 RNo w 2R , where w 1;1 1;Nm ; w2;1 ; . . . ; wNs ;Nm gives the identified components of the system mode shape of the ith mode (i = 1, . . ., Nm) at the No measured DOF (degrees of freedom) from the rth data segment (r = 1, . . ., Ns). The No measured DOF in a structure are usually a smaller subset of the Nd DOF of an appropriate structural model. It has been
39
Y. Huang et al. / Structural Safety 64 (2017) 37–53
shown to be advantageous [3,10,12–15] to introduce the system T frequencies x2 ¼ x21 ; . . . ; x2Nm 2 RNm 1 and system T T mode shapes / ¼ /1 ; . . . ; /TNm 2 RNd Nm 1 to represent the actual underlying modal parameters of the linear dynamics of the structural system at all DOFs corresponding to those of the structural model. The system modal parameters are therefore distinct from the corresponding modal parameters of the structural model. Based on a defined stochastic model class M (see ^ to update ^ 2 and w [21]), one can use the available modal data x
natural
the system modal parameters x2 and / together with structural model parameters h, without using mode matching [2,3,10,12– 15]. 2.1. Structural model class
As part of the definition of the stochastic model class M, we choose a set of linear structural models with Nd DOFs and with classical normal modes (so a damping matrix need not be explicitly modeled) where each model has the same known mass matrixM 2 RNd Nd and an uncertain stiffness matrix K parameterized as a linear combination of (Nh + 1) substructure stiffness matrices Kj, j = 0, 1, . . ., Nh: Nh X KðhÞ ¼ K0 þ hj Kj
ð1Þ
j¼1
where K0 is any part of the global stiffness matrix that is considered known and that corresponds to parts of the structure that are believed to be not vulnerable to damage, and Kj 2 RNd Nd ; j ¼ 1; . . . ; N h ; is the substructure stiffness matrix representing the nominal contribution of the jth substructure to the global stiffness matrix K from a structural model (e.g., specified by a finite element model of a structure). The corresponding stiffness scaling parameter hj, j = 1, . . ., Nh, is a factor that allows modification of the nominal jth substructure stiffness so as to be more consistent with the real structure behavior and any reduction of hj, j = 1, . . .Nh, is assumed to correspond to damage in the jth substructure. The parameter vectorh ¼ ½h1 ; . . . ; hn 2 RNh ; therefore specifies the structural model in M: 2.2. Prior distribution for structural model parameters and system modal parameters We do not assume that the Nm system natural frequencies x2 and Nm system mode shapes / satisfy the eigenvalue problem corresponding to any structural model specified by h because there will always be modeling errors e ¼ eT1 ; . . . ; eTNm 2 RNm Nd ; where ðKðhÞ x2i MÞ/i ¼ ei ; i ¼ 1; . . . ; N m ; for the ith system mode and the structural model specified by h. The uncertain eigenvalue equation error e therefore provides a bridge between the behaviors of the real system and a deterministic structural model. As described in [39], the following prior PDF conditional on the equation-error precision hyper-parameter b is chosen:
pðx2 ; /; hjbÞ ¼ c0 ð2p=bÞNm Nd =2 ( ) Nm bX 2 2 exp kðKðhÞ xi MÞ/i k 2 i¼1
can deduce the prior PDF for h conditional on system modal parameters [x2, /] as: 1
1
pðhjx2 ;/;bÞ ¼ pðx2 ;/;hjbÞ=pðx2 ;/jbÞ ¼ NðhjðHT HÞ HT b;ðbHT HÞ Þ ¼ ð2p=bÞNh =2 jHT Hj1=2 T 1 1 b exp ðh ðHT HÞ HT bÞ HT Hðh ðHT HÞ HT bÞ 2 ð3Þ where the marginal prior PDF for the system modal parameters [x2, /] and the parameter vector b and matrix H are defined by:
pðx2 ; /jbÞ ¼ c0 ð2p=bÞðNh Nd Nm Þ=2 jHT Hj1=2 1 b T T exp ðb b b HðHT HÞ HT bÞ 2 3 ðx21 M K0 Þ/1 6 7 .. 7 b¼6 . 4 5 ðx2Nm M K0 Þ/Nm N
ð4Þ
2
2
K1 /1 6 . H¼6 4 .. K1 /Nm
ð5Þ m N d 1
3 KNh /1 7 .. 7 . 5 KNh /Nm Nm N .. .
ð6Þ d N h
We model the uncertainty in b by the following Exponential hyper-prior: 1
1
pðbjb0 Þ ¼ b0 expðb0 bÞ
ð7Þ
where b0 is the scale parameter of the distribution. The Exponential prior for b is the maximum entropy distribution with the constraint of E(b|b0) = b0 and is a special case of the widely used Gamma prior. To penalize values of b0 that are too large, we assign an Exponential hyper-prior for the scale parameter b0 as the last stage of the hierarchical model:
pðb0 jjÞ ¼ Expðb0 jjÞ ¼ j expðjb0 Þ
ð8Þ
Remark 2.1. The joint prior in (2) is defined based on the Gaussian probability model for uncertain equation errors e, which is the maximum entropy PDF subject to the first two moments as constraints. Those values of {x2, /, h} that make the equation errors e smaller in norm are considered more plausible a priori; i.e. they have larger values for their prior PDF. However, the joint PDF in (2) has an undesirable peak where all the components of the system mode shapes / are zeros since the corresponding equation errors e are then also zero. For real damage detection when the model is expected to have large modeling errors, the scale parameter b0 is important for scaling the equation-error precision parameter b appropriately to prevent equation error e, and so the components of /, becoming zero. This point is discussed further at the end of the Appendix. 2.3. Likelihood function for structural model parameters and system modal parameters
ð2Þ
where c0 is a normalizing constant and |||| denotes the Euclidean vector norm. Note that the finite value of the equation-error precision parameter b in (2) provides a soft constraint for the eigenequation and it allows for the explicit control of how closely the system and model modal parameters agree. Following [39], we
Since the onset of damage is typically in a small number of locations in a structure in the absence of structural collapse, the poten^u can tial change in the parameter vector from damage Dh ¼ h h be considered as a sparse vector with relatively few non-zero components, where h and hu are the stiffness scaling parameters for the current state (possibly damaged) and the calibration (undamaged) state, and ^ hu is the unique MAP estimate of hu from applying Baye-
40
Y. Huang et al. / Structural Safety 64 (2017) 37–53
sian updating using a large amount of time-domain vibration data from the undamaged structure. Following [39], we choose the MAP value ^ hu from the calibration state as pseudo-data for h to define a likelihood function as:
pð^hu jh; aÞ ¼ Nð^hu jh; AÞ ¼
Nh Y Nð^hu;j jhj ; aj Þ
ð9Þ
j¼1
where A ¼ diagða1 ; . . . ; aNh Þ. The idea here is to learn each aj from hu;j , which is interpreted the modal data and if aj ? 0, then hj ! ^ as the jth substructure being undamaged. The choice of likelihood function in (9) is motivated by the closely-related SBL framework. Although the conventional strategy in SBL is to use an automatic relevance determination (ARD) Gaussian prior PDF [37] to model sparseness, here we incorporate the ARD concept in the likelihood function along with the prior on h in (3). Following the same approach as in [39], the likelihood function for the system modal parameters x2 and / is:
^ x2 ; /; hÞ ¼ pðx ^ ^ ^ 2 ; wj ^ 2 jx2 Þpðwj/Þ ^ 2 jLx2 ; EÞNðwjC/; pðx ¼ Nðx CÞ ð10Þ where the selection matrixC 2 RNo Nm Ns Nm Nd with ‘‘1s” and ‘‘0s” picks ^ the observed DOF of the whole ‘‘measured” mode shape data set w from the system mode shapes /; L ¼ ½INm ; . . . ; INm T 2 RN0 Nm Nm is the matrix between the vector of Ns sets of Nm identified natural fre^ 2 and the Nm system natural frequencies x2; quencies x E ¼ block diagðE1 ; . . . ; ENs Þ is a block diagonal covariance matrix 1 with the diagonal block Er ¼ diag q1 1 ; . . . ; qN m ; r ¼ 1; . . . ; N s ; C ¼ g1 INo Ns Nm ; INm andINo Ns Nm denote the identity matrices of
corresponding size; and q ¼ ½q1 ; . . . ; qNm T and g are prescribed precision parameters for the predictions of the identified natural fre^ from the system modal ^ 2 and mode shapes w quencies x parameters, respectively. In a hierarchical manner, exponential priors are placed on the parameters qi and g:
pðqi jsi Þ ¼ Expðqi jsi Þ ¼ si expðsi qi Þ; i ¼ 1; . . . ; Nm pðgjmÞ ¼ ExpðgjmÞ ¼ m expðmgÞ
ð11Þ
which are the maximum entropy priors with support [0, 1) for and m1 of qi and g, respectively. Then the given mean values s1 i prior PDF for the parameter vector q is given by:
pðqjsÞ ¼
Nm Nm Nm Y Y X pðqi jsi Þ ¼ si exp si qi i¼1
i¼1
!
ð12Þ
i¼1
where s ¼ ½s1 ; . . . ; sNm . T
2.4. Joint PDF for hierarchical Bayesian model By combining all stages of the hierarchical Bayesian model, the joint PDF for all of the stochastic variables is:
^ ^hu Þ ¼ ^ 2 ; w; pðx2 ; q; s; /; g; m; h; a; b; b0 ; j; x 2 2 ^ ^ ^ jx ; qÞpðwj/; gÞpðhu jh; aÞ pðx pðx2 ; /; hjbÞpðqjsÞpðgjmÞpðaÞpðsÞpðmÞpðbjb0 Þpðb0 jjÞpðjÞ
ð13Þ
^ ^ 2 jx2 ; qÞpðwj/; where PDFs pðx gÞ; pð^hu jh; aÞ; pðx2 ; /; hjbÞ; pðqjsÞ; pðgjmÞ; pðbjb0 Þ and p(b0|j) are defined in (10), (9), (2), (12), (11), (7) and (8), respectively. The prior PDFs for aj, si, m and j are taken as uniform over large open intervals that start at zero. The acyclic graph of the hierarchical Bayesian model is shown in Fig. 1, where each arrow denotes the conditional dependencies used in the joint probability model. The bidirectional arrow in the graph of the hierarchical Bayesian model represents the statistical interaction between the system modal parameters x2 and /, which comes from
Fig. 1. Acyclic graph representing the information flow in the hierarchical Bayesian model.
the joint prior p(x2, /|b) in Eq. (4). Compared with non-hierarchical Bayesian models where b, q and g are selected a priori [15], our hierarchical Bayesian model allows all sources of uncertainty and correlation to be learned from the data, and hence potentially produces more accurate damage identification. Remark 2.2. The hierarchical Bayesian model in Fig. 1 is similar to that of [39]; however, the exponential hyper-prior PDF pðajkÞ for a is replaced by a simple uniform distribution since the learning of this hyper-prior using real-data with large modeling errors may produce over-sparse models with too few stiffness reductions. In addition, an exponential hyper-prior is assigned for scale parameter b0 for appropriately scaling the prior distribution of the equation-error precision parameter b when learning from data where the model is expected to have large modeling errors.
3. Bayesian inference To facilitate the goal for a rigorous SBL formulation, we collect all uncertain parameters except stiffness scaling parameter h in the vector d = [(x2)T, qT, sT, /T, g, m, aT, b, b0, j]T and the full posterior uncertainty in h is explicitly incorporated when finding the MAP estimates of all parameters in d using the evidence maximization procedure. 3.1. Approximation of the full posterior PDF using Laplace’s method As widely known, Bayesian inference is based on the full poste^ ^ ^ 2 ; w; rior PDF over all uncertain parameters pðh; djx hu Þ which is given by Bayes’ theorem, but which is almost always analytically intractable because of the high-dimensional normalization integral. Nevertheless, we can turn to a hierarchical Bayesian procedure combined with Laplace’s asymptotic approximation [40] to find an effective approximation of the full posterior PDF ^ ^ ^ 2 ; w; pðh; djx hu Þ: Using the probability product rule, we rewrite the full posterior as:
^ ^hu Þ ¼ pðhjd; x ^ ^hu Þpðdjx ^ ^hu Þ ^ 2 ; w; ^ 2 ; w; ^ 2 ; w; pðh; djx
ð14Þ
The first factor in (14) is the posterior PDF of the stiffness scaling parameter vector conditional on d, which is given by:
41
Y. Huang et al. / Structural Safety 64 (2017) 37–53
^ ^hu Þ ¼ pðhjd; ^hu Þ ¼ pð^hu jh; dÞpðhjdÞ=pð^hu jdÞ ^ 2 ; w; pðhjd; x
ð15Þ
^ ^ 2 and w where the first equality holds because h is independent of x when conditional on x2 and / in d (see Fig. 1) and the second equality is due to Bayes’ Theorem. As a consequence of combining a Gaussian prior (3) and a linear model within a Gaussian likelihood (9), the conditional posterior is also Gaussian:
^ ^hu Þ ¼ Nðhjl; RÞ ^ 2 ; w; pðhjd; x
ð16Þ
stituting the MAP estimates of m into that of g and letting all other ~ g ~ parameters be fixed at their MAP values, the MAP estimates /; ~ and m are derived as (using (A3), (A5) and (A10) in (A1) in the Appendix):
1 ~T F ~ þ bdiag ~¼ b ~F ~ ~ Nm N Þ þ g ~ 1 Þ; . . . ; trðRT ~ CT C / trðRT d ^ b ~ trðRU ~ 1 Þ; . . . ; trðRU ~ Nm N Þ T ~ CT w g d
ð22Þ
with mean and covariance matrix:
N N N 2
l ¼ RðbHT b þ A1 ^hu Þ
ð17aÞ
1 R ¼ bHT H þ A1
g~ ¼ s^ m o~ 2 ¼ ~ m kw C/k
ð17bÞ
~ a ~ in H in (6), ~ is given by (17b) using b, ~ and also / where R
The normalizing constant for the posterior PDF in (15) is called the pseudo-evidence (or pseudo-marginal likelihood) for the model classMðdÞ given by pseudo-data^ hu and it can be analytically evaluated:
pð^hu jdÞ ¼
Z
pð^hu jh; dÞpðhjdÞdh ¼
Z
T
pð^hu jh; aÞpðhjx ; /; bÞdh ð18Þ
1
The second factor in (14) is given by Bayes’ Theorem:
ð19Þ
We assume that the model class is globally identifiable based on ^ h ^u jdÞ in ^ 2 ; w; the modal data, meaning here that the likelihood pðx (19) has a unique global maximum over d and then so does the posterior at ~ d (the MAP value of d). We then treat d in (14) as a ‘nuisance’ parameter vector and integrate it out using Laplace’s asymptotic approximation [40]:
^ ^hu Þpðdjx ^ ^hu Þdd pðhj~d; x ^ ^hu Þ ^ 2 ; w; ^ 2 ; w; ^ 2 ; w; pðhjd; x ð20Þ
^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; where ~ d ¼ arg max pðdjx hu Þ ¼ arg max½pðx hu jdÞ pðdÞ because the denominator in (19) is independent of d. In the next ^ ^ ^ 2 ; w; section, we will focus on the task of maximizing pðdjx hu Þ for ~ finding the MAP value d. ^ ^ 3.2. MAP estimate of d from PDF pðdjx ; w; hu Þ ^2
Finding the MAP value ~ d is equivalent to maximizing log J1(d) with respect to d where:
^ ^hu jdÞpðdÞ ¼ pðx ^ gÞpð^hu jx2 ; /; b; aÞ ^ 2 ; w; ^ 2 jx2 ; qÞpðwj/; J 1 ðdÞ ¼ pðx pðx ; /jbÞpðqjsÞpðgjmÞpðbjb0 Þpðb0 jjÞ 2
3
0 .. .
7 7 5
~Þ x Kðl
~ 21 M
ð24Þ
N m N d Nm Nd
~ a ~ x ~ 2 in b and H in (5) ~ is given by (17a) using b, ~ and also /, where l and (6), and for q = 1, . . ., NmNd:
Tq ¼ PTq Pq
Uq ¼ PTq
NX m Nd
! ~ s Ps /
s–q
^2 ^ ^ ^ ^hu Þ ¼ pðx ; w; hu jdÞpðdÞ pðdjx ; w; ^ h^u Þ ^ 2 ; w; pðx ^2
Z
~Þ x ~ 21 M Kðl .. .. . . 0
where D = A + b (H H) and we have utilized the hierarchical structure exhibited in Fig. 1. This is a Gaussian distribution over the Nh-dimensional pre-damage parameter vector hu evaluated at ^u , and is readily evaluated for arbitrary values of d. its MAP value h
^ ^hu Þ ¼ ^ 2 ; w; pðhjx
6 ~¼6 F 4
ð23Þ
2
1 ¼ Nð^hu jðHT HÞ HT b; DÞ 1
2
1
ð21Þ
where we have used the hierarchical structure exhibited in Fig. 1. Substituting equations (10), (18), (4), (12), (11), (7) and (8) into (21), it is shown in the Appendix that explicit expressions can be obtained for iterative maximization of log J1(d) where we update all of the parameters in d successively and repeat until convergence. By maximizing log J1(d) with respect to the system mode shapes / and the associated hyper-parameters g and m successively, sub-
Pq ¼
@H @/q
ð25Þ
~ requires repeated evaluation of the Determining the MAP value/ ~T F ~þg ~F ~ ~ 1 Þ; . . . ; trðRT ~ Nm N ÞÞÞ for dif~ CT C þ bdiagðtrð inverse of X ¼ ðb RT d ferent trial values for the other parameters. The dimension NmNd of / will often be large in applications (e.g., NmNd > 100). To facilitate the goal of efficient computation of the inverse of the matrix X, we rewrite it in the block diagonal form: X ¼ diagðX1 . . . ; XNm Þ; where the Nd Nd matrices Xi for i = 1, . . ., Nm are given by: ~ ~ ~ ði1ÞN þ1 Þ;. . . ;trðRT ~ iN ÞÞ þ g ~ 2 MÞ2 þ bdiagðtrð ~ N m CT Cs ; Xi ¼ bðKð l~ Þ x RT i
~ is given by: and then the MAP value /
d
d
s
~ ¼ diag X1 ; . . . X1 g ^ b ~ trðRU ~ 1 Þ; . . . ; trðRU ~ Nm N Þ T ~ CT w / 1 Nm d ð26Þ where Cs is the sub-matrix which contains the elements of the first ~ now involves Ns rows and Nd columns of C. The estimation of / inversions of Nm square matrices Xi (i = 1, . . ., Nm) of order Nd, ~ in (22) involves an inversion whereas the original estimation of / of one square matrix of order (NmNd). The latter requires a factor of N2m more computational effort than (26). Similarly, the MAP estimates of the system natural frequencies x2 and their associated hyper-parameters q and s are also derived as (using (A3), (A5) in (A1) in the Appendix):
h
i1 ~ T ~c ~ 1 x ~G ^2 þb LT E
~T G ~ ~ 1 L þ b ~G ~ 2 ¼ LT E x
and for each i = 1, . . ., Nm:
ð27Þ
42
Y. Huang et al. / Structural Safety 64 (2017) 37–53
Ns 2
q~ i ¼ P Ns r¼1
1 2 ¼ ~ si 2
ð28Þ
^ 2r;i x ~i x
6 ~¼6 G 4
~1 M/ .. . 0
3
0 .. .. 7 7 . . 5 ~N M/ m N d N m Nm
ð29Þ
ð30Þ
N d N m 1
Similarly, the MAP estimates of the hyper~ j ; j ¼ 1; . . . ; N h ; are derived as (using (A3) and (A11) parametersa in (A1) in the Appendix): 2 ~ þ ð^hu l ~ Þj a~ j ¼ ðRÞ jj
ð31Þ
~0 and j ~ b ~ are obtained by taking the Finally, the MAP values b, derivatives of log J1(d) with respect to b, b0 and j and substituting the MAP estimates of j into that of b and b0 (see (A13), (A15) and (A16) in the Appendix):
~¼ b
Nd Nm
PNh j¼1
~ ~ 1 1a j Rjj
~ þ 2b ~1 ~l ~ bk kH 0
~0 ¼ b=2 ~ ¼ 1=j ~ b
¼
log½pð^hu jx2 ; /; b; aÞpðhj^hu ; x2 ; /; b; aÞdh log½pðhjx2 ; /; bÞpð^hu jh; aÞ= pðhj^hu ; x2 ; /; b; aÞpðhj^hu ; x2 ; /; b; aÞdh log½pð^hu jh; aÞpðhj^hu ; x2 ; /; b; aÞdh Z log½pðhj^hu ; x2 ; /; b; aÞ=pðhjx2 ; /; bÞ pðhj^hu ; x2 ; /; b; aÞdh
~1 ~ Þ/ Kðl 7 6 . 7 ~c ¼ 6 .. 5 4 ~N ~ Þ/ Kðl m
Z
Z
3
2
Z
¼
where:
2
log½pð^hu jdÞ ¼
2
ð32Þ
ð33Þ
~ is inversely dependent on the A nice feature of (32) is that b square of the scaling of the mode shape and the units of the stiffness and mass matrices, which leads to the desirable property that the inferred stiffness reduction results are scale-invariant, since the terms bHTH and bHTb in (17a) and (17b) are independent of these scale selections. The MAP estimate ~ d can be calculated by going through the sequence of equations (26), (23), (27), (28), (31), (32) and (33) and then updating the Gaussian posterior in (20) for h by substituting ~ d into (17a,b) before repeating this process. For damage detection with real data where there are large modeling errors, it was sometimes observed during iterative solution for the MAP estimates ~ d that the components of the estimated sys~ kept decreasing towards zero while the tem mode shape / ~ tended to infinity. To avoid equation-error precision parameterb ~0 for the iterative procedure this occurring, a good initial choice of b
~0 is is required, which is discussed at the end of the Appendix, and b updated separately in an outer loop until convergence of the optimization occurs for the other parameters in d. 3.3. Information-theoretic interpretation of producing sparseness in the inferred stiffness reductions
For maximization of the pseudo-evidencepð^ hu jdÞ with respect to d, there is an information-theoretical interpretation of the trade-off between data fitting and model complexity [21,36] that can be demonstrated as follows, where we use Bayes’ Theorem and the hierarchical model structure in Fig. 1:
ð34Þ
Eq. (34) shows that the log pseudo-evidence, which is to be maximized, is the difference between the posterior mean of the log likelihood function for h (the first term) and the relative entropy (or Kullback-Leibler information) of the posterior with respect to the prior distribution for h (the second term). The first term quantifies the ability of the hierarchical model to match the MAP vector ^ hu which is obtained from the calibration state; it is
hu;j in the posterior maximized if all hj are tightly clustered around ^ for h (i.e. all aj ? 0). The second term reflects the amount of information extracted from the pseudo-data ^ hu ; it is minimized if h is
independent of ^ hu ; conditional on (x2, /, b) (i.e. all aj ? 1). This term penalizes models that have fewer parameter components hj ^u ; and therefore forces the model updating differing from those in h
^u and relatively more from the to extract less information fromh system modal parameters x2 and /, and so from the ‘‘measured” ^ as implied by the observations of the hierar^ 2 and w, modal data x chical model structure exhibited in Fig. 1. However, overextraction of information from the system modal parameters x2 and / will produce a structural model vector h with too large of a difference from ^ hu that is overly sensitive to the details of the
information in the specified system modal parameters x2 and /, ^ In ^ 2 and w: and therefore in the ‘‘measured” modal parameters x other words, the sensor noise and other environmental effects may not be ‘‘smoothed out” so they may have a large effect on the damage detection performance. To summarise, (34) demonstrates that the evidence maximization procedure over hyper-parameter vector d automatically involves a trade-off between the average data-fit of the model class defined by d and the information it extracts from the data. This is the Principle of Model Parsimony or the Bayesian Ockham razor [21,37] at work. The Bayesian procedure is effectively implementing Ockham’s Razor by assigning lower probabilities to a structural model whose parameter vector h has too large or too small differ^u obtained at the calibration state (too few or too ences from h
many aj ? 0), therefore suppressing the occurrence of false and missed damage detection alarms.
4. Proposed algorithms for identification of stiffness loss The above formulation leads to the following algorithms for probabilistic inference of structural stiffness loss. There are two variants of the algorithm. For the calibration stage, model sparseness is not expected and hence Algorithm 1 is used without optimization of the hyper-parameter vector a; instead, we fix all components aj with some large values. For the monitoring stage, Algorithm 2 is used where all hyper-parameters in d are estimated to ensure sparse stiffness loss inference.
43
Y. Huang et al. / Structural Safety 64 (2017) 37–53
4.1. Algorithm 1: Bayesian inference for the structural model in the calibration stage
structural stiffness reductions in the monitoring stage, which is summarized in Algorithm 2.
Algorithm 1
Algorithm 2
INPUTS:
INPUTS:
^ ¼w ^ u for Nm ^2 ¼ x ^ 2u and w . Identified Ns sets of modal data x modes from the calibration stage (N s P 3) . A chosen nominal value h0 of hu (e.g. from a finite-element model) and associated prior covariance matrix R0 ¼ r20 INh with large variance
r20 .
^ 2; ~ ¼ 10ðN s N m N o 2Þ=kwk . Initial parameter values for g PN s 4 ^ r;i ; i ¼ 1; . . . ; N m , q~ i ¼ 10ðNs 2Þ= r¼1 x PN s ~ ~ , where w ~ ¼ 20b ^ r;1 k2 and b ^ 21 MÞ b0 ¼ 2N s N o = r kðKðh0 Þ x 0 K and M are the reduced stiffness and mass matrices that correspond to the measured DOFs ALGORITHM: P s 2 ^ r =N s ~ 2 ¼ Nr¼1 1. Initialize l = h0, R = R0 and x x ~ ~ from input values ~ ; b0 and b ~; q 2. Initialize g ~ with large values (e.g. a ~ j ¼ 109 ) 3. Fix all components in a 4. While convergence criterion on ^ hu is not met (Outer loop) ~ is not met (Inner 5. While convergence criterion on b loop) ~ using (26), and then update g ~ using (23) 6. Update MAP / ~ 2 using (27), and then update q ~ using 7. Update MAP x (28) ~ using (32) 8. Update MAP b 9. Calculate the conditional posterior mean ^ hu ¼ l and covariance matrix Ru = R for h = hu using (17a) and (17b) ~ between the current 10. End while (the ratio of change of b iteration and the previous iteration of the inner loop is sufficiently small (e.g. smaller than 0.001)) ~0 using (33) 11. Updateb 12. End while (the ratio of change in all components of ^ hu between the current iteration and the previous iteration of the outer loop is sufficiently small (e.g. smaller than 0. 01)) OUTPUTS: . MAP estimates of all uncertain parameters in d . Posterior mean ^ hu and covariance matrix Ru of the stiffness scaling parameters hu at the calibration stage
^ ¼w ^ d for Nm ^2 ¼ x ^ 2d and w . Identified Ns sets of modal data x modes from the monitoring stage (N s P 3) . MAP estimate ^ hu and the associated posterior covariance matrix Ru from the calibration stage . Initial parameters values for ~ ¼ 2N s N o =PNs kðKð ~ , and ^ w ~ ¼ 20b ^ r;1 k2 and b ^ 21 MÞ b hu Þ x 0 0 r ~ i are the same as in Algorithm 1 ~ and q initial values for g ~ u of the system natural frequencies ~ 2 and / . MAP estimates x u
and mode shapes from the calibration stage ALGORITHM: ~ j ¼ N 2h , j = 1, . . ., Nh 1. Initialize a ~ as their MAP estimates from the cal~ 2 and / 2. Initialize x ibration stage ~0 and b ~ from input values ~; b ~; q 3. Initialize g 4. While convergence criterion on ~ hd is not met (Outer loop) ~ is not met (inner 5. While convergence criterion on a loop) ~ j ¼ 0 if aj (j = 1, . . ., Nh) becomes smaller than amin 6. Set a (chosen as 109 in the example later) ~ j ¼ 0 if ~ 7. Set a hj > ^ hu;j ðj ¼ 1; . . . ; N h Þ when the iteration number l > 2 8. Calculate the conditional posterior MAP value ~ hd ¼ l and covariance matrix Rd = R using (17a) and (17b) 9. Do steps 6–8 in Algorithm 1 ~ using (31) 10. Update a ~ 11. End while (the ratio of change in all components of a between the current iteration and the previous iteration of the inner loop is sufficiently small (e.g. smaller than 0.001)) ~0 using (33) 12. Update b 13. End while (the ratio of change in all components of ~ hd between the current iteration and the previous iteration of the outer loop is sufficiently small (e.g. smaller than 0.01)) OUTPUTS: . MAP estimates of all uncertain parameters in d . Posterior mean ~ hd and covariance matrix Rd of the stiffness scaling parameters hd at the monitoring stage
4.2. Algorithm 2: Sparse Bayesian learning for inference of structural stiffness reductions in the monitoring stage When real data from a structure is used, the modeling uncertainty in the structural stiffness may be too large to have the stiffness reduction reliably detected by most existing methods. To regularize a possibly severe ill-posed inverse problem, Ching and Beck [12,14] imposed a constraint in their Expectation–Maximiza tion algorithm that the values of the stiffness parameters for a damaged state cannot exceed those for the calibration state. However, we have not found an analytically tractable solution for SBL when imposing this additional prior information in the formulation. Therefore, we present an alternative strategy: a0j s for those substructures that have their MAP estimate ~ hj > ^ hu;j for the current iteration will be set to zero directly (the first two iterations will be discarded to make sure that the parameters in d are not directly influenced by the initial setting). Then the prior information that any stiffness increase from the calibration model is unacceptable is effectively utilized in our SBL procedure for inference of
Remark 4.1. Algorithm 2 is performed by iterating between the two groups of parameters, hd and d. It starts by considering all substructures as possibly damaged (the reason for the choice a~ j ¼ N2h ; j = 1, . . ., Nh, at the first iteration) and then causes the ~ j < 109 ; to be exactly ‘‘inactive” components hj of hd, which have a ^u;j from the calibration stage when optimizing over the equal to h hyper-parameters aj. After convergence, there are only a few ‘‘active”h0j s that are changed from their calibration values ^ hu;j and their corresponding substructures are considered to be damaged.
Remark 4.2. After obtaining the posterior PDFs of the stiffness scaling parameters hd and hu for both the current monitoring stage (where damage is possible) and the undamaged calibration stage, respectively, the posterior mean and standard deviations for hd and hu can be used to construct an approximation for the integrals involved (see (54) of [39]), allowing computation of the probability
44
Y. Huang et al. / Structural Safety 64 (2017) 37–53
that each substructure stiffness parameter has been reduced by more than a specified fraction of its value in the initial calibration stage [2,3,12,15,39].
explicitly addressed in the new SBL algorithm. We discuss these issues in this section. 5.1. Theory
Remark 4.3. The chosen Gaussian likelihood function for the stiffness scaling parameter h in (9) allows the stiffness to increase and ^u with equal chance. For many decrease from the calibration state h structures of interest, a prior judgment can be made that the structure is expected to lose stiffness where it is damaged. For example, structural damage in steel members tends to be at connections with local buckling or weld fracture, while in R/C members, damage usually occurs because of spalling of concrete and breaking of the bond between the reinforcing steel and the surrounding concrete. In these cases, laboratory tests of structural members have shown that there is a reduction in stiffness. For such cases, we impose the constraint in Step 7 of Algorithm 2 to alleviate the ill-conditioning in real damage inversions, that is, the solution for the MAP parameter~ d in Algorithm 2 is obtained by solving the following constrained optimization:
^ ^hu Þ; subject to ^ 2 ; w; max pðdjx
lj 6 ^hu;j ðj ¼ 1; . . . ; Nh Þ
where the posterior mean l of hd is a function of d. Even though ~ j ¼ 0 for any ~ it is a ‘‘greedy” optimization if we set a hj > ^ hu;j ðj ¼ 1; . . . ; N h Þ during any iteration (except the first two, once an aj = 0 in a certain iteration, it stays at zero), the strategy suppresses those local maxima that correspond to models having some substructures with stiffness increases from the calibration state; this significantly increases the probability of finding the global maximum with respect to d. In addition, under the condition that the modeling errors are zero-mean, there is a possible compensation between the inferred stiffness increases and decreases for different substructures that result from modeling errors alone. Therefore, by suppressing stiffness increases, the strategy helps to reduce the occurrence of inaccurate stiffness reductions and also induces a higher degree of sparseness in the identified stiffness reduction results. Remark 4.4. In our new formulation, the Gaussian and Exponential PDFs have a scale invariant property and their associated hyper-parameters are also properly learned. This leads to the desirable property that the algorithms are automatically scaleinvariant: the identification of the stiffness scaling parameters is independent of any scaling of the measured mode shapes and of the units chosen for the mass and stiffness matrices. Our algorithm in [39] does not have this scale-invariance property and so care must be taken in the selection of the scale parameter b0 for the equation-error precision b.
5. Comparison of new SBL algorithm with that in [39] In [39], the hyper-parameters a for the likelihood function in (9) are learned by maximizing the evidence function R ^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; pðx hu jaÞ ¼ pðx hu jh; n; aÞpðh; njaÞdhdn along with the stiffness scaling parameter h and other parameters n = {(x2)T, qT, sT, /T, g, b} to produce sparseness in the inferred stiffness reductions. For tractable estimation of the evidence function, two ^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; approximations pðx hu jaÞ pðx hu j~ n; aÞ and ~ n^ n are uti2 ~ ¼ argmax pðnjx ^ h ^ is obtained by maxi^u Þ and n ^ ; w; lized, where n ^ h ^u Þ. The excellent performance of the ^ 2 ; w; mizing the PDF pðn; hjx proposed SBL method in [39] has been verified with synthetic data. However, there are several important issues that needed further exploration for reliable application to real data and they have been
The
approximation
of
the
evidence
function
^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; pðx hu jaÞ pðx hu j~ n; aÞ in [39] is valid only if ^ h ^u Þ has a unique maximum at ~ ^ 2 ; w; n with a sharp peak at its pðnjx mode. Furthermore, for real structural damage detection where there are large modeling errors, it is challenging to confirm the accuracy of the approximation ~ n^ n employed in the formulation in [39]. Therefore, these approximations to estimate the evidence function may possibly be unreliable. The new formulation presented in Sections 3.1 and 3.2 establishes the probability of localized stiffness reductions without these possibly unreliable approximations and is therefore more rigorous. Another benefit of the new formulation is that all parameters in d are learned by an evidence maximization procedure and therefore the inferred sub-structural stiffness reductions are expected to be sparser than in [39]. Note that due to the hierarchical nature ~ 2 andq ~; x ~ i are estimated of the model in Fig. 1, the MAP values of g in the same way as those in [39]. However, the new terms ~ RU ~ diagðtrðRT ~ Nm N ÞÞ; ~ 1 Þ; . . . ; trðRU ~ Nm N Þ ~ 1 Þ; . . . ; trðRT b½trð and b d d PNh 1 ~ ~ ~ ~ R in (22) and (32) for the MAP values / and b are 1 a jj j¼1 j important for the algorithm to produce more sparse substructural stiffness reductions that help to suppress the possible occurrences of false damage detections. In addition, the information-theoretic explanation for the sparseness in the inferred stiffness reductions in Section 3.3 is important to interpret how the proposed algorithm suppresses the occurrences of false and missed damage detections, which was not discussed in [39]. 5.2. Computation for real applications Although the formulation in [39] does not involve solving an eigenvalue problem during each iteration, the dimension (NmNd) of the system mode shape vector / is a concern for computational efficiency since its MAP estimation has a computational effort of OðN 3m N 3d Þ, which is computationally demanding when the dimension NmNd is not small. In the new algorithm, we enhance the computational efficiency by using (26) to reduce the computational cost by a factor of N 2m . The second computational improvement is that our new algorithm, unlike our algorithm in [39], has the desirable property of scale-invariance, i.e., the identified stiffness reductions are independent of the choice of units for both the stiffness and mass matrices, as well as the scaling of the measured mode shapes. This property is a consequence of choosing the value of b0 adaptively to scale the prior distribution of the equation-error precision b. In [39], we suggested an a priori selection of b0 based on an initial empirical study but, unfortunately, this can lead to unreliable results for both the undamaged and possibly damage states because they are somewhat sensitive to b0. In the new algorithm, b0 is estimated by its MAP value, which is given by (33), and it is updated separately in an outer loop initialized from a reasonable value as specified in (A18). The numerical and experimental results given later support the effectiveness of this choice of b0. For application to real structural response data, the modeling error sometimes dominates over the useful data information for constraining the parameter updating and then it becomes challenging to implement the damage detection algorithm in [39]. For example, an increase in the identified substructure stiffness can sometimes occur. In the new algorithm, a more effective regu-
Y. Huang et al. / Structural Safety 64 (2017) 37–53
larization strategy is proposed based on our prior expectations in order to further alleviate the ill-conditioned stiffness inversion by suppressing mean stiffness increases during the optimization to find the MAP estimates of the parameters. 6. Illustrative applications The applicability of the proposed methodology to identify substructure stiffness reductions is illustrated with the IASC-ASCE Phase II SHM benchmark problems, first using the simulated data [41], and second using the experimental data [42]. The first application serves to verify the proposed methodology while the second one illustrates its application to real data from a damaged structure. The benchmark structure is a four-story, two-bay by twobay steel braced-frame structure and a diagram for the analytical model is depicted in Fig. 2 along with its dimensions, in which the x-direction is the strong direction of the columns. A detailed description of the benchmark problem, including detailed nominal properties of the structural elements in the analytical model, can be found in [12]. 6.1. Simulated Phase II benchmark problem Simulated Phase II benchmark problem is defined in Bernal et al. [41] and consists of both brace damage and unbraced beam-column joint damage cases (see [12] for the various damage pattern cases considered here). For instrumentation scenarios, only the results of the challenging partial-sensor scenario are presented, where only the measurements at the third floor (mid-height) and the roof are available. For generation of simulated time-domain data and the extraction of the modal properties, the reader can refer to [12] or [14] for detailed information. To locate the faces sustaining brace damage, the identification model class is defined with the stiffness matrix K parameterized as
KðhÞ ¼ K0 þ
XX s
v
hsv Ksv
ð35Þ
where s = 1, . . ., 4 refers to the story number and v = +x, x, +y, y indicates the direction of the outward normal of each face at each floor, yielding a stiffness scaling parameter vector h with sixteen components, corresponding to sixteen substructures (four faces of four stories). The K0sv s are the ‘nominal’ stiffness matrices to make the nominal value of each stiffness scaling parameter to be 1.0.
Fig. 2. The diagram of the benchmark structure [12,14].
45
For the braced structure, they are computed based on shearbuilding assumptions by projecting the full stiffness matrix for the 120-DOF benchmark structural model [41] onto 12 DOFs for a 3-D shear-building model with rigid floors (three DOFs per floor: translations parallel to the x- and y-axes and rotation about the zaxis). K0 is the nominal stiffness matrix contribution from the columns and beams. For the unbraced structure, we project the 120-DOF benchmark structural model [41] onto 20-DOFs for a 3-D shear-building model that assumes rigid floors in the x–y plane and allows rotation along the x- and y-axes (five DOFs per floor: translations parallel to the xand y-axes and rotation about the x-, y- and z-axes). It is assumed that the rotational stiffnesses of all beam-column connections in the same story along the x-axis (or along the y-axis) are identical. Two parameters are used for the rotational stiffness in each story to give eight stiffness parameters in total:
KðhÞ ¼ K0 þ
XX hsu Ksu s
ð36Þ
u
where s = 1, . . ., 4 refers to the story number and u = x, y indicates the axis along which the rotational stiffness is active. The Ksu are the nominal rotational stiffness matrices computed based on the model assumptions for the original undamaged structure; K0 is the nominal stiffness matrix contribution from the columns and beams. Ideally, we would like to treat each structural member (e.g., brace, beam or column) as a substructure so that we can infer from the dynamic data which, if any, members have been damaged by the event. However, the information available from the structure’s local network of sensors will generally be insufficient to support a member-level resolution of stiffness loss from damage, so larger substructures consisting of assemblages of structural members may be necessary in order to reduce the number of model parameters in h. A tradeoff is therefore required between the number of substructures (and hence the resolution of the damage locations) and the reliability of the probabilistically-inferred damage state. In this simulated data example, we only investigate the partialsensor scenario. During the calibration stage, we utilize Algorithm 1 based on 100 independent sets of identified modal parameters from the undamaged structure to find the MAP value of the structural stiffness scaling parameters^ hu and the corresponding c.o.v. (coefficients of variation), which are shown in Tables 1 and 2, for the undamaged braced structure (column RB.ps) and unbraced structure (column RU.ps), respectively. It is evident from the tables that the MAP estimates of all stiffness scaling parameters hj are accurate, that is, the errors are smaller than 3.7% and 2.5% for the braced and unbraced structure, respectively. ^u from During the monitoring stage, we choose the MAP value h the calibration stage as pseudo-data for h and perform Algorithm 2 based on ten sets of identified modal parameters as the primary data. The identified stiffness ratios of the MAP stiffness parameters with respect to those for the undamaged case (RB) are tabulated in Tables 1 and 2, for the four brace and three joint damage patterns, respectively, along with their corresponding c.o.v. The actual damaged locations are made bold for comparison. It is found that the damage patterns are reliably detected in both qualitative and quantitative ways. The identified stiffness ratios for the actual damaged substructures are close to their actual values (0.943 and 0.887 for the brace damage cases and 0.500, 0.667, 0.750 and 0.833 for the joint damage cases, see [12]) and all components with non-bold font have stiffness ratios that are exactly equal to one. In addition, the corresponding c.o.v. for the exactly identified non-damaged components is equal to zero, which means these substructures have no stiffness reduction with full confidence compared with that of the calibration stage. This is a benefit of
46
Y. Huang et al. / Structural Safety 64 (2017) 37–53
Table 1 Brace damage cases in Simulated Phase II benchmark: Identification results from Algorithms 1 and 2 for the partial-sensor scenario. Parameter
h1,+x h2,+x h3,+x h4,+x h1,+y h2,+y h3,+y h4,+y h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
RB.ps
DP1B.ps
DP2B.ps
DP3B.ps
MAP value ð^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
DP3Bu.ps c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
0.984 1.010 1.026 0.963 0.981 0.992 1.005 1.012 0.997 0.985 0.976 1.021 0.998 0.995 0.988 1.006
0.192 0.094 0.094 0.075 0.235 0.122 0.100 0.098 0.325 0.162 0.171 0.125 0.221 0.117 0.093 0.097
1.000 1.000 1.000 1.000 0.865 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.896 1.000 1.000 1.000
0.000 0.000 0.000 0.000 0.203 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.196 0.000 0.000 0.000
1.000 1.000 1.000 1.000 0.941 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.936 1.000 1.000 1.000
0.000 0.000 0.000 0.000 0.183 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.179 0.000 0.000 0.000
1.000 1.000 1.000 1.000 0.858 1.000 0.952 1.000 1.000 1.000 1.000 1.000 0.886 1.000 0.954 1.000
0.000 0.000 0.000 0.000 0.208 0.000 0.087 0.000 0.000 0.000 0.000 0.000 0.201 0.000 0.081 0.000
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.867 1.000 0.957 1.000
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.183 0.000 0.069 0.000
⁄.ps denotes the partial-sensor scenario. ⁄The true ratio value for h1,+y for DP1B and DP3B and h1,y for DP1B, DP3B and DP3Bu are 88.7% of the undamaged values. ⁄The true ratio value for h1,+y for DP2B, h3,+y for DP3B, h1,y for DP2B and h3,y for DP3B and DP3Bu are 94.3% of the undamaged values.
Table 2 Joint damage cases in Simulated Phase II benchmark: Identification results from Algorithms 1 and 2 for the partial-sensor scenario. Parameter
RU.ps
h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
DP1U.ps
DP2U.ps
DP1Uu.ps
MAP value ð^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
0.984 1.010 1.026 0.963 0.981 0.992 1.005 1.012
0.192 0.094 0.094 0.075 0.235 0.122 0.100 0.098
1.000 1.000 1.000 1.000 0.528 0.615 1.000 1.000
0.000 0.000 0.000 0.000 0.902 1.346 0.000 0.000
1.000 1.000 1.000 1.000 0.737 1.000 1.000 1.000
0.000 0.000 0.000 0.000 0.633 0.000 0.000 0.000
1.000 1.000 1.000 1.000 0.797 0.799 1.000 1.000
0.000 0.000 0.000 0.000 0.595 1.085 0.000 0.000
⁄The true ratio values for h1,y and h2,y for DP1U.ps are 50% and 66.7% of the undamaged values, respectively. ⁄The true ratio values for h1,y for DP2U.ps are 66.7% of the undamaged values. ⁄The true ratio values for h1,y and h2,y for DP1U.ps are 75% and 83.3% of the undamaged values, respectively.
the proposed sparse Bayesian formulation, which reduces the uncertainty of the unchanged components. It produces sparse ~ sv ! 0 models by learning the hyper-parameter a, where a ~ su ! 0Þ implies that R(sv)(sv) ? 0 (or R(su)(su) ? 0) and ðor a ~ hu Þ ðor ~ hsu ! ð^ hu Þ Þ. hsv ! ð^ sv
su
In order to further portray the identification performance of the proposed algorithms, the identified stiffness reduction ratios of Algorithm 2 are compared with those of our previous Bayesian inversion algorithm in [39] for various damage patterns, and the results are presented in Fig. 3, where the stiffness reduction ratio is computed as the difference between each MAP value ^u Þ ðor ðh ^u Þ Þ of the stiffness scaling parameters from the caliðh sv
su
~su ) from the current monbration stage and the MAP value ~ hsv (or h ^u Þ (or ð^ h itoring stage, normalized by ðh u Þsu ). A variant of Algorithm sv 2 with the stiffness increase constraint removed (Step 7 is omitted) is also studied in order to investigate the effect of this constraint on the identification performance. The results for the modified Algorithm 2 without constraints and for Algorithm 2 are presented in Fig. 3(a),(b) in the columns (iii) and (iv), respectively, while for the method in [39], the results with the Laplace-prior parameter k = 0 and k optimized are given in Fig. 3(a),(b) in the columns (i) and (ii), respectively. Note that a uniform hyper-prior over a is utilized for Algorithm 2, which is equal to the Laplace prior case with k ¼ 0: In Fig. 3, the 16 (or 8) components hsv (or hsu) of h are presented in the order h1, . . ., h16 (or h1, . . ., h8) that corresponds to the order listed in the first column of Table 1 (or Table 2).
The results for Algorithm 2 without constraints in Fig. 3(b)(iii) are more sparse in comparison with the results of the method in [39] with k ¼ 0 (Fig. 3(b)(i)), which shows the advantage of the new formulation in Section 3 as discussed in Sections 5.1 and 5.2. The results for Algorithm 2 are even better with no false damage detections observed, and hence an appropriate level of sparseness of stiffness reductions is inferred by imposing the constraint in Step 7 to suppress stiffness increases. The algorithm in [39] with the optimal Laplace prior parameterk has the same perfect damage assessment performance as Algorithm 2 for the brace damage patterns (compare Fig. 3(a)(ii),(iv)) but, as shown in Fig. 3(b)(ii),(iv), this is no longer true for the joint damage patterns; in fact, Fig. 3 (b)(ii) shows that totally over-sparse stiffness reduction models with no nonzero components are produced for the algorithm with optimalk in the DP2U.ps and DP1Uu.ps cases. This is because when finding the MAP value of a by maximiz^ ^ ^ ^ ^ 2 ; w; ^ 2 ; w; ingpðajx hu ; kÞ / pðx hu jkÞpðajkÞ, the Laplace hyperpriorpðajkÞ in the method of [39] will dominate over the evi^ h ^u jaÞ in these cases. ^ 2 ; w; dencepðx A more complete picture of substructure damage is given by the dam probabilities of damage, Pdam sv ðf Þ (or P su ðf Þ), with different fractional stiffness losses f e [0, 1], calculated using the theory presented in [39]. In Fig. 4, the probability of damage is shown for the sixteen h0sv s (or eight h0su s) for various damage patterns from applying Algorithm 2. It is clear in all cases that the substructures with real stiffness losses have damage with a probability of unity.
47
Y. Huang et al. / Structural Safety 64 (2017) 37–53
DP1B.ps
0.2
0.1
0.1
0.1
0
0
0
0.2 0.1
DP2B.ps 0.2 DP2B.ps 0.1
0 0.2
0.2
DP3B.ps
0 4 8 12 16
DP2B.ps
0.2
DP3B.ps
DP3B.ps
0.1
0.1
0.1
0.1
0
0
0
0
0 4 8 12 16
0 4 8 12 16
0 4 8 12 16
DP3B.ps
0 4 8 12 16
0.2
0.2 DP3Bu.ps 0.2 0.2 DP3Bu.ps DP3Bu.ps DP3Bu.ps 0.1 0.1 0.1 0.1 0
0 0 4 8 12 16 (i)
0
0 0 4 8 12 16 (ii)
0 4 8 12 16 (iii)
0 4 8 12 16 (iv)
(a)
DP 1U.ps
0.8
DP1U.ps
DP 1U.ps
0
0 0 2 4 6 8
0.8 0.4
0.4
0 0 2 4 6 8
0.6 0.4 0.2 0
0.8 0.4
0
0 4 8 12 16 0.2
DP1U.ps
0.4
0 0 4 8 12 16
0.2
0.8
DP2B.ps
0.1
0 0 4 8 12 16
0.2
0.2 DP1B.ps 0.1 0
0.1
0 0 4 8 12 16
DP1B.ps
0 4 8 12 16
0 4 8 12 16
0 4 8 12 16 Stiffness reduction ratios
0.2
DP1B.ps
Stiffness reduction ratios
0.2
0 2 4 6 8
0 2 4 6 8
0.6 0.6 0.6 DP2U.ps DP2U.ps 0.4 DP2U.ps 0.4 DP2U.ps 0.4 0.2 0.2 0.2 0 0 0 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8
0.6 0.6 0.6 0.6 DP1Uu.ps DP1Uu.ps DP1Uu.ps DP1Uu.ps 0.4 0.4 0.4 0.4 0.2 0.2 0.2 0.2 0 0 0 0 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8
(i)
(iii)
(ii)
(iv)
(b)
Fig. 3. Simulated Phase II benchmark for (a) Brace damage cases and (b) joint damage cases: Comparison of stiffness reduction ratios from the two Bayesian inversion methods for the substructures for the investigated damage patterns using the data of the partial-sensor scenario. Column (i) Algorithm in [39] with k ¼ 0; (ii) Algorithm in [39] with k optimized in each iteration; (iii) Algorithm 2 with no constraints; (iv) Algorithm 2. (The bars colored red are ratios that agree with the actual damage locations. Also, the negative bars in (a), (b) (i) and (iii) correspond to small stiffness gains). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
For the substructures that are undamaged, most of the median fractional stiffness losses are close to zero, and the probability is very small that the damage extent exceeds 2%, which is within the tolerable uncertainty level for applications.
6.2. Experimental Phase II benchmark problem In the second example, the proposed methodology is applied to both the brace and joint damage patterns in the IASC-ASCE experimental Phase II benchmark problem [42]. Nine configurations are considered, in which Configs. 1–6 are five brace damage patterns with Config. 1 being the calibration (undamaged) case, and Configs. 7–9 are two joint damage patterns for the unbraced structure with Config. 7 being the calibration case. The five brace damage patterns considered are as follows: 1) Config. 2: removal of all braces (100% stiffness loss) on the –y face; 2) Config. 3: removal of the left-handside brace in each story on the –y face; 3) Config. 4: removal of the left-hand-side braces in the first and fourth stories on the –y face; 4) Config. 5: removal of the left hand-side brace in the first story on the –y face; 5) Config. 6: removal of two braces in the second story on the + x face. For the two joint damage patterns, we have: 1) Config. 8: loosen both ends of the right-hand-side beam at each floor on the y face in Fig. 2; 2) Config. 9: loosen both ends of the right hand-side beam at the first and second floors on the y face. For each configuration, experimental acceleration data from the benchmark structure measured after impact of a sledgehammer are utilized. Sensor data were collected from five Kinemetrics EPI sensors placed near the base and floor centers (Nodes 5, 14, 23, 32, and 41 in Fig. 2) that sensed the accelerations in the +y direction, and from ten force–balance accelerometers mounted on the +y and y faces of the base and all floors (Nodes 2, 8, 11, 17, 20, 26, 29, 35, 38, and 44 in Fig. 2) that sensed accelerations in the +x direction. For each configuration, the time-domain data from each sensor is divided into three segments and three independent sets of
experimental modal parameters are extracted by the modal identification procedure MODE-ID [43,44]. For each time segment, five modes (Nm = 5), consisting of the first and second translation modes in the x and y directions and the first torsion mode are identified for the brace damage cases, and eight modes (Nm = 8), consisting of the first, second and third translation modes in the x and y directions and the first and second torsion modes are identified for the joint damage cases. These modal identification results can be found in [12]. The same 3D 12-DOF and 20-DOF structural models adopted in the simulated data study are used for the brace and joint damage cases, respectively. We first run Algorithm 1 to get the MAP estimates of the stiffness scaling parameter ^ hu in the calibration stage. The available information of three independent sets of identified ^ u is insufficient to support an accurate ^ 2u and w modal parameter x inference of the stiffness scaling parameters for both cases. We therefore assume that the stiffness scaling parameters of +x and x faces for all floors, and +y and y faces for all floors, are identical in order to reduce the number of stiffness model parameters to two for the braced cases, while for the unbraced cases, we assume all stiffness scaling parameters are identical. During the monitor^u from the calibration stage is used as ing stage, the MAP value h pseudo-data for h and Algorithm 2 is implemented based on three ^ d . As in Section 6.1, ^ 2 and w sets of identified modal parameters x d
for brace and joint damage cases, all 16 and 8 components hj of h, respectively, are updated according to the order listed in the first columns of Tables 1 and 2. Based on the results from Algorithm 2, the stiffness ratios and the associated c.o.v. are tabulated in Tables 3 and 4 for the brace and joint damage cases, respectively. For the brace damage shown in Table 3, it is seen that all damage is effectively detected although there is a false detection for h16 = h4,y of Config. 5 with a 3.6% stiffness reduction; however, in this case, the difference in value for no damage (1.000) is less than one standard deviation (5.2%) from the MAP value of 0.964 for h4,y and so the probability of no damage is
48
Y. Huang et al. / Structural Safety 64 (2017) 37–53
1
0.8
1,+y
0.6
1,-y
Probability of damage
Probability of damage
1
0.4 0.2 0 -0.1
-0.05
0
0.05 0.1 Damage extent, f
0.15
0.2
0.8
1,+y
0.6
1,-y
0.4 0.2 0 -0.1
0.25
-0.05
0
0.05 0.1 Damage extent, f
(a)
0.6
1,-y
0.4
3,+y
0.2
3,-y
-0.05
0
0.05 0.1 Damage extent, f
0.15
0.2
Probability of damage
Probability of damage
1,+y
0.8 1,-y
0.6 0.4 0.2
3,-y
0 -0.1
0.25
-0.05
0
0.05 0.1 Damage extent, f
(c)
0.15
0.2
0.25
(d)
1
1 1,y
0.8
Probability of damage
Probability of damage
0.25
1
0.8
0.6 0.4 0.2 0 -0.4
0.2
(b)
1
0 -0.1
0.15
-0.2
0
0.2 0.4 0.6 Dama ge extent, f
0.8
1
0.8
1,y
0.6
2,y
0.4 0.2 0 -0.4
-0.2
0
0.2 0.4 0.6 Dama ge extent, f
(e)
0.8
1
(f) Probability of damage
1 0.8
1,y
0.6
2,y
0.4 0.2 0 -0.4
-0.2
0
0.2 0.4 Damage extent, f
0.6
0.8
1
(g) Fig. 4. Simulated Phase II benchmark: Probability of substructure damage exceeding f calculated using Algorithm 2 for the partial-sensor scenario: (a) DP1B.ps; (b) DP2B.ps; (c) DP3B.ps; (d) DP3Bu.ps;(e) DP1U.ps; (f) DP2U.ps and (g) DP1Uu.ps.
quite high. The MAP estimate of stiffness ratios for a particular face with removal of one and two braces are underestimated as approximately 50% and 30%, respectively, compared to the expected stiffness ratios of 77.4% and 54.9%, respectively [12]. For the quantification of the posterior uncertainty, the c.o.v. is generally much higher compared with those values in Table 1 for the simulated data benchmark due to larger modeling errors in this real
case, especially for those components corresponding to real damage locations. For the joint damage cases shown in Table 4, although most damage except for h3,y is detected, four false detections are found, and some of them indicate significant stiffness reduction (e.g., h4,x and h4,y of Config. 9). It is believed that the relative poor performance compared with the brace damage cases is due to the fact
49
Y. Huang et al. / Structural Safety 64 (2017) 37–53 Table 3 Brace damage cases in Experimental Phase II Benchmark: Identification results from Algorithms 1 and 2. Parameter
h1,+x h2,+x h3,+x h4,+x h1,+y h2,+y h3,+y h4,+y h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
Config. 1
Config. 2
Config. 3
Config. 4
Config. 5
Config. 6
MAP value ^u Þ ðh
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
MAP ratio ~ h ^u Þ ðh=
c.o.v. (%)
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658 4.509 9.658
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.302 0.370 0.268 0.160
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 36.148 58.243 32.306 45.287
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.443 0.646 0.386 0.293
0.000 0.000 0.000 0.000 1.593 0.000 0.000 0.000 0.000 0.000 0.000 0.000 28.930 34.911 25.493 33.986
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.767 1.000 1.000 0.608
0.000 0.000 0.000 0.000 0.000 1.981 0.000 0.000 0.000 0.000 0.000 0.000 22.028 0.000 1.359 22.014
1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 1.000 1.000 0.677 1.000 1.000 0.964
0.000 0.000 0.000 0.000 0.000 1.345 0.000 0.000 0.000 0.000 0.000 0.000 13.541 0.000 0.000 5.245
1.000 0.835 1.000 1.000 0.994 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
0.000 15.719 0.000 0.000 3.458 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
⁄The true ratio values for h1,y, h2,y, h3,y and h4,y for Config. 2 and h1,y for Config. 6 are 54.9% of the undamaged values. ⁄The true ratio values for h1,y, h2,y, h3,y and h4,y for Config. 3, h1,y and h4,y for Config. 4, h1,y for Config. 5 are 77.4% of the undamaged values.
Table 4 Joint damage cases in Experimental Phase II Benchmark: Identification results from Algorithms 1 and 2. Parameter
h1,x h2,x h3,x h4,x h1,y h2,y h3,y h4,y
Config. 7
Config. 8
Config. 9
MAP value ð^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
MAP ratio ð~ h=^ hu Þ
c.o.v. (%)
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
8.281 8.281 8.281 8.281 8.281 8.281 8.281 8.281
0.925 1.000 1.000 1.000 0.715 0.766 1.000 0.747
0.832 0.000 0.000 2.100 0.599 0.931 0.000 2.163
0.957 1.000 1.000 0.724 0.736 0.870 1.000 0.862
0.934 0.000 0.000 1.916 0.676 1.157 0.000 2.141
⁄The true ratio values for h1,y, h2,y, h3,y and h4,y for Config. 8 are 83.3% of the undamaged values. ⁄The true ratio values for h1,y and h2,y for Config. 9 are 83.3% of the undamaged values.
that the stiffness matrix of the structure is dominated by the columns and the installed sensors only measure the accelerations in the x and y directions, rendering the identification of rotational stiffness difficult. However, it is seen that the damage assessment performance here is better than in [14], since less false damage detections are produced and the identified stiffness ratios for the damaged substructures are closer to their actual values (0.833). In Fig. 5, the modified Algorithm 2 without constraints (Fig. 5 (c)) and Algorithm 2 (Fig. 5(d)) are compared with the method in [39] by presenting the stiffness reduction ratios for various damage cases, where hyper-parameterk ¼ 0 in Fig. 5(a) andk is optimized for the assigned Laplace prior in Fig. 5(b). It is interesting to see totally over-sparse stiffness reduction models with no nonzero components are produced when k is optimized. In comparison with the method in [39] with k ¼ 0; Algorithm 2 with no constraints has better performance (Fig. 5(c)), although false detections are still observed in h12 = h4,x for Config. 4, h16 = h4,y for Config. 6, and h4 = h4,x, h8 = h4,y for Config. 9 that correspond to significant identified stiffness reductions. In addition, some undamaged substructures show a significant increase in stiffness that is unrealistic. However, incorporating the additional constraint to suppress stiffness increases in Algorithm 2 results in exact localization of all stiffness reductions (Fig. 5(d)) for the brace damage patterns, although the modeling error seems to be too large to have the small effect of rotational stiffness damage reliably detected for
the joint damage patterns (especially Config. 9), which is an inherent difficulty of this particular damage identification case. Another reason for the better performance of Algorithm 2 compared with the method in [39] for the experimental data is that it automatically chooses the optimal scale parameter b0 whereas in [39] we select b0 by trial-and-error but the damage detection results can be sensitive to the value of b0. The damage probability curves for the stiffness scaling parameters h are plotted in Fig. 6 for Configs. 2–6 for brace damage patterns and Configs. 8–9 for joint damage patterns. The non-zero probabilities for negative fractional stiffness losses f are a result of the Gaussian approximations used in the calculation of the curves. For brace damage, all the actual damaged substructures are clearly shown to have a large damage probability when their fractional stiffness loss is close to the actual fraction. For the undamaged substructures, the posterior medians of the estimated damage corresponding to a damage probability of 0.5 are around zero, although in some cases small probabilities of about 0.2 are estimated for damage fractions exceeding 10%, which is due to the larger posterior uncertainty for the stiffness scaling parameters in the calibration and monitoring stages (Table 3). For the joint damage cases, h1,y, h2,y and h4,y for Config. 8 and h1,y, h2,y, h4,y and h4,x for Config. 9 are clearly shown to have a large damage probability, where false detection occurs in h4,y and h4,x for Config. 9 and missed detections occurs in h3,y for Config. 8 (see also Table 4).
50
Y. Huang et al. / Structural Safety 64 (2017) 37–53
1 0.5 0 -0.5 -1
Stiffness reduction ratios
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
Config2
0 4 8 12 16 Config3
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1
0 4 8 12 16 Config4
1 0.5 0 -0.5 -1
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
0
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
Config3
0 4 8 12 16 Config4
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
Config9
Config9
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
0 4 8 12 16 Config4
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 1 0.5 0 -0.5 -1 0.5
0 4 8 12 16 Config3
0 4 8 12 16 Config4
0 4 8 12 16 Config5
0 4 8 12 16 Config6
0 4 8 12 16 Config8
-0.5 0 2 4 6 8
0.5
Config2
0
Config9
0
-0.5
-0.5
1 0.5 0 -0.5 -1
Config3
-0.5
0
0
1 0.5 0 -0.5 -1
0 4 8 12 16 1 0.5 0 -0.5 -1
0 2 4 6 8 0.5
Config2
0
-0.5 0 2 4 6 8
1 0.5 0 -0.5 -1
0 4 8 12 16
0
-0.5 0.5
Config2
0 2 4 6 8 0.5
Config9
0 -0.5
-0.5
0 2 4 6 8
0 2 4 6 8
0 2 4 6 8
0 2 4 6 8
(a)
(b)
(c)
(d)
Fig. 5. Experimental Phase II benchmark: Comparison of stiffness reduction ratios from the two Bayesian inversion methods for the all substructures for five brace damage and two joint damage scenarios: (a) The algorithm in [39] with k ¼ 0; (b) The algorithm in [39] with k optimized in each iteration; (c) Algorithm 2 with no constraints; (d) Algorithm 2. (Note that a negative value corresponds to a stiffness gain).
7. Concluding remarks In the paper, we examine the problem of inferring damageinduced stiffness reductions in a structure from structural vibration response by taking a hierarchical Bayesian model perspective. For each chosen substructure, not only the most probable estimate of the spatially-sparse stiffness changes is inferred from the identified modal parameters but also the associated posterior uncertainties of the stiffness parameters are quantified, including the probability of substructure damage of various amounts. Our approach leads to an improved sparse Bayesian learning (SBL) algorithm for probabilistic damage detection compared to our previously proposed SBL algorithm in [39]. Apart from the three appealing features in the original method listed in [39], that there is no need to solve the nonlinear eigenvalue problem nor to match model and experimental modes explicitly, as well as the automated selection of the algorithmic parameters, the new method has the following four additional advantages. First, two possibly unreliable approximations in the original method in [39] have been eliminated. Second, the computation is more efficient in finding the MAP values of the system mode shapes. Third, the identified stiffness reductions are independent of the units of the stiffness and mass matrices, as well as the scaling of the identified
mode shapes. Fourth, a regularization strategy during hyperparameter estimation has been developed to further alleviate the ill-conditioning of the damage inversion for the case of large modeling errors when applied to real data. This strategy is based on the expectation that damaged structures actually lose stiffness at the damaged locations, so we incorporate an additional constraint during optimization for the MAP (maximum a posteriori) values of the parameters other than the stiffness scaling parameters, that the latter cannot have mean values in a damaged state that are greater than their mean values in the calibration state, which reduces the occurrences of incorrect identifications of stiffness changes. The improved effectiveness and robustness of the proposed method compared with our previously-published Bayesian method [39] was demonstrated through the analysis of the IASC–ASCE Phase II simulated and experimental benchmark problem to identify brace damage and damage to beam-column connections in the steel-frame benchmark structure. The results show that our new damage assessment method is reliable, despite significant error in the structural model in the real experiments. For future studies, it would be useful to explore the uncertainties in the identified system modal parameters and stiffness parameters from the calibration stage, since these are effectively ignored when only their MAP values are utilized.
51
Y. Huang et al. / Structural Safety 64 (2017) 37–53
1 4,-y
0.6
3,-y
0.4
2,-y
0.2
1,-y
0 -0.3
0
0.3 0.6 Damage extent, f
0.9
1.2
1
0.8
4,-y
0.6
3,-y
0.4
2,-y
0.2
1,-y
0 -0.3
0
0.3 0.6 Damage extent, f
(a)
1.2
0.6 0.4 1,-y
0.2 0 -0.3
0
0.3 0.6 Damage extent, f
1.2
(c)
1,-y Probability of damage
0.8 0.6 0.4 0.2 0
0.3 0.6 Damage extent, f
0.9
1.2
0.8 0.6 0.4 2,+x 0.2 0 -0.3
0
0.3 0.6 Damage extent, f
(d)
0.9
1.2
(e) 1
0.8
1,y
0.6
2,y
0.4
Probability of damage
1
4,y
0.2 0 -0.3
0.9
1
0 -0.3
Pro babili ty o f damage
4,-y
0.8
(b) 1
Pro babili ty o f damage
0.9
Probability of damage
0.8
P ro ba bili ty o f damage
Probability of damage
1
0
0.3 0.6 Damage extent, f
0.9
1.2
0.8
1,y
0.6
4,x
0.4
2,y 4,y
0.2 0 -0.3
0
0.3 0.6 Damage extent, f
(f)
0.9
1.2
(g)
Fig. 6. Experimental Phase II benchmark: Probability of substructure damage exceeding f using Algorithm 2 for the brace damage cases: (a) Config. 2; (b) Config. 3; (c) Config. 4; (d) Config. 5 and (e) Config. 6, and joint damage cases: (f) Config. 8 and (g) Config. 9.
Acknowledgements The first author is partially supported by the George W. Housner Earthquake Engineering Research Fund of the California Institute of Technology. The authors gratefully acknowledge this support. The first author is also partially supported by a grant from the National Natural Science Foundation of China (NSFC grant no. 51308161) and the International Postdoctoral Exchange Fellowship Program 2014 of the Office of China Postdoctoral Council (No. 20140018), and this support is also gratefully acknowledged.
Appendix A. MAP estimation of d Iterative equations for calculating the MAP estimate of d = [(x2)T, qT, sT, /T, g, m, aT, b, b0, j]T are obtained by differentiation of the logarithm of J1(d) defined in (21) with respect to d and then setting these derivatives to zero to find the stationary points of J(d) = log J1(d). Ignoring terms in the logarithm that are independent of d, the objective function is defined as:
T g ^ T ^ C/Þ 1 ðx ^ 2 Lx2 Þ E1 ðx ^ 2 L x2 Þ JðdÞ ¼ ðw C/Þ ðw 2 2 1 1 þ ðN d Nm Nh Þ log b log jHT Hj 2 2 1 T b T 1 T T ðb b b HðH HÞ H bÞ þ No Ns N m log g 2 2 Nm 1 1 X 1 þ Ns ðlog qi Þ log jA þ b1 ðHT HÞ j 2 i¼1 2 T 1 1 1 1 ðh^u ðHT HÞ HT bÞ ðA þ b1 ðHT HÞ Þ 2 Nm X 1 ð^hu ðHT HÞ HT bÞ þ log m mg þ ½log si si qi i¼1
1 b0 b
log b0 þ log j jb0
ðA1Þ
Following [39], we can simplify the following two terms in (A1): 1
log jA þ ðbHT HÞ j ¼ Nh logb log jHT Hj þ logjAj þ logjR1 j ðA2aÞ
52
Y. Huang et al. / Structural Safety 64 (2017) 37–53 1
T
1 1 1 ð^hu ðHT HÞ HT bÞ ðA þ ðbHT HÞ Þ ðh^u ðHT HÞ HT bÞ T
1
Then: NX NX m Nd m Nd @ðHT HÞ ¼ PTq H þ HT Pq ¼ PTq /s Ps þ /s PTs Pq @/q s¼1 s¼1 ! NX m Nd ¼ 2/q PTq Pq þ PTq /s Ps
1
¼ ð^hu ðHT HÞ HT bÞ ðbHT H bHT HðA1 þ bHT HÞ bHT HÞ 1
ð^hu ðHT HÞ HT bÞ T
1
1
¼ bðl ðHT HÞ HT bÞ HT Hðl ðHT HÞ HT bÞ T þð^hu lÞ A1 ð^hu lÞ
þ
1 1
where R = (bH H + A ) is the covariance matrix of h defined in (17b) and it is readily evaluated for arbitrary values of parameters /, a and b. Notice that (A2a) shows that the term with log |HTH| in (A1) cancels out. Now consider the sum of two of the terms in (A1): 1
T
@ðlogjR1 jÞ @ðHT HÞ ¼btr R @/q @/q
1
1
ðA3Þ
T
HT Hðl ðHT HÞ HT bÞ þ ð^hu lÞ A1 ð^hu lÞ T 2 ¼ bkHl bk þ ð^hu lÞ A1 ð^hu lÞ
Note that when differentiating (A3) with respect to uncertain parameters /, x2, b and a, the terms involving derivatives of l will drop out since:
bðdlÞ H ðHl bÞ þ bðHl bÞ Hdl þ ðdlÞ A ð^hu lÞ T
T
T
T
1
T
þð^hu lÞ A1 ðdlÞ ¼ ðdlÞ ðbHT H þ A1 Þl þ lT ðbHT H þ A1 Þdl bðdlÞ HT b T
T
ðA4Þ
bb Hdl dlT A1 ^hu ^hTu A1 dl T
T
T
¼ ðdlÞ R1 l þ lT R1 dl ðdlÞ R1 l lT R1 dl ¼ 0 Then for the term b||Hl b||2 in (A3), the derivatives with respect to /, x2 and b are expressed as:
@ðbkHl bk2 Þ=@/ ¼ @ðbkF/k2 Þ=@/ ¼ 2bFT F/ @ðbkHl bk2 Þ=@ x2 ¼ @ðbkGx2 ck2 Þ=@ x2 ¼ 2bðGT Gx2 GT cÞ
ðA5Þ where F is defined in (24) while G and c are defined in (29) and (30), respectively. The derivatives of log|R1| in (A2a) with respect to /, a and b are now derived. To differentiate log|R1| with respect to /, we first rewrite H as the following form:
6 H¼6 4
K1 /1 .. . K1 /Nm
3
KNh /1 7 .. .. 7 . . 5 KNh /Nm Nm N
¼
NX m Nd
/q Pq
ðA6Þ
q¼1 d N h
where /q are the NmNd components of / ¼ /T1 ; . . . ; /TNm and
@H Pq ¼ ; q ¼ 1; . . . ; Nm Nd @/q
as defined in (25).
!
ðA9Þ
Consequently, the derivative of log|R1| with respect to vector / is given by:
@ðlogjR1 jÞ ¼ 2b diag trðRT1 Þ; . . . ; trðRTNm Nd Þ / @/ T þ 2b trðRU1 Þ; . . . ; trðRUNm Nd Þ
ðA10Þ
Using the determinant derivative formula again, the derivatives of log|R1| with respect to aj and b are given by:
@ðlogjR1 jÞ @ðA1 Þ ¼ tr R @ aj @ aj
!
¼ a2 j Rjj
ðA11Þ
Nh X @ðlogjR1 jÞ ¼ trðRHT HÞ ¼ trðb1 RR1 b1 RA1 Þ ¼ b1 ð1 a1 j Rjj Þ @b j¼1
ðA12Þ where Rjj is the jth diagonal element of the posterior covariance matrix of R in (17b). Finally, setting the derivatives of J(d) in (A1) to zero leads to the update formulae for the parameters in d that are given in Section 3.2. We give some more details regarding (32) for the MAP estimation of equation error precision b, which is sensitive to the sparseness level of the identified stiffness reductions. Setting the ~ is estiderivative of J(d) with respect to b to zero, the MAP value b mated as:
@ðbkHl bk2 Þ=@b ¼ kHl bk2
2
¼ 2b/q trðRT q Þ þ 2btrðRU q Þ
¼ bðb b b HðHT HÞ HT bÞ þ bðl ðHT HÞ HT bÞT 1
Nm Nd s–q /s Ps
¼ 2b/q trðRTq Þ þ btrðRUq Þ þ btrðRUTq Þ
1
T
P
Therefore, using the Jacobi’s determinant derivative formula:
1
1 1 ðA þ b1 ðHT HÞ Þ ðh^u ðHT HÞ HT bÞ T
Pq ¼ 2/q Tq þ Uq þ UTq
where Tq ¼ PTq Pq and Uq ¼ PTq
bðb b b HðHT HÞ HT bÞ þ ð^hu ðHT HÞ HT bÞT T
! /s PTs
s–q
ðA2bÞ T
NX m Nd
ðA8Þ
s–q
T
2 RNm Nd 1
ðA7Þ
~¼ b
P h P h ~ ~ ~ 1 ~ 1 Nd N m Nj¼1 ð1 a N d Nm Nj¼1 ð1 a j Rjj Þ j Rjj Þ ¼ P 2 N 2 m ~ þ 2b ~1 ~1 ~l ~ i k þ 2b ~ 2 MÞ/ ~Þ x ~ bk kðKðl kH 0
i¼1
i
ðA13Þ
0
~0 and j ~ from the We continue by estimating the MAP values b derivatives of J(d) as follows:
~0 ¼ b
~ 2b qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ ~b 1 þ 1 þ 4j
j~ ¼ 1=b~0
ðA14Þ
ðA15Þ
~0 into (A 14), we get the MAP estimate ~ ¼ 1=b If we substitute j of b0 given in (33):
~0 ¼ b=2 ~ b
ðA16Þ
As a final point, we note that for real damage detection when the model has large modeling errors, it was found initially that during optimization for the MAP parameters, the components of
Y. Huang et al. / Structural Safety 64 (2017) 37–53
system mode shape / sometimes tended to zero along with the corresponding equation errors e, as mentioned in Remark 2.1. We found that to avoid this problem, we needed to carefully ~0 ~ and its associated scale parameterb choose the initial values of b in the iterative procedure in Algorithms 1 and 2. First, based on (A13) and (A16) but simplifying some terms in (A13):
~¼P b Ns r
4Ns No w ^ r;1 k2 ~ 21 MÞ kðK x
~0 ¼ P b Ns r
2Ns No w x ^ r;1 k2 ~ 21 MÞ kðK
ðA17Þ
ðA18Þ
^ r;1 2 RNo gives the identified mode shape of the 1st mode where w from the rth data segment. (A17) and (A18) are established by P m ~ i k2 =N m ~ 1 k2 ~ 2i MÞ/ ~ 21 MÞ/ replacing Ni¼1 kðK x withkðK x and P 2 N 2 m 2 ^ x ~ 21 MÞ/1 k =Nd with ~ kðK x kð K k =4N N . Here, K MÞ w r;1 s o i¼1 i are the reduced No No stiffness and mass matrices, respecand M tively, that correspond to the measured DOFs. Then, to achieve a good performance in the optimization for the MAP values in Algo~0 initially. Presumably, ~ ¼ 20b rithms 1 and 2 in Section 4, we take b
~ works because it is better if we make the larger initial value for b the model updating extract a larger amount of information from the system modal parameters x2 and / in the beginning, and there^ before pro^ 2 and w, fore from the ‘‘measured” modal parametersx ducing sparseness of the inferred stiffness reductions. References [1] Doebling SW, Farrar CR, Prime MB. A summary review of vibration-based damage identification methods. Shock Vib Digest 1998;30(2):91–105. [2] Vanik MW, Beck JL, Au SK. Bayesian probabilistic approach to structural health monitoring. J Eng Mech 2000;126(7):738–45. [3] Beck JL, Au SK, Vanik MW. Monitoring structural health using a probabilistic measure. Comput-Aid Civ Infrastruct Eng 2001;16(1):1–11. [4] Ou JP, Li H. Structural health monitoring in mainland China: review and future trends. Struct Health Monitoring 2010;9(3):219–31. [5] Farrar CR, Worden K. Structural health monitoring: A machine learning perspective. John Wiley & Sons; 2013. [6] Lynch JP. An overview of wireless structural health monitoring for civil structures. Phil Trans R Soc A 2007;365:345–72. [7] Li H, Huang Y, Ou J, Bao Y. Fractal dimension-based damage detection method for beams with a uniform cross-section. Comput-Aid Civ Infrastruct Eng 2011;26(3):190–206. [8] Moaveni B, Stavridis A, Lombaert G, Conte JP, Shing PB. Finite-element model updating for assessment of progressive damage in a 3-story infilled RC frame. J Struct Eng 2013;139:1665–74. [9] Bernal D. Load vectors for damage localization. J Eng Mech 2002;128(1):7–14. [10] Vanik MW. A Bayesian probabilistic approach to structural health monitoring, Technical Report EERL 97-07. Caltech, Pasadena, CA: Earthquake Engineering Research Laboratory; 1997. [11] Sohn H, Law KH. A Bayesian probabilistic approach for structure damage detection. Earthquake Eng Struct Dyn 1997;26(12):1259–81. [12] Ching J, Beck JL. Two-step Bayesian structure health monitoring approach for IASC-ASCE phase II simulated and experimental benchmark studies, Technical Report EERL 2003-02. California Institute of Technology, Pasadena, CA: Earthquake Engineering Research Laboratory; 2003. [13] Yuen KV, Au SK, Beck JL. Two-stage structural health monitoring approach for phase I benchmark studies. J Eng Mech 2004;130(2):16–33. [14] Ching J, Beck JL. Bayesian analysis of the phase II IASC–ASCE structural health monitoring experimental benchmark data. J Eng Mech 2004;130(10):1233–44. [15] Yuen KV, Beck JL, Katafygiotis LS. Efficient model updating and health monitoring methodology using incomplete modal data without mode matching. Struct Control Health Monitoring 2006;13(1):91–107. [16] Yuen KV. Bayesian methods for structural dynamics and civil engineering. Wiley; 2010.
53
[17] Yin T, Lam HF, Chow HM. A Bayesian probabilistic approach for crack characterization in plate structures. Comput-Aid Civ Infrastruct Eng 2010;25 (5):375–86. [18] Sankararaman S, Mahadevan S. Bayesian methodology for diagnosis uncertainty quantification and health monitoring. Struct Control Health Monitoring 2013;19:88–106. [19] Figueiredo E, Radub L, Worden K, Farrar CR. A Bayesian approach based on a Markov-chain Monte Carlo method for damage detection under unknown sources of variability. Eng Struct 2014;80(1):1–10. [20] Arangio S, Bontempi F. Structural health monitoring of a cable-stayed bridge with Bayesian neural networks. Struct Infrastruct Eng 2015;11(4):575–87. [21] Beck JL. Bayesian system identification based on probability logic. Struct Control Health Monitoring 2010;17(7):825–47. [22] Girolami M, Calderhead B, Vyshemirsky V. System identification and model ranking: The Bayesian perspective. In: Learning and Inference in Computational Systems Biology. MIT Press; 2010. [23] Green PL, Cross EJ, Worden K. Bayesian system identification of dynamical systems using highly informative training data. Mech Syst Sig Process 2015;56–57:109–22. [24] Worden K, Hensman JJ. Parameter estimation and model selection for a class of hysteretic systems using Bayesian inference. Mech Syst Sig Process 2012;32:153–69. [25] Au SK, Zhang FL. Fundamental two-stage formulation for Bayesian system identification. Part I: General theory. Mech Syst Sig Process 2016;66– 67:31–42. [26] Zhang FL, Au SK. Fundamental two-stage formulation for Bayesian system identification. Part II: Application to ambient vibration data. Mech Syst Sig Process 2016;66–67:43–61. [27] Cox RT. Probability, frequency and reasonable expectation. Am J Phys 1946;14 (1):1–13. [28] Jaynes ET. Information theory and statistical mechanics. Phys Rev 1957;106:620–30. [29] Jaynes ET. Probability theory: the logic of science. Cambridge, U.K.: Cambridge University Press; 2003. [30] Zhou H, Ni Y, Ko J. Eliminating temperature effect in vibration-based structural damage detection. J Eng Mech 2011;137(12):785–96. [31] Sohn H, Allen DW, Worden K, Farrar CR. Structural damage classification using extreme value statistics. J Dyn Sys Meas Control 2005;127(1):125–32. [32] Tipping ME. Sparse Bayesian learning and the relevance vector machine. J Mach Learning Res 2001;1:211–44. [33] Wipf D, Palmer J, Rao B. Perspectives on sparse Bayesian learning. Adv Neural Inf Process Syst 2004;16. [34] Ji S, Xue Y, Carin L. Bayesian compressive sensing. IEEE Trans Sig Process 2008;56(6):2346–56. [35] Huang Y, Beck JL, Wu S, Li H. Bayesian compressive sensing for approximately sparse signals and application to structural health monitoring signals for data loss recovery. Probabilistic Eng Mec. 2016. doi: http://dx.doi.org/10.1016/j. probengmech.2016.08.001. [36] Huang Y, Beck JL, Wu S, Li H. Robust Bayesian compressive sensing for signals in structural health monitoring. Comput-Aid Civ Infrastruct Eng 2014;29 (3):160–79. [37] Mackay DJC. Bayesian methods for adaptive models (Ph.D. Thesis in Computation and Neural Systems). Pasadena, California: California Institute of Technology; 1992. [38] Huang Y, Beck JL. Novel sparse Bayesian learning for structural health monitoring using incomplete modal data. In: Proceedings of 2013 ASCE International Workshop on Computing in Civil Engineering, Los Angeles, CA; June 2013. [39] Huang Y, Beck JL. Hierarchical sparse Bayesian learning for structural health monitoring with incomplete modal data. Int J Uncertain Quantif 2015;5 (2):139–69. [40] Beck JL, Katafygiotis LS. Updating models and their uncertainties. I: Bayesian statistical framework. J Eng Mech 1998;124(4):455–61. [41] Bernal D, Dyke SJ, Lam HF, Beck JL. Phase II of the ASCE benchmark study on SHM. In: Proceedings of 15th eng. mechanics conf., ASCE; 2002. p. 1048–55. [42] Dyke SJ, Bernal D, Beck JL, Ventura C. Experimental Phase II of the Structural Health Monitoring Benchmark Problem. In: Proceedings of 16th Eng. Mechanics Conf., ASCE, Reston, VA; 2003. [43] Beck JL, Jennings PC. Structural identification using linear models and earthquake records. Earthquake Eng Struct Dyn 1980;8(2):145–60. [44] Beck JL, May BS, Polidori DC. Determination of modal parameters from ambient vibration data for structural health monitoring. In: Proceedings of first world conference on structural control, Pasadena, California; 1994. p. TA3:3–TA3:12.