Damage detection inevitably involves uncertainties originated from measurement noise and modeling error. It may cause incorrect damage detection results if not appropriately treating uncertainties. To this end, vibration-based Bayesian model updating (VBMU) is developed to utilize vibration responses or modal parameters to identify structural parameters (e.g., mass and stiffness) as probability distribution functions (PDF) and uncertainties. However, traditional VBMU often assumes that mass is well known and invariant because simultaneous identification of mass and stiffness may yield an unidentifiable problem due to the coupling effect of the mass and stiffness. In addition, the posterior PDF in VBMU is usually approximated by single-chain based Markov Chain Monte Carlo (MCMC), leading to a low convergence rate and limited capability for complex structures. This paper proposed a novel VBMU to address the coupling effect and identify mass and stiffness by adding known mass. Two vibration data sets are acquired from original and modified systems with added mass, giving the new characteristic equations. Then, the posterior PDF is reformulated by measured data and predicted counterparts from new characteristic equations. For efficiently approximating the posterior PDF, Differential Evolutionary Adaptive Metropolis (DREAM) Algorithm are adopted to draw samples by running multiple Markov chains parallelly to enhance convergence rate and sufficiently explore possible solutions. Finally, a numerical example with a ten-story shear building and a laboratory-scale three-story frame structure are utilized to demonstrate the efficacy of the proposed VBMU framework. The results show that the proposed method can successfully identify both mass and stiffness, and their uncertainties. Reliable probabilistic damage detection can also be achieved. Keywords: Vibration-based Bayesian model updating; Uncertainties; The coupling effect of mass and stiffness; Added mass; Differential Evolutionary Adaptive Metropolis (DREAM). 1. Introduction Civil infrastructures are inevitably exposed to different types of structural damages (e.g., aging, deterioration, shrinkage, corrosion, and degradation, etc.) over their lifespan due to internal factors like deficient structural design, imperfect construction, and material defects; and external factors like natural disasters, environmental change, and man-made disasters [1]. Visual inspection is a common way to evaluate structural health and integrity, usually conducted every two years to comply with the federal and local regulations. However, visual inspection results predominantly rely on inspectors’ experience and knowledge, yielding a subjective decision. In addition, visual inspection is sometimes inaccessible to large and complex structures; its implementation is costly and inefficient due to traffic closure [2]. Therefore, developing a robust and efficient damage detection strategy is demanding and practically significant. Vibration-based damage detection (VDE) has been extensively investigated due to the advancement of the measurement and acquisition of vibration signals at a low cost [3]. The main principle of vibration-based SHM is that damage in structural parameters (e.g., mass and stiffness) can be reflected on the change in vibration responses (e.g., acceleration and displacement) and modal parameters (e.g., natural frequency, damping ratio, and mode shape) [4]. Generally, VDE can be categorized into two groups: response-based and model-based methods. Response-based methods directly or indirectly interpret dynamic responses without computer-simulated models (nonparametric approach). Damage is identified by changes in modal parameters between healthy and damaged conditions [5]. However, response-based methods have a critical drawback of only detecting damage location and cannot quantify damage severity. Also, measured 2 responses from the field are always limited. Model-based methods use measured data to update the parameterized computer-simulated models so that damage can be detected, localized, and qualified by the variation in structural parameters [6]. A comprehensive overview of VDE can be found in [7]. Traditional VDE is a deterministic estimation that only yields a single damage quantification solution and ignores uncertainties originating from measurement noise and modeling error. A false alarm in condition evaluation may be conducted if underlying uncertainties are not correctly estimated. Vibration-based Bayesian model updating (VBMU) has been widely applied in various structures for damage detection due to its simple mathematical theorem and powerful capability of handling uncertainties. Successful applications of VBMU includes buildings [8, 9], bridges [10, 11], and lab-scaled structures [12, 13]. Beck et al. [14] established the foundation of the Bayesian model updating approach. Further derivatives and modifications [10, 15, 16] extended the Bayesian approach’s capability and efficacy. Bayesian approach characterizes to-be-updated structural parameters as random variables and formulates parametric model updating function within the Bayes’ theorem. The posterior PDF is explicitly built using prior knowledge from engineering judgment and likelihood function consisting of measured data. In short, VBMU utilizes vibration data to identify structural parameters as probability distribution functions (PDF). The PDF naturally serves as a measure of parameter uncertainty, such as mean, and variance can be easily determined. The key strengths of the Bayesian approach are as follows: 1) rationally and reliably handling incomplete experimental data; 2) using Bayes’ theorem, physical model parameters are characterized by the PDF; 3) only repeating straightforward model evaluations to avoid the most inverse problem’s challenges of unidentifiability, ill-posedness, and non-uniqueness [11]. However, classical VBMU usually assumes that mass is well known and invariant; only stiffness is identified [10, 15], believing that mass is less critical. However, this is not always valid when noticeable variation in mass occurs. Furthermore, the classical characteristic equation is given by: (𝐊−𝝀𝐌)𝝓=0 (1) where 𝐌 is a mass matrix; 𝐊 is a stiffness matrix; 𝝀 are eigenvalues (square of natural frequencies); 𝝓 are eigenvectors (mode shapes). Eq. (1) shows that structural parameters in mass and stiffness are coupled concerning the natural frequency and mode shape. Therefore, simultaneous identification of mass and stiffness can be defined as an unidentifiable problem because an infinite combination of mass and stiffness exists and gives the same natural frequency (herein, the coupling effect of mass and stiffness) [17]. To avoid this issue, mass is usually well-estimated or exactly known for updating stiffness in traditional VBMU due to the availability of the mass information in a deterministic manner (e.g., dimensions). The coupling effect is recently successfully addressed by updating mass and stiffness with an acceptable level. However, it still requires certain prior information and remains a challenge to quantify parameter uncertainties. Xu et al. [18] proposed a time-domain nonlinear restoring force to identify mass and stiffness; however, an external force is required. Zhang and Li [19] presented a loop substructure identification method for mass and stiffness, while mass magnitude at sensor location should be known. Do and Gül [20] established a time series-based model to identify mass and stiffness features. Using incomplete measured data, Lei et al. [21] employed an extended Kalman filter (EKF) to determine the mass-stiffness coupled coefficient. Xu et al. [22] used cross-modal energy sensitivity with a cross-model cross-mode method to identify mass and stiffness changes; Khanmirza et al. [23] applied soft-computing methods to identify mass-damping-stiffness for multistory shear building. In addition to these methods, Rezaiee-Rajand et al. [24] developed a sensitivity-based method with modal kinetic and modal strain energy to update mass and stiffness. Nevertheless, the aforementioned methods are deterministic estimations and cannot provide uncertainties of structural parameters. Recently, Zeng and Kim [16] proposed a new VBMU framework to identify both mass and stiffness and their uncertainties for 2D and 3D numerical shear frames and demonstrated that classical VBMU cannot update both mass and mass stiffness, simultaneously. The coupling effect of mass and stiffness is successfully addressed using two sets of vibration data acquired from two systems: original and modified with added known mass. The new prior distribution, incorporating the new characteristic equation with mass addition, is derived, giving a new posterior PDF. Mass addition strategies are practically achievable, 3 such as moving vehicles on the bridge [25, 26] and adding stationary weights to the structure [27]. Zeng and Kim [16] employed the asymptotic approximation method [14] to circumvent high-dimensional integrals involved in the posterior PDF for Bayesian inference. To be specific, the objective functions are obtained by taking the negative logarithm of the posterior PDF. The analytical formulations of optimal model parameters are derived by the linear optimization method; associated uncertainties are quantified by an inverse Hessian matrix of the objective function. Finally, structural parameters are iteratively updated using modal parameters, e.g., frequency and mode shape. However, the asymptotic approximation method assumes that parameters have unimodal and Gaussian distribution that does not necessarily guarantee an actual physical model when a high level of modeling error and measurement noise occur in practice, especially for multi-modal and non-Gaussian posterior [9, 11]. Also, an insufficient amount of data and complex model class may lead to an unidentifiable problem. One promising way to solve multi-modal and unidentifiable problems is using Markov Chain Monte Carlo (MCMC) to generate samples to approximate the posterior PDF. The high-dimensional integrals in the Bayesian approach can be reasonably calculated. Another attractive feature in MCMC does not require the assumption on the physical model and accurately represents the posterior PDF. Various MCMC techniques have been developed for posterior distribution sampling, such as Metropolis-Hastings (MH) algorithms [28], Gibbs sampling [29], Hamiltonian Markov chains [30], and delayed rejection adaptive Metropolis (DRAM) [8, 11]. These methods adopt a single Markov chain to draw samples, it has demonstrated a limited capability to treat high-dimensional, multi-modal, and flat manifold PDFs. Therefore, they have a relatively low convergence rate and cannot guarantee adequate exploration in parameter space for a target PDF [31]. This work proposed using the Differential Evolution Adaptive Metropolis (DREAM) [31] algorithm to proceed with the distribution estimate. DREAM is essentially a multi-chain MCMC that runs different paths in parallel to target the posterior PDF. It combines different powerful strategies, including a genetic algorithm for population evolution [32], self-adaptive randomized subspace sampling, and outlier chain detection [33], to quickly achieve convergence and seek the best solution by running multiple Markov chains. A wide range of applications has shown that DREAM exhibits excellent performance for complex problems with high-dimensionality, nonlinearity, numerous peaks, and large uncertainties in different research fields, including hydrology [34, 35], chemistry [36, 37], geophysics [38, 39] and renewable energy technique [40], etc. However, to the authors’ best knowledge, DREAM has not been investigated in SHM for civil infrastructures. The current study attempts to explore the efficacy of DREAM in VBMU. This paper presents a new VBMU by treating both the mass and stiffness as equally important. The proposed VBMU intrinsically addressed the coupling effect of mass and stiffness by two sets of vibration data measured from the original and modified system with added known mass. The new characteristic equations are constructed to substitute the traditional one in Eq. (1). The posterior PDF is reformulated by measured modal data and predicted counterparts from the new characteristic equations. The DREAM algorithm is then employed to generate samples for approximation of the posterior PDF. The proposed VBMU simultaneously identifies the mass and stiffness; their uncertainties are also straightforward provided by the estimated PDF. A numerical study on a ten-story shear building and an experimental study on a three-story aluminum frame small-scale model are used at intact and damaged structural states to verify the accuracy and feasibility of the proposed method. The outline of the presented work is listed as follows. The background of classical VBMU and a brief review of the Bayesian approach’s theory are first described in Section 2, along with model parameterization. Section 3 presents the methodology of the proposed VBMU, in which the new characteristic equations, strategies of mass addition, and DREAM algorithm are introduced explicitly. Section 4 shows one illustrative example to validate the methodology using a numerical example, followed by the validation of laboratory-scale testing. Probabilistic damage detection is also performed. Finally, conclusions and summaries are provided in Section 5. 4 2. Classical vibration-based Bayesian model updating approach In this section, the structural model's parameterization is first introduced in Section 2.1. The background of classical vibration-based Bayesian model updating is then presented in Section 2.2. The challenge of updating mass and stiffness using the classical Bayesian approach is also discussed. 2.1. Model parametrization Assuming a structural model is linear and can be parameterized by model parameters with Degree-of-freedoms (DOFs), 𝑁!, and defined model class, 𝐶. A commonly used parameterization of stiffness matrix, 𝐊(𝜽), and mass matrix, 𝐌(𝜷), could be described as follows [10]: 𝐊(𝜽)=𝐊"+∑𝜃#𝐊#$!#%& 𝐌(𝜷)=𝐌"+∑𝛽'𝐌'$"'%& (2) where 𝜽=[𝜃&,𝜃(,⋯, 𝜃$!]) and 𝜷=7𝛽&,𝛽(,⋯, 𝛽$"8)are stiffness and mass parameters vectors, respectively, which will be updated to let the FE model agrees with the real structure using vibration measurements. 𝑁* and 𝑁+ are the number of stiffness and mass parameters, respectively; The 𝑙th stiffness parameter forms the 𝑙th elemental stiffness matrix, 𝐊#=𝜕𝑲/𝜕𝜃#; similarly, the 𝑚th mass parameter forms the 𝑚th elemental mass matrix, 𝐌'=𝜕𝑴/𝜕𝛽'. Note that 𝐊" and 𝐌" are defined as constant matrices that are not dependent on model parameters. For the sake of simplicity, 𝐊" and 𝐌" are set as zero. The parameterization in Eq. (2) allows to update the model at an elemental level and capture local variations in stiffness and mass, as each parameter characterizes an elemental matrix. This also leads to a more reasonable and accurate damage detection with different levels of deterioration. 2.2. Background of classical VBMU The strength of the Bayesian model updating approach lies in that it uses both the prior information (existing structural knowledge) and measured data (new structural knowledge) to estimate the posterior PDF. In other words, the Bayesian approach updates the prior PDF by measured data, yielding the posterior PDF. Consider a model class 𝐶 which represents the relationship between measured data 𝐷 and parameters of interest 𝜴, containing stiffness parameter 𝜽 and mass parameter 𝜷. Then, by following Bayes’ theorem, the posterior PDF of uncertainty parameter 𝜴 is given by [10]: 𝑃(𝜴|𝐷,𝐶)=𝑃(𝜴|𝐶)𝑃(𝐷|𝜴,𝐶)𝑃(𝐷|𝐶) (3) where 𝑃(𝜴|𝐶) is the prior PDF, reflecting the probability of uncertainty parameters without observed data. In many cases, the selection of prior PDF depends on engineers’ judgment and physical meaning. The uniform distribution is widely used as the uninformative prior PDF to ensure the measured data entirely dominates Bayesian inference and minimizes the effect of prior information. The term of 𝑃(𝐷|𝐶) is a normalizing constant so that the posterior PDF can be integrated to unity over the parameter space, which is given by, 𝑃(𝐷|𝐶)=C𝑃(𝜴|𝐷)𝑃(𝐷|𝜴,𝐶) 𝑑𝜴 (4) The likelihood function, 𝑃(𝐷|𝜴,𝐶), describes how likely the measurements are reproduced from a model parameterized by a set of 𝜴. Considering an uninformative prior PDF, the posterior PDF is proportional to the likelihood function: 𝑃(𝜴|𝐷,𝐶)=𝑐"𝑃(𝐷|𝐶) (5) where 𝑐" represents a constant value to reflect both 𝑃(𝐷|𝐶) and 𝑃(𝜴|𝐷). Generally, for vibration-based system identification, the common measured data in the likelihood function consists of measured natural frequencies and mode shapes. Then, two error functions (EF) of a 5 given one mode, 𝑚, are adopted to formulate the likelihood function, namely frequency EF and mode shape EF [15]. Frequency EF, 𝜀,,', is defined as: 𝜀,,'=𝑓H'−𝑓'(𝜴) (6) where 𝑓H' is the 𝑚th measured frequency, 𝑓'(𝜴) is the 𝑚th calculated frequency in a model given a set of 𝜴. Mode shape EF, 𝜀'.,', is defined as: 𝜀'.,'=𝝓I'−𝑳"𝝓'(𝜴) (7) where 𝝓I' and 𝝓'(𝜴) are measured mode shape and calculated one of the 𝑚th mode, respectively. 𝑳" consists of '1s' or '0s' to match measured partial mode shapes with theoretical counterparts. Note all mode shapes are normalized to unity norm to map them in the same context. With the assumption that 𝜀,,' and 𝜀'.,' follow zero-mean Gaussian distribution, then the posterior PDF in Eq. (3) is rewritten as follows: 𝑃(𝜴|𝐷,𝐶)=𝑐"exp N−12𝜅(𝐽(𝜴)S (8) The objective function, Eq. (9), can evaluate the accuracy of predicted natural frequency and mode shape obtained from new characteristic equations against the measured data. 𝐽(𝜴)= TUV𝑓H'−𝑓'(𝜴)W(+NX𝝓I'−𝑳"𝝓'(𝜴)Y)X𝝓I'−𝑳"𝝓'(𝜴)YSZ/'%& (9) where 𝜅 is an uncertainty parameter of prediction error. In the current study, the variances of the measured frequency and mode shape are used as 𝜅(. 𝜅 consists of 𝜎,,' and 𝜎'.,'; 𝜎𝑓,𝑚 and 𝜎'.,' are the standard derivation of the 𝑚th measured frequency and mode shape, respectively. These two weighting factors can be identified by either Bayesian modal analysis [41] or stochastic subspace identification (SSI) based uncertainty analysis [42], rather than manually tunning. For avoiding intractable high-dimensional integrals, MCMC is employed to approximate the posterior PDF in Eq. (8) without any assumption on a model by iteratively drawing samples from the target distribution. Classical VBMU calculates theoretical frequency, 𝑓'(𝜴), and mode shape, 𝝓'(𝜴) in Eq. (9), given a set of 𝜴 using the classical characteristic equation in Eq. (1). Understandably, simultaneous updating stiffness and mass yield an unidentifiable problem due to the coupling effect of mass and stiffness. The infinite sets of mass and stiffness derive the same frequency so that correct model updating cannot be achievable. The new characteristic equations with added mass will substitute classical ones and address the coupling effect in the next section. 3. Formulation of a new vibration-based Bayesian model updating approach New characteristic equations with added known mass are first presented in Section 3.1 to address the coupling effect of mass and stiffness for updating mass and stiffness. The mass-adding strategies are discussed in Section 3.2, including the number, location, and magnitude of added mass. The DREAM algorithm, a multi-chain MCMC to approximate the posterior PDF, is presented in section 3.3. 3.1. New characteristic equations with added mass Classical VBMU adopts a classical FE model, e.g., (𝐊−𝝀𝐌)𝝓=0 to update stiffness with known mass. In the proposed method, two groups of measured data acquired from original and modified systems with added mass ∆𝒎, are used to derive the new characteristic equations. This aims to address the coupling effect of mass and stiffness. The original and modified systems with added mass, ∆𝒎, are merged into one equation based on the fundamentals of structural dynamics. The core idea of addressing the coupling effect of mass and stiffness 6 is to eliminate either mass or stiffness when updating each of them. For example, first, characteristic equations for the original and modified systems are expressed as: 𝑲𝝓=𝑴𝝓𝝀 (10) 𝑲𝝓2=(𝑴+∆𝒎)𝝓2𝝀′ (11) where 𝝀 and 𝝓 are eigenvalues (square of natural frequencies) and mode shapes before modification; 𝛌′ and 𝝓2 are eigenvalues and mode shapes after modification. Second, the new eigen-equation error with added mass when updating mass is derived as: (𝝀−𝝀2)3&𝛌2𝝓2)∆𝒎𝝓−𝝓2)𝑴𝝓=𝟎 (12) For the sake of simplicity, more details can be found in [16]. Eq. (12) can be rewritten as: V𝛌2𝝓2)∆𝒎−𝝀𝝓2)𝑴+𝛌2𝝓2)𝑴W𝝓=𝟎 (13) Define 𝐀=𝛌2𝝓2)∆𝒎+𝛌2𝝓2)𝑴 , 𝐁=𝝓2)𝑴, then Eq. (13) is expressed as: (𝐀−𝝀𝐁)𝝓=𝟎 (14) Similarly, when updating stiffness, the new eigen-equation error is shown as [16]: (𝝀23𝟏−𝝀3𝟏)3𝟏𝝓2)∆𝒎𝝓−𝝓2)𝑲𝝓 (15) Eq. (15) can be rewritten as: V𝝀23𝟏𝝓2)𝑲−𝝀3𝟏𝝓2)𝑲−𝝓2)∆𝒎W𝝓=𝟎 (16) Define 𝐄=𝝀23𝟏𝝓2)𝑲−𝝓2)∆𝒎 , 𝐅=𝝓2)𝑲, then Eq. (16) is expressed as: (𝐅−𝝀𝐄)𝝓=0 (17) Eqs. (14) and (17) are defined as the new characteristic equations to replace the classical ones in Eq. (1). It is noted that the two new characteristic equations have the same formats as the generalized eigenvalue problem, 𝝀 and 𝝓 can be easily solved in mathematics or solver in the computer program, such as ‘eig’ function in MATLAB. Two new characteristic equations eliminate the coupling effect of mass and stiffness. For example, mass updating by using Eq. (14) does not require any stiffness information. Likewise, updating stiffness does not require any mass information by using Eq. (17). For output-only modal analysis, the mode shapes are not mass-normalized, and only unscaled mode shapes are identified because of unknown excitation forces. Before the model updating, the measured mode shapes have to be normalized by either the mass-change scaling method [27] or stiffness-change scaling method [43] to ensure measured and predicted mode shapes are comparative. In this paper, the mass-change scaling method is adopted to calculate scaled mode shapes. In addition, because only a few DOFs are available related to the sensor location, limited sensors in practice usually lead to incomplete measured mode shapes. Therefore, mode shape expansion techniques [44] can expand measured mode shapes to complete mode shapes of full DOFs. 3.2. Strategy of adding mass The optimized mass-change strategy has been comprehensively discussed in [27], including mass magnitude, number of added mass, and locations of added mass. Generally, two criteria for creating a modified system with added known mass are required: Step 1) noticeable frequency change is observed between the original system and modified system; Step 2) mode shapes after modification change slightly. The frequency change and mass addition are correlated by natural frequencies in the original and modified systems. Considering a structure with multiple DOFs, the relation between added mass and frequency shift can be expressed as [27]: ∆𝑓𝑓=1−e11+∆𝑀𝑀∗ (18) where ∆𝑓 is the frequency change after adding mass (= 𝑓2−𝑓); 𝑓 and 𝑓′ are the natural frequencies in the original and modified systems, respectively; ∆𝑀=𝝍)∆𝒎𝝍; 𝑀∗=𝝍)𝑴𝝍, where 𝝍 is unscaled mode shape in the original system; ∆𝒎 is a diagonal matrix with main diagonal are added mass; 𝑴 is a mass 7 matrix in the original system. Eq. (18) allows us to determine expected added mass when the terms of ∆𝑓, 𝑓, and 𝑀∗ are known; 𝑓 and 𝝍 are identified by modal analysis, the analytical mass is used as 𝑀∗. Based on [27], we select a frequency ratio, 𝑓/𝑓′ to determine ∆𝑀 in Eq. (18). Note that the selection of ratio depends on the expected accuracy in modal analysis and mode shape normalization. Finally, the magnitude of added mass can be estimated using Eq. (18). The number of added mass depends on the number of modes to identify in modal analysis. Ideally, the added or attached mass should be as many as possible. López-Aenlle et al. [27] recommended that the number of added mass should be at least the number of peaks and valleys of each mode shape. To optimize the location of added mass, the most significant frequency shift can occur when the mass is attached to the peaks and valleys of mode shape, while the frequency shift is minimal when mass is attached to the nodal positions. 3.3. DREAM algorithm The posterior PDF needs high-dimensional integrals that is impractical for complex structures. In the present work, the DREAM algorithm proposed by Vrugt et al. [33] is used to approximate the posterior PDF by generating samples based on a differential evolutionary algorithm. Compared to other single-chain MCMC methods, the DREAM has the appealing feature of running multiple chains simultaneously to explore global solutions. The DREAM algorithm uses randomized subspace sampling to automatically tune the mean and variance of the proposal distribution. Therefore, it is highly robust to the selection of the prior distribution. The removal of outlier chain and crossover schemes are also used to expedite convergence to a target distribution. Practical applications exhibited high efficiency and accuracy in the sampling for the problems having high-dimensionality, nonlinearity, numerous peaks, and local optima. Theoretical background and detailed MATLAB procedures in DREAM can be found in [31, 33]. The flowchart of the DREAM algorithm is also shown in Fig. 1. The main implementation steps of the DREAM algorithm are summarized as follows: Fig. 1. Flowchart of DREAM algorithm Step 1: Initialize the problem dimension 𝑁, the number of Markov chains 𝑃, unknown parameter vector, 𝜴67 (𝑖=1,2,3⋯,𝑁; 𝑗=1,2,3⋯,𝑃), and the maximum iteration, 𝐼'89. 𝑁 individual 8 samples for each chain are randomly generated from the selected prior distribution as initial values, such as 𝜴6&, 𝜴6(, ⋯, 𝜴6:. Step 2: A mutation operation is performed to generate candidate samples at each parameter sample of each iteration for the 𝑘th Markov chain. Crossover operation is then used to iteratively update current candidate samples from the mutation process based on crossover probability 𝐶𝑅 within the range of [0,1]. Step 3: Calculate the posterior probability and acceptance rate of updated candidate samples at the 𝑠th iteration: 𝛼V𝜴6,.7,𝑤6,.;&7W=⎩⎨⎧minz𝑝V𝑤6,.;&7|𝐷W𝑝V𝜴6,.7|𝐷W,1} ; 𝑝V𝜴6,.7|𝑦W>01 ; 𝑝V𝜴6,.7|𝑦W<0 (19) where 𝜴6,.7 and 𝑤6,.;&7 are the samples at the 𝑠th iteration and (𝑠+1)th iteration, respectively; 𝛼V𝜴6,.7,𝑤6,.;&7W is the acceptance rate; 𝑝V𝑤6,.;&7|𝑦W and 𝑝V𝜴6,.7|𝑦W are the posterior probability of 𝑤6,.;&7 and 𝜴6,.7, respectively. 𝐷 is the measured data. Step 4: Determine whether accepting or rejecting the samples of 𝑤6,.;&7. If 𝛼V𝜴6,.7,𝑤6,.;&7W>𝑢, 𝑢 is randomly generated from a uniform distribution 𝑈(0,1). Then, accept a new sample of 𝑤6,.;&7, otherwise reject and keep the iteration. Step 5: Remove the outlier chain using the Inter-Quartile-Range (IQR) statistical method [33]. Specifically, ℋ is firstly defined as the mean of the logarithm of the posterior distribution of the last half samples in each chain. ℋ=𝑄<−𝑄& is calculated, where 𝑄& and 𝑄< are the lower and upper quartile of the 𝑃 chains. Chains with ℋ<𝑄&−2∙IQR are detected as aberrant ones. Note removal of outlier chain is necessary, as outlier chains will impair the distribution estimate and slow down the evolution so that reaching a good convergence is impossible. In addition, outlier chains frequently present in high-dimensional problems and tend to be stuck in local optima, resulting in a biased distribution [33]. Step 6: The iteration process stops when Markov chains converge to the target posterior distribution. Otherwise, repeat steps 2-5. DREAM algorithm uses Gelman-Rubin statistics, scale reduction factor 𝑅.=8= [45], as a convergence criterion to determine whether the calculation terminates or not. DREAM algorithm stipulates that if 𝑅.=8=<1.2 for all unknown parameters, a stable posterior PDF is achieved. Note the value of 1.2 has been demonstrated as a robust indication to officially declare stationary and reliable convergence [33]. 𝑅.=8= has an expression as follows: 𝑅.=8== e𝛾−1𝛾+𝑃+1𝑃∙𝑍𝐵𝛾 (20) where 𝛾 is the number of iteration samples of each chain; 𝑃 is the number of Markov chains used for sampling; 𝑍 is the mean of the variance of total 𝑃 Markov chains; the ratio of 𝐵/𝛾 is the variance of the mean of 𝑃 parallel Markov chains. In summary, the proposed VBMU addresses the coupling effect of mass and stiffness by using two sets of measurements from the original and modified system with added mass. Two new characteristic equations (herein, Eq. 14 and 17) substitute the classical one (Eq. 1). Fig. 2 shows the flowchart of the proposed method. First, the natural frequencies and mode shapes of the original and modified system are identified using the output-only modal analysis method. Note that mode shapes need to be normalized by the mass-change scaling method before updating mass and stiffness (described in Section 3.1). Second, the objective functions in Eq. (9) with measurements in the original system are used to measure the accuracy of analytical frequencies and mode shapes satisfying with new characteristic equations, e.g., Eqs. (14) and (17). Third, the DREAM is used to approximate the posterior PDF and estimate the quantity of interests (PDF, mean, 9 and variance) (described in Section 3.3.). The procedures of updating mass and stiffness are independent and individually implemented. Therefore, the coupling effect has been removed in the entire updating process. Fig.2. The flowchart of the proposed Bayesian model updating 3.4. Probabilistic damage detection The quantities of interest obtained from samples are useful to infer the probability of variation in mass and stiffness parameters. The mean of updated parameters and associated uncertainties can be used to quantify the damage probability given a specific damage fractional level, 𝑑, compared with the structural healthy condition. Based on asymptotic Gaussian approximation, the probability of damage in the 𝑙th structural parameter can be given by [14]: 𝑃#!8'(𝑑)=𝑃V𝜴#>!<(1−𝑑)𝜴#?!|𝐶W =C𝑃V𝜴#>!<(1−𝑑)𝜴#?!|𝜃#?!,𝐶W𝑝(𝜴#?!@3@|𝐶)𝑑𝜴#?! (21) » Φ(&3!)𝜴#∗%&3𝜴#∗'&D(&3!)((E#%&)(;(E#'&)( where Φ(∙) represents the standard Gaussian cumulative distribution function for random variables; 𝜴# and 𝜎# represents the optimal value of the 𝑙=F structural parameter and its standard derivation, respectively. Superscripts, 𝑢𝑑 and 𝑝𝑑, represent the undamaged and possibly damaged structural state, respectively. 4. Illustrative examples The efficacy of the proposed VBMU is evaluated by a numerical example in Section 4.1, followed by an experimental test with a laboratory-scale three-story shear frame in Section 4.2. 4.1. Numerical example: a ten-story shear building The ten-story shear building sketch is shown in Fig. 3, modeled as a ten-DOFs structure. Assume the connection between column and floor is rigid; mass and stiffness at each floor are uniformly distributed. 10 Also, suppose one sensor is installed on each floor to measure all modal displacements in each mode shape. Lumped mass is used and taken as M6=25 kg, 𝑖=1,2,⋯,10. While the inter-story stiffness at each floor is taken as K6=1.5×10G N/m, 𝑖=1,2,⋯,10. We define stiffness coefficient (SC) as 𝜃6=K68/K6, and mass coefficient (MC) 𝛽6=M68/M6, where K68 and K are the 𝑖th actual and theoretical stiffness, respectively; M68and M are the 𝑖th actual and theoretical mass, respectively, resulting in a total of 20 coefficients to be updated. The FE model of this shear building is constructed based on fundamental structural dynamics using MATLAB. The natural frequencies and mode shapes for the original structure can be obtained by the eigenvalue problem so that the first six natural frequencies are 5.827, 17.350, 28.486, 38.985, 48.613, and 57.156 Hz. To create a modified structure, we first select a frequency ratio, 𝑓/𝑓′, of 1.02. Using Eq. (18), the magnitude of added mass can be estimated as 1 kg. For the sake of simplicity, each floor has the equivalent mass addition by concrete blocks with the weight of 1 kg, as shown in Fig. 3 (b). Gaussian white noise with zero-mean and 1% coefficient of variation (COV) is considered and added to the exact frequency and mode shape for all the modes of interest. Note that the value of 1% could reflect realistic noise level in practice based on modal test by Bayesian modal analysis [46]. Mass and stiffness are updated by two sets of simulated measured data acquired from original and modified systems. (a) (b) Fig. 3. Ten-story shear building: (a) original structure; (b) modified structure with mass addition (concrete block) 4.1.1. Finite model updating In the first case, no modeling error is assumed between the actual structure and the FE model. Also, the structure is healthy by setting all 𝜃 and 𝛽 as unity. The first six modes are assumed as available measured data. The DREAM algorithm is used to generate samples for estimating of the posterior PDF. Every sample will yield the analytical frequencies and mode shapes using new characteristic equations. Initial settings in DREAM are defined as: ten Markov chains are run parallelly with 6,000 samples per chain; initial values for 10 SCs and 10 MCs are set as a range of [0.5 1.5]. (a-i) (a-ii) 11 (b-i) (b-ii) Fig. 4. Trace plots and 𝑅.=8= in healthy scenario: (a-i): Trace plots of ten MCs; (a-ii) Trace plots of ten MCs SCs; (b-i): 𝑅.=8= of MCs (b-ii): 𝑅.=8= of SCs Fig. 4 shows the results of updated coefficients. Fig. 4 a-i and a-ii are trace plots of one Markov chain that show how each MC and SC are updated with samples, respectively. All the parameters achieved a stable state. Fig. 4 b-i and b-ii display the variation of the convergence diagnosis for mass and stiffness updating, respectively. The scale reduction factor, 𝑅.=8=, assesses whether the Markov chain converges or not. The 𝑅.=8= of each parameter quickly decays below 1.2, satisfying DREAM's convergence criterion and attains the stationary posterior distribution. The last 30,000 samples herein of ten Markov chains are used to calculate the quantity of interests of all parameters, such as mean and standard derivations. The results of updated coefficients are listed in Table 1, including mean and standard derivation (S.D.). The identified mean values exhibit an excellent agreement with actual counterparts. The errors and standard derivations for all coefficients are small; the maximum error of 2.08% is observed. The histograms of the marginal distribution of ten SCs and MCs are shown in Fig. 5; red curves represent a fitted distribution based on mean and standard derivation. Each histogram has a clear peak and is well-approximated by Gaussian distribution. Overall, each parameter is reasonably identified as the correct values and has a fairly good convergence. (a) (b) Fig. 5. Histograms of updated coefficients: (a) MCs; (b) SCs 12 Table 1. Results of updated coefficients Coefficients Actual Updated Mean S.D. (%) Error (%) 𝛽& 1.0000 0.9831 1.31 1.69 𝛽( 0.9922 1.27 0.78 𝛽< 0.9997 1.81 0.03 𝛽H 1.0013 1.71 0.13 𝛽I 1.0132 2.01 1.32 𝛽G 1.0102 1.49 1.02 𝛽J 1.0023 1.68 0.23 𝛽K 0.9795 2.02 2.05 𝛽L 1.0051 1.56 0.51 𝛽&" 1.0016 1.29 0.16 𝜃& 0.9842 1.01 1.58 𝜃( 0.9911 1.66 0.89 𝜃< 0.9956 1.63 0.44 𝜃H 0.9899 1.48 1.01 𝜃I 0.9972 1.87 0.28 𝜃G 0.9963 2.19 0.37 𝜃J 0.9830 1.49 1.70 𝜃K 0.9914 2.79 0.86 𝜃L 1.0132 1.90 1.32 𝜃&" 0.9792 1.69 2.08 Table 2. Results of updated frequencies and MAC values Mode No. Actual Frequency (Hz) MAC Updated Error (%) 1 5.827 5.803 0.41 1.0000 2 17.350 17.275 0.43 1.0000 3 28.486 28.490 0.02 0.9999 4 38.985 38.920 0.17 0.9998 5 48.613 48.580 0.07 0.9996 6 57.156 56.854 0.53 0.9985 7 64.422 64.171 0.39 0.9973 8 70.248 70.225 0.03 0.9979 9 74.506 74.521 0.02 0.9987 10 77.099 76.936 0.21 0.9988 The updated frequencies and MAC values are summarized in Table 2. It is observed that updated frequencies are almost the same as actual ones; the relative error is less than 1%. The values of the Modal assurance criterion (MAC) [47] that reflect the similarity of updated and actual mode shapes are also close to unity. It is worth mentioning that the higher modal parameters from the 7th to 10th order are not used in the updating process, but they are still successfully identified. 4.1.2 Probabilistic damage detection In the second case, the probabilistic damage detection is performed to detect simulated damage location and extent by the proposed VBMU. Damage extent is defined as the change in mass/stiffness coefficients at a specific floor. The damage scenario considered in this example is shown in Table 3. The negative sign 13 denotes the reduction of mass/stiffness. The model at healthy condition is assumed known. The unity of MCs and SCs represents a healthy structural state. Table 3. Damage scenario Parameters Mass Stiffness Reduction in percent (location) -10% (2nd floor), -20% (6th floor), -30% (9th floor) -10% (3rd floor), -20% (6th floor), -30% (9th floor) The same modified system is created by adding mass as described in Section 4.1.1, and the same measurement points and vibration data are selected to identify the damage, e.g., the first six frequencies and mode shapes. In addition, the initial settings in DREAM are the same as the healthy example in Section 4.1.1. Fig. 6 shows the results of damage detection. Fig. 6 a-i and a-ii are trace plots of one of ten Markov chains for MCs and SCs, respectively, visualizing that all the coefficients stably converge. In Fig. 6 b-i and b-ii, the convergence criterion, 𝑅.=8= is less than 1.2, indicating the stationary Markov chains are achieved. The last 30,000 samples are used to calculate the mean and standard derivation of all coefficients. (a-i) (a-ii) (b-i) (b-ii) Fig. 6. Trace plots and 𝑅.=8= in damages scenario: (a-i): Trace plots of ten MCs; (a-ii) Trace plots of ten MCs SCs; (b-i): 𝑅.=8= of MCs (b-ii): 𝑅.=8= of SCs Table 4 lists identified coefficients and their standard (S.D.) derivations in damage scenarios. It is observed that all updated MCs and SCs are almost identical to actual values. The maximum errors for all coefficients are less than 2% except 𝜃&" with the error of 2.48%, revealing an outstanding performance in damage localization and quantification on both mass and stiffness. Fig. 7 shows the histograms of MCs and SCs; the red curves are fitted Gaussian distribution based on samples. The Gaussian distribution can desirably approximate the posterior PDF. It is also found that the fitted curves in SCs are relatively wider spreading than those in MCs, demonstrating that identified stiffness has larger uncertainty than mass. Table 4. Results of updated coefficients for damage scenario Coefficients Actual Updated Mean S.D. (%) Error (%) 𝛽& 1.0000 1.0140 1.31 1.40 𝛽( 0.9000 0.8885 1.11 1.28 𝛽< 1.0000 1.0133 1.26 1.33 𝛽H 1.0000 1.0026 1.30 0.26 14 𝛽I 1.0000 1.0130 1.15 1.30 𝛽G 0.8000 0.7932 0.98 0.85 𝛽J 1.0000 0.9876 1.28 1.24 𝛽K 1.0000 1.0126 1.34 1.26 𝛽L 0.7000 0.6887 0.87 1.61 𝛽&" 1.0000 0.9847 1.33 1.53 𝜃& 1.0000 1.0171 1.45 1.71 𝜃( 1.0000 1.0196 1.48 1.96 𝜃< 0.9000 0.8891 1.30 1.21 𝜃H 1.0000 1.013 1.52 1.30 𝜃I 1.0000 1.0147 1.35 1.47 𝜃G 0.8000 0.7865 1.16 1.69 𝜃J 1.0000 1.0175 1.45 1.75 𝜃K 1.0000 1.0177 1.40 1.77 𝜃L 0.7000 0.6943 0.99 0.81 𝜃&" 1.0000 1.0245 1.42 2.45 The updated frequencies and MAC values in the damage scenario are derived using updated MCs and SCs, as shown in Table 5. All errors are less than 1% indicating the efficacy of damage detection. Although incomplete modal information, e.g., only the first six modes were used, all the frequencies and MAC values are identified accurately. (a) (b) Fig. 7. Histograms of updated coefficients for damage scenario: (a) MCs; (b) SCs Table 5. Results of updated frequencies and MAC values for damage scenario Mode No. Actual Frequency (Hz) MAC Updated Error (%) 1 5.937 5.976 0.65 1.0000 2 17.149 17.223 0.44 1.0000 3 27.968 28.083 0.41 0.9999 4 37.528 37.646 0.31 0.9998 15 5 48.866 48.992 0.26 0.9998 6 56.971 57.204 0.41 0.9993 7 66.449 66.822 0.56 0.9992 8 71.327 71.519 0.27 0.9987 9 74.604 74.831 0.30 0.9984 10 77.088 77.434 0.45 0.9985 (a) (b) Fig. 8. Probabilistic damage curves: (a) MCs; (b) SCs The probabilistic damage curves are also plotted using identified coefficients and uncertainty information by Eq. (21), as displayed in Fig. 8. It is found that curves at the damaged location are distinguishable from those at a healthy location by observing the curve’s distance from healthy cases. Furthermore, some quantities can be interpreted from curves. For example, mass on the sixth floor (𝛽G) and stiffness on the ninth floor (𝜃L) have a possible reduction of 20% and 30% with a high probability of 83.2% and 81.6%, respectively. Thus, the proposed VBMU exhibits excellent performance on damage detection on mass and stiffness; both damage location and severity are successfully identified. 4.2. Experimental test: a three-story shear frame The experimental test was performed to verify the accuracy and efficacy of the proposed VBMU for both mass and stiffness identification. A shear building, made of aluminum, has a height and width of 914 and 305 mm., respectively. All the plates and columns have the same geometric properties. The length and width of a plate are both 305 mm with a 25 mm thickness. The column has the length, width, and thickness of 254, 25, and 6 mm, respectively. The initial Young’s modulus and mass density of the aluminum are estimated as 69 GPa and 2,700 kg/m3, respectively. The shear building is modeled as a three-DOF structure using the MATLAB program based on the dimensions and material properties. Free vibration test was performed by inducing the excitation using a rubric hammer. The hammer impacted the structure on the top floor. Horizontal responses were measured by the three IMI 603C01 accelerometers fixed with magnets in the middle of the left side of each floor plate; the associated LabVIEW data acquisition software was used to process the measured signal. In the measurement, ten-second data were recorded with a sampling frequency of 2,000 Hz. The acceleration at each floor is also preprocessed by a low-pass filter with a cut-off frequency of 50.2 Hz, and down sampled to 200 Hz to identify the frequencies of interest and remove noise from high frequencies. The automated stochastic subspace identification (SSI) [48] is used to identify modal parameters, e.g., natural frequencies and mode shapes, and associated uncertainties. Uncertainties on modal parameters measure modal parameters’ accuracy and can be used as weighting factors, such as 𝜅 in Eq. (8). Fig. 9 (a) shows the experimental setup in the laboratory for the original system at the Civil and Environmental Engineering at the University of Louisville. To create the modified system, the ratio of frequency in the original system to that in the modified system is assumed to be 1.04; the magnitude of mass addition is then estimated as 0.545 kg using Eq. (18). Therefore, the concrete block with a weight of 0.545 kg is added to each floor, as shown in Fig. 9 (b). The same measurement and modal identification are carried out for the modified system. Before model updating and damage detection, mode shapes in two systems are normalized 16 by the mass-change scaling method. Similar to the numerical example, it is convenient to use MCs and SCs as updating indices. Each floor has a representative value of MC and SC, giving a total of six coefficients to be updated, e.g., 𝛽&,𝛽(, and 𝛽< (MCs) and 𝜃&,𝜃(, and 𝜃< (SCs) with labeling a subscript number from the bottom floor (1) to the top floor (3). (a) (b) Fig. 9. Test setup of the shear building: (a) original system; (b) modified system with concrete block 4.2.1. Finite model updating In the first case, the natural frequencies and mode shapes in the original and modified system under the healthy state are used to update the model. The initial settings in DREAM are as follows: ten Markov chains are simultaneously run to generate a total of 20,000 samples (2,000 samples per chain); all MCs and SCs have initial values ranging from 0.5 to 1.5. The number of samples designed for the experimental test is less than that in the numerical example, because we have fewer coefficients to update in this test. (a) (b) (c) Fig. 10. Results of updated mass: (a) trace plot; (b) square: the sample mean of each chain, error bar: ±2 standard derivations; (c) convergence diagnosis, 𝑅.=8= (a) 17 (b) (c) Fig. 11. Results of updated stiffness: (a) trace plot; (b) square: the sample mean of each chain, error bar: ±2∙S.D.; (c) convergence diagnosis, 𝑅.=8= Figs. 10 (a) and 11 (a) are the trace plots of the variation of MCs and SCs values, respectively, as samples increases in one Markov chain. The stable convergence of each coefficient is visually observed. The rest of the figures in Figs. 10 and 11 give the updating results over ten Markov chains and convergence diagnosis. The sample mean of MCs and SCs in each Markov chain and ±2∙S.D. are shown in Fig .10 (b) and 11 (b), respectively. The mean value of each coefficient is almost identical to one another among ten chains, indicating that the updating results are reliable and accurate. The convergence diagnosis, 𝑅.=8=, shown in Fig. 10 (c) and 11 (c), are a useful graphical tool to evaluate convergence state. The resulting plots of 𝑅.=8= that quickly decrease below 1.2, indicating that the sampling process is performed to achieve the stationary Markov chain. Herein, the last 10,000 samples are used to calculate the quantities of interest, mean and standard derivation. Fig. 12. Density distribution of updated coefficients The density distribution of each coefficient estimated by the Gaussian kernel estimator (GKE) is also presented in Fig. 12. It can be seen that the density distributions of 𝛽&, 𝛽<, 𝜃&, and 𝜃< are non-Gaussian and multi-modal, indicating in practice, stiffness and mass do not necessarily follow the Gaussian distributions. The estimated distributions illustrate that DREAM is appropriate to approximate the distribution with non-normal shape and multi-peaks. Furthermore, except 𝛽( and 𝜃( which are distributed over a narrow region, the coefficients have a wide-ranging distribution, suggesting they have larger uncertainties (see S.D. in Table 6). Table 6. Results of updated coefficients under healthy condition Coefficients Initial Updated Mean S.D. (%) Change (%) 𝛽& 1.0000 1.2341 10.36 23.41 𝛽( 1.1864 8.91 18.64 𝛽< 1.1237 11.52 12.37 𝜃& 0.7995 16.44 -20.05 𝜃( 1.2452 10.20 24.52 𝜃< 1.0680 17.81 6.80 18 Table 7. Results of updated frequencies and MAC values under healthy condition Mode No. Frequency (Hz) MAC Actual FE model Initial Updated Initial Error (%) Updated Error (%) 1 7.657 8.667 13.20 7.78 1.63 0.9882 0.9979 2 22.47 24.213 7.76 22.23 1.08 0.9938 0.9972 3 33.77 34.856 3.22 34.00 0.68 0.9954 0.9965 Table 6 shows updated MCs and SCs and their S.D. The updated frequencies and MAC values are tabulated in Table 7. All the MCs increase but 𝜃& decreases. The model updating aims to match measured responses with analytical counterparts. In this case, measured frequencies are overall smaller than those in the FE model (see Table 7). From fundamental structural dynamics, frequency is proportional to stiffness but inversely proportional to mass. Therefore, stiffness and mass have to decrease and increase, respectively, to match measured frequencies with those in the FE model. The frequency errors of all modes are significantly reduced and MAC values are updated to be close to 1.0. These values demonstrate satisfactory updating model results. 4.2.2. Probabilistic damage detection Two damage scenarios are intentionally introduced with increasing severity in the shear building by reducing the thickness of column and increasing the weight of the floor, as shown in Table 8, the positive/negative sign denotes the increasing/reduction. The thickness of one column at the second and third floor is reduced by 50%, resulting in a 21.8% stiffness reduction in the corresponding floor; A concrete block with the weight of 1.54 kg is added to the second and third floor to mimic mass change due to damage, which produces 21.5 % mass increase in the corresponding floor. Table 8. Damage scenarios Notation Mass change Stiffness change D1 +21.5% (3rd floor) -21.8% (3rd floor) D2 +21.5% (2nd floor), +21.5% (3rd floor) -21.8% (2nd floor), -21.8% (3rd floor) (a) (b) (c) (d) Fig. 13. Test setup for two damage scenarios: (a) D1; (b) D1 with added mass; (c) D2; (d) D2 with added mass. The concrete block with a weight of 0.545 kg (the same as Section 4.2.1) is added to each floor to construct the modified structure for both damage scenarios. The same measurement procedures were performed as described in Section 4.2 for two damage scenarios. Fig. 13 shows the experimental setup of two damage scenarios. Modal analysis is also implemented by automated SSI to extract natural frequencies and mode shapes of the original and modified system in two damage scenarios. The proposed method is then performed to identify MCs and SCs based on the updated FE model (healthy condition), as described 19 in Section 4.2.1. In the updating process, the DREAM algorithm generates samples to target the posterior PDF by the same initial settings as before. Figs. 14 -17 show the updated results of mass and stiffness in two damage scenarios. Figs. (a)s in Fig. 14 - 17 are trace plots that show the iteration of each coefficient with samples increasing; stable convergence is achieved in trace plots. The updating results of MSc and SCs over ten Markov chains are presented in Figs. (b)s in Fig. 14 - 17. It is seen that all coefficients in both damage scenarios are identified as consistent with each other among ten Markov chains, indicating reliable and accurate identification. In addition, the convergence diagnosis, 𝑅.=8= quickly decays below 1.2 for all coefficients, demonstrating that the stationary convergence is reached. In damage detection, the last 10,000 samples are used to calculate the mean values and standard derivations. Note the mean values under the healthy condition in Section 4.2.1 are used as baselines, so the undamaged floor has the MC and SC with unity value. (a) (b) (c) Fig. 14. Results of updated mass in D1: (a) trace plot; (b) square: the sample mean of each chain, error bar: ±2∙S.D.; (c) convergence diagnosis, 𝑅.=8= (a) (b) (c) 20 Fig. 15. Results of updated stiffness in D1: (a) trace plot; (b) square: the sample mean of each chain, error bar: ±2∙S.D.; (c) convergence diagnosis, 𝑅.=8= (a) (b) (c) Fig. 16. Results of updated mass in D2: (a) trace plot; (b) square: the sample mean of each chain, error bar: ±2standard derivations; (c) convergence diagnosis, 𝑅.=8= (a) (b) Fig. 17. Results of updated stiffness in D2: (a) trace plot; (b) square: the sample mean of each chain, error bar: ±2standard derivations; (c) convergence diagnosis, 𝑅.=8= 21 (a) (b) Fig. 18. Density distribution of updated coefficients: (a) D1; (b) D2 Fig. 18 shows the density distribution estimated by GKE. It is observed that some coefficients exhibit multi-modal features and are non-Gaussian shaped. Especially, 𝜃& and 𝜃< in damage scenario No. 1 (D1) and 𝛽&, 𝛽( and 𝜃< in damage scenario No. 2 (D2). It indicates that structural parameters do not always follow Gaussian distribution, the asymptotic optimization method may not suitable to estimate the posterior PDF with a non-Gaussian shape. However, the proposed method is able to approximate non-Gaussian distribution with an accurate level. It is also found that the larger uncertainties are revealed in some coefficients, such as 𝛽( and 𝜃& in both damage scenarios (S.D. ranges from 7.9% to 16.2%). Their distributions are flatter and spread across a relatively wider region. While the distributions of other coefficients are concentrated in a narrow region and have pronuounced peaks, meaining these coefficients are more certain (S.D. ranges from 1% to 10.1%). The identified damage severities of two damage scenarios is shown in Fig. 19. The identified damage severities of 𝛽< and 𝜃< in damage scenario No. 1 (D1) are 23.35% and 24.72%, respectively, which is close to actual values of 21.5% for 𝛽 and 21.8% for 𝜃; The identified damage severities of 𝛽(, 𝛽<, and 𝜃(, 𝜃< in damage scenario No. 2 (D2) are 24.55%, 23.14%, 19.18% and 20.18% respectively, which also agree well with actual values of 21.5% for 𝛽 and 21.8% for 𝜃. The false damage detection is only observed with less than 4%. These results demonstrate the efficacy of the proposed method in both mass and stiffness updatings and achieves damage localization and quantification. (a) (b) Fig. 19. Identified damage severity: (a) D1; (b): D2 Table 9. Results of updated frequencies and MAC values in two damage scenarios Damage scenario Mode No. Measured Frequency (Hz) MAC Updated Error (%) D1 1 7.140 7.234 1.31 0.9942 2 19.800 19.420 1.92 0.9936 3 32.311 32.680 1.14 0.9976 D2 1 6.905 6.922 0.25 0.9941 2 19.641 20.244 2.98 0.9957 3 29.043 29.512 1.59 0.9999 22 (a) (b) Fig. 20. Probabilistic damage curves: a: D1; b: D2 The updated frequencies and MAC values are calculated using identified MCs and SCs, as shown in Table 9. It is observed that all modal parameters in both damage scenarios are in accordance with measured counterparts, indicating the FE model is successfully updated by the proposed method. Based on idetified mean values of MCs and SCs and their uncertainties under healthy and damaged state, the probabilitisc damage curves can be plotted, as shown in Fig. 20. It is worth mentioning that the negative, 𝑑 represents mass/stiffness increase, and vice versa. For damage scenario No.1 (D1), 𝛽< and 𝜃< have a probability (63.3% and 67.5%, respectively) of having possible damage 21.5% and 21.8%. For damage scenario No. 2 (D2), 𝛽( and 𝜃( have a possible change of 21.5% and 21.8% with a probability of 62.1% and 67.7%, respectively. In practice, damage can be detected by probabilistic curves, because the curves related to damage location are generally easily distinguished from the ones related to healthy location. For example, the curves of 𝛽< and 𝜃< in D1 are clearly separated from others, indicating the location corresponding to 𝛽< and 𝜃< (herein is the third floor) may have certain damage. Similar observation is found in D2. The proposed method can detect damage in mass and stiffness along with location and severity in a probabilistic manner. The engineers can be informed that some repairing work may be necessary at certain location. 5. Conclusions This paper proposed a novel vibration-based Bayesian model updating approach to simultaneously identify structural mass and stiffness. In this work, the coupling effect of mass and stiffness is successfully addressed using two sets of vibration data from original and modified system with added known mass. Following conclusion and contributions from numerical examples and experimental tests are summarized as follows: • The results in numerical example and experimental test illustrate that the proposed approach can simultaneously identify structural mass and stiffness with an accurate level and their uncertainties by addressing the coupling effect of mass and stiffness. • In experimental test, some mass coefficients exhibit larger uncertainties, indicating the effect of mass on structural integrity cannot be ignored, and the assumption of mass is known and invariant in classical Bayesian approach may be questionable when noticeable change in mass is observed, such as 21.5% mass increase in damage scenarios in this test for mimicking the mass change due to unknown damages. • The results in experimental test reveal the structural parameters, e.g., mass and stiffness, do not always follow Gaussian distribution. Thus, the asymptotic approximation method may not be suitable for this situation. The DREAM algorithm runs multiple Markov chains in parallel and sufficiently seek all possible solutions, resulting in high capability to treat the posterior PDF with high-dimensionality, multi-modality, and numerous peaks. To the best of authors’ knowledge, it is the first attempt to apply DREAM for Bayesian model updating in the field of civil engineering. • The probabilistic damage detection is also implemented by the proposed Bayesian approach. The results in experimental test demonstrate that the proposed approach enables to reliably and accurately identify damage location and severity. In addition, the probabilistic damage curves allow engineers to quickly localize damage, indicating the proposed approach is practically valuable. 23 Although the proposed vibration-based Bayesian approach exhibits potential applications, the efficacy of the approach should be verified in more complex structures. In addition, in the current study, the stationary mass is added to the structure, which may not be always practical in real-world settings. Therefore, moving mass, e.g., vehicles on bridge or elevator in building, can be used to create modified systems that can be considered in the future study. CRediT authorship contribution statement Jice Zeng: Methodology, Software, Validation, Formal analysis, Writing-original draft & editing; Young Hoon Kim: Conceptualization, Supervision, Review & editing. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement The authors would like to express gratitude and sincere appreciation for the partial financial support by the University of Louisville.