On choice and effect of weight matrix for response sensitivity-based damage identification with measurement and model errors

Mechanical Systems and Signal Processing 114 (2019) 1–24

Contents lists available at ScienceDirect

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

On choice and effect of weight matrix for response sensitivity-based damage identification with measurement and model errors Zhong-Rong Lu, Junxian Zhou, Li Wang ⇑ Department of Applied Mechanics and Engineering, Sun Yat-sen University, Guangzhou, PR China

a r t i c l e

i n f o

Article history: Received 3 June 2017 Received in revised form 17 April 2018 Accepted 1 May 2018

Keywords: Optimal weight matrix Structural damage identification Measurement and model errors Enhanced response sensitivity approach Hybrid data

a b s t r a c t This paper aims to present a thorough view on the choice and effect of the weight matrix for response sensitivity-based damage identification with measurement and/or model errors. The derivation of the optimal weight matrix is mainly twofold. On the one hand, when only measurement errors are involved, the optimal weight matrix is found to be inverse proportional to the measurement error covariance by minimizing the expectation of squares error of the whole identification results. On the other hand, if model errors are additionally considered, the optimal weight matrix then depends not only on the measurement error covariance, but also on the model error covariance. Further analysis reveals that the optimal weight matrix can also make the ‘relative error’—square-root of expectation of squares error in every individual damage parameter minimized. Then, the effect of the proposed optimal weight matrix with measurement and/or model errors is studied on two typical examples—a plane frame and a simply-supported plate. Results show that when hybrid types of measurement data—accelerations, displacements and/or eigenfrequencies are used or when the response data is sensitive to model errors, the optimal weight matrix should be invoked to get reasonably good identification results and the improvements brought by the optimal weight matrix are substantial. The whole work shall be instructive for damage identification when different types of measurements are available and when model errors are non-negligible. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction In-service structures are always working under environmental aggressions, possibly unexpected loads and material creep etc, and as a consequence, damage or degradation of structural stiffness may occur which is harmful for practical use. To avoid this, the structural health monitoring (SHM) system [1] is necessarily introduced for the sake of preventive monitoring, detection and maintenance. Indeed, the central ingredient of the SHM system resides in the structural damage identification. Generally, damage is embodied as the reduction of structural stiffness and therefore, the main task of structural damage identification is to inversely identify the damage locations and intensities (or stiffness reductions) from the measured data, e.g., acceleration responses, displacement responses, eigenfrequencies and/or eigen modes of a structure. Typically, such a procedure belongs to a class of inverse problems and is often modeled as a nonlinear (weighted) least-squares optimization ⇑ Corresponding author. E-mail address: [email protected] (L. Wang). https://doi.org/10.1016/j.ymssp.2018.05.007 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

2

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

problem [2,3]. Various methods including the genetic-like algorithms [4–6] and the gradient-based approaches [2,3,7] have been developed to solve such a nonlinear least-squares optimization problem, among which the hessian-free gradient-based minimization techniques such as the Gauss-Newton method, the Levenberg-Marquardt method [19] and the sensitivitybased approaches [4] are often preferred. In this paper, the sensitivity-based approaches [3,4,7–11] are used to solve the nonlinear least-squares optimization problem under the following considerations,
Structural damage identification (or finite element model updating) is often accomplished with a large number of unknowns and therefore, the genetic-like algorithms with random search nature [5,6] are often not suitable for this problem due to the prohibitively high computational cost;
Among the gradient-based approaches, the commonly used Newton method would require more complex second-order sensitivity analysis and more computational efforts [3] than the sensitivity-based approaches for which only the (firstorder) sensitivity analysis should be conducted; this indicates that the Newton method is also inappropriate for the optimization problem;
Recently, the sensitivity-based approaches along with proper regularizations or equivalent trust-region constraints have been shown to be weakly convergent [3]. In principle, the success of the sensitivity-based approaches lies in the sensitivity analysis and usage of regularization techniques [3,14]. On the one hand, different types of the measured data often correspond to different models [12] of the structure; for instance, time-domain data corresponds to the dynamic equation while eigenfrequencies and eigen modes are derived from modal analysis. On the other hand, different models would lead to quite different sensitivity analysis [2] and therefore, the sensitivity-based approaches should be adjusted with respect to the measurement types. Over the years, considerable work has been focused on the development and application of the sensitivity-based approaches for damage identification in various kinds of structures and with various types of the measured data. Classically, the modal data including eigenfrequencies and eigen modes is most used for damage identification. Farhat and Hemez [7] proposed to use eigen modes and the sensitivity analysis was found to be conducted in an element-by-element manner. Later, Bakir et al. [9] well identified the stiffness parameters from the modal data by incorporating the trust-region algorithm into the sensitivity approach. As is noteworthy, complete data of an arbitrary mode is required for the above two approaches which is often unavailable in many test cases. To this end, Chen and Maung [8] developed a regularized sensitivity approach for finite element model updating with incomplete modal data. In practice, the modal data including eigenfrequencies and incomplete modes is easily accessed and this explains why substantial attention was paid to damage identification or finite element model updating with the modal data in the literature. However, high-frequency information is often lost when experimental modal data is utilized for damage identification and this means that the modal data is not very suited for detecting and quantifying localized small size damage [12]. In contrast, as pointed out by Link and Weiland [12], time domain response data from impact tests carries high-frequency information and is sensitive to damage even if the damage is local [15]; this indicates that time domain response data is beneficial for the detection of local damage. Nevertheless, few publications [3,10– 12] were found until recently addressing time domain response data in conjunction with the sensitivity analysis for damage identification. Lu and Law [10] proposed to use response sensitivity analysis for inverse damage identification, and later, Lu and Wang [3] enhanced the response sensitivity approach so that even large damage can be identified; Link and Weiland [12] conducted an elaborate comparative study for sensitivity-based damage identification using modal data and using time domain response data, for which the latter was shown to be superior in damage identification. In this paper, the time domain response sensitivity approach [3] for damage identification is followed up for further study. The eigenfrequency data which is easily available may be additionally used to explore the possible improvement in damage identification with hybrid data [13]. As is mentioned above, the objective function of damage identification (or the optimization problem) is in the weighted least-squares form for which the residual between the measured data and the implicitly derived data and a weight matrix are involved. The sensitivity analysis is conducted regarding the residual term, while limited attention has been paid to the choice and effect of the weight matrix for the weighted nonlinear least-squares optimization problems. Notwithstanding, the weight matrix shall be carefully selected. On the one hand, different kinds of data—displacement, velocity, acceleration or even eigenfrequency correspond to different dimensions, e.g., m/s2 for acceleration and m for displacement where m represents length in ‘meter’ and s time in ‘second’. Without properly defined weight matrix, equivalent replacement of dimensions, e.g., m replaced by 103 mm (or millimeter) may lead to quite different objective functions and different identification results thereof. To avoid this, the scaling effect behind the weight matrix should be implicated. On the other hand, when the measured data is at different levels of errors, the weight matrix should be chosen such that large error corresponds to less weight. Primitively, there have already been some reasonable choices of the weight matrix [2,16–18], such as being the reciprocal of the covariance matrix of the measured data. Particularly, in sensitivity-based model updating, proper choice of the weight matrix in each updating step may even pose the equivalent effect of regularization [2,17]. In this paper, the main aim is to present a thorough view on the optimal choice and the specific effect of the weight matrix for response-sensitivity-based damage identification. Distinct to the work in [2,17], the weight matrix herein for the nonlinear least-squares optimization problem does not change during the iterative identification procedure and does not intend to pose the regularization effect; rather, it tries to render the expectation of the squares error of the final identified results minimized. Moreover, in addition to the measurement errors, the model errors in the external load, the damping and even the predefined stiffness of the joint are also considered and how to properly choose the weight matrix in such cases is also addressed.

3

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

The remainder of the paper is structured as follows. Section 2 briefly states the damage identification procedure as a nonlinear least-squares optimization problem. In Section 3, the idea behind the enhanced response sensitivity approach [3] is followed up and generalized to cope with time-domain response and eigenfrequency data. How to choose the optimal weight matrix with measurement and/or model errors is specifically analyzed in Section 4. Moreover, how to correct the model parameters with errors is also presented. In Section 5, two numerical examples are studied in detail to see the effect of the proposed optimal weight matrix and final conclusions are drawn in Section 6. 2. Problem statement Consider a structure with n degrees-of-freedom (DOFs) governed by the following dynamic equations,

€ ðtÞ þ CuðtÞ _ Mu þ KuðtÞ ¼ f ðtÞ; _ ¼ u_ 0 uð0Þ ¼ u0 ; uð0Þ

t > 0;

ð1Þ

where M; C; K are the time-invariant mass, damping and stiffness matrices of a discrete-mass structure or of a continuous _ u € represent the displacement, velocity and acceleration structure after finite element discretization. u ¼ ½u1 ; u2 ; . . . ; un T ; u; responses of the structure, respectively and u0 ; u_ 0 are initial displacements and velocities. f denotes the external load vector. Usually, the stiffness matrix is represented as an affine function of the stiffness or damage parameters a 2 A # Rm where A P denotes the admissible parametric space, i.e., K ¼ KðaÞ :¼ K0 þ m k¼1 ak Kk . The damping matrix C is assumed to be the Rayleigh damping, that is, C ¼ a1 M þ a2 K with a1 ; a2 being two positive coefficients. In this way, C is also an affine function of the damage parameters a. Remark 1. The Rayleigh damping used in this paper is a simple as well as widely-used model to describe the damping in structures. Such a damping arises naturally after finite element modeling of the visco-elastic structural components where the two Rayleigh coefficients a1 ; a2 are used to describe the linear viscosity of the material. Generally, the structural stiffness would degrade after damage and then, by the Rayleigh damping definition C ¼ a1 M þ a2 K, the Rayleigh damping would also decrease after damage; this is the practice only in some limited cases such as when the damage is embodied as a portion reduction of a certain structural section. In most damage cases such as when plasticity or cracks appear somewhere, the damping would increase with the appearance of the damage and under these circumstances, new damping models shall be introduced. In this paper, the Rayleigh damping is used for simplicity, nevertheless, other damping models can be analyzed and tackled likewise. Based on the Rayleigh damping assumption, the damage parameters a would arise linearly in the stiffness matrix K and the damping matrix C, while for the model parameters b with possible errors to be discussed later in this paper, they may be implicated in the mass matrix M, the damping matrix C, the stiffness matrix K and even the load vector f . It is noteworthy that the model error to be considered in this paper is restricted to be caused by the aleatoric uncertainty in the specific model parameters, rather than the epistemic uncertainty of the used model. The dynamic Eqs. (1) can be numerically solved by the Newmark-b method. In addition, with the mass and stiffness matrices, one can also get the eigenfrequencies and eigen modes ðxk ; /k Þ; k ¼ 1; 2; . . .with x1 6 x2 6 . . .such that

K/k ¼ x2k M/k :

ð2Þ

In damage identification, the main goal is to identify the damage parameters a from the measured data; the data includes but is not limit to the acceleration (resp. displacement) responses of the kth DOF at a series of l sampling time points ^ € ¼ ½u ^ ðt 1 Þ; u ^ ðt2 Þ; . . . ; u ^ ðt ÞT (resp. d ^ ¼ ½u € € € ^ ðt1 Þ; u ^ ðt2 Þ; . . . ; u ^ ðt ÞT ) and the first 0 < t1 < t2 < . . . < t that are represented as d l

k

k

k

k

l

k

k

k

k

l

^ ¼ ½x ^ 21 ; x ^ 22 ; . . . ; x ^ 2q T in the increasing order where the over hat symbol ^ means the measured data of the corq eigenvalues K ^ and a possible form of R ^ shall be responding quantity. For brevity, all the measured data is stored in a column vector R

1 ^ K C B^ ^ ¼ B di ; i 2 S u C R A @ ^€ d € j ; j 2 Su 0

ð3Þ

where S u ; S u€ are two sets of DOFs where respective displacement and acceleration responses are measured. Theoretically, given the damage parameters a, the results of the measured quantities corresponding to the measured data can be acquired by resorting to the dynamic Eqs. (1) and the modal analysis (2) and therefore, they are denoted by RðaÞ. In other words, given the damage parameters a and the two models (1) and (2), the measured quantities RðaÞ can be implicitly computed as

0

KðaÞ

1

0

½x21 ; x22 ; . . . ; x2q 

T

1

B C B C T C RðaÞ ¼ @ di ðaÞ; i 2 S u A :¼ B @ ½ui ðt1 Þ; ui ðt 2 Þ; . . . ; ui ðt l Þ ; i 2 S u A: € T dj ðaÞ; j 2 S u€ € j ðt1 Þ; u € j ðt 2 Þ; . . . ; u € j ðt l Þ ; j 2 S u€ ½u where the first q eigenvalues x21 ; x22 ; . . . ; x2q are also in the increasing order.

ð4Þ

4

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

^ defined and obtained, a general way to tackle the damage identification is to find the damage With both RðaÞ and R parameters a that solve the following weighted nonlinear least-squares optimization problem,

n

a ¼ arg min gðaÞ :¼ kR^  RðaÞk2W

o

ð5Þ

a2A

where W is a positive-definite weight matrix and k  kW :¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðÞT WðÞ. Different weight matrix W could result in different

damage identification results and how to properly choose such a weight matrix will be addressed in Section 4. In what follows, how to solve the damage identification problem (5) with the measured eigenvalue and response data (3) via the sensitivity analysis [3] is presented. 3. Damage identification based on enhanced sensitivity analysis 3.1. Sensitivity analysis Problem (5) is a typical nonlinear optimization problem and therefore, it should be solved iteratively. To this end, the key  such that gða  þ daÞ lies in how to quickly get a proper update da from the previously obtained or initialized parameters a becomes as small as possible. For such a nonlinear least-squares problem (5), a general way to handle the nonlinearity [19] is ^  Rða  þ daÞ as follows, to linearize R

^  Rða  þ daÞ  dRða  Þ  Sða  Þda R ^  Þ :¼ R  Rða  Þ; dRða 0h

@K @ a1

@K ; @@K a2 ; . . . ; @ am

i

1

C Bh i C B @d @d @di i i C  Þ ¼ ra Rða  Þ :¼ B Sða B @ a1 ; @ a2 ; . . . ; @ am ; i 2 S u C A @ € € € @d @d @d ½@ a1j ; @ a2j ; . . . ; @amj ; j 2 S u€

ð6Þ

^ and the derived data Rða  Þ is called the residual between the measured data R  Þ, and Sða  Þ the (first-order) senwhere dRða sitivity matrix. More specifically, the residual term dRðaÞ can be easily obtained by directly solving the forward problems (1) and (2), while the sensitivity matrix SðaÞ requires sensitivity analysis on the dynamic Eqs. (1) and the modal Eqs. (2) for which the detailed procedure reads,
Response

sensitivity

analysis.

Sensitivity of the displacement and acceleration T € € € € di ¼ ½ui ðt 1 Þ; ui ðt 2 Þ; . . . ; ui ðt l Þ ; i 2 S u ; dj ¼ ½uj ðt 1 Þ; uj ðt2 Þ; . . . ; uj ðtl Þ ; j 2 S u€ to the damage parameters

responses

T

0 @u ðt Þ 1 i

0 @ u€j ðt1 Þ 1

1

@ ak

B @u ðt Þ C B i 2 C @ ak C @di B C; ¼B C @ ak B B ... C A @ @ui ðt l Þ @ ak

@a

B € k C B @ uj ðt2 Þ C € @ dj B @ ak C C; k 2 f1; 2; . . . ; mg ¼B C @ ak B B .. C . A @ € j ðt l Þ @u @ ak

should pertain to the following dynamic sensitivity equations [10],

€ @u @ u_ @u _ ðtÞ þ C ðtÞ þ K ðtÞ ¼ Kk ðuðtÞ þ a2 uðtÞÞ; @ ak @ ak @ ak @u @ u_ ð0Þ ¼ 0; ð0Þ ¼ 0 @ ak @ ak

M

t > 0;

ð7Þ

_ which are in the same form with the forward problem (1) except that the load vector f ðtÞ is replaced by Kk ðuðtÞ þ a2 uðtÞÞ and the initial states are set to zeros herein.
Eigenvalue sensitivity analysis. Sensitivity of the eigenvalues to the damage parameters

0 @ðx2 Þ 1 1

B @ ak C B @ðx2 Þ C 2 C @K B B @a C ¼ B k C; @ ak B .. C B . C A @

k 2 f1; 2; . . . ; mg

@ðx2q Þ @ ak

is performed as follows [2],

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

@ðx2j Þ @ ak

¼

/Tj



@K @ ak

  x2j @@M a /j k

/Tj M/j

¼

/Tj Kk /j /Tj M/j

;

j ¼ 1; 2; . . . ; q

5

ð8Þ

where /j ; xj are the eigen pair obtained by solving the modal Eqs. (2). Now, with the linearization (6), a linear least-squares approximation of the original goal function gðaÞ, designated as  ; daÞ is obtained, i.e., g~ða

 Þ  Sða  Þdak2W  ; daÞ ¼ kdRða g~ða

ð9Þ

and then, the update da is determined regarding the easily solvable linear least-squares problem (9); this will be elaborated in the next subsection. 3.2. Enhanced response sensitivity approach The enhanced response sensitivity approach operates on Eq. (9) to get the proper update and thereafter, to iteratively solve the nonlinear least-squares optimization problem. In doing so, two aspects are considered in order to acquire a proper update da. On the one hand, the linear least-squares problem (9) can be easily solved, however, it may be ill-posed. To circumvent the possible ill-posedness, the Tikhonov regularization [20] shall be introduced, yielding the following update 1

 Þ  Sða  Þdak2W þ kkdak2 ¼ ½ST ða  ÞWSða  Þ þ kI ST ða  ÞWdRða Þ dak ¼ arg minkdRða da

ð10Þ

where k P 0 is the regularization parameter, k  k denotes the usual L2 norm of vectors and I is the identity matrix.  þ daÞ into the On the other hand, the key for transforming the original nonlinear least-squares objective function gða  ; daÞ agree well with  ; daÞ lies in the linearization (6). Thus, to make g ~ ða approximate linear least-squares form g~ða  ; daÞ agree with  þ daÞ, the update da should be reasonably small. To indicate how well the approximate goal function g ~ ða gða  þ daÞ, an agreement indicator qðda; a  Þ is defined, the original nonlinear function gða

 Þ  gða  þ daÞ  Þk2W  kdRða  þ daÞk2W gða kdRða qðda; a Þ ¼ ~  : ¼ 2  ; daÞ kdRða g ða; 0Þ  g~ða  ÞkW  kdRða  Þ  Sða  Þdak2W

ð11Þ

Good agreement requires that the agreement indicator should be larger than a critical number qcr [19,3], i.e., the following agreement condition should be satisfied,

qðda; a Þ P qcr 2 ½0:25; 0:75:

ð12Þ

Finding a proper update that verifies the agreement condition is just called the trust-region constraint. An implicit perspec ; 0Þ along with the  ; daÞ < g~ða tive behind the trust-region constraint is that descent of the approximate goal function g~ða  þ daÞ < gða  Þ. agreement condition (12) means descent of the original nonlinear goal function gða At the first glance, the above two aspects—the Tikhonov regularization and the trust-region constraint seem quite different, nevertheless, choice of a reasonable regularization parameter k and treatment of the trust-region constraint can be considered at the same time to enhance the sensitivity-based approach. In fact, the Tikhonov regularization is to some extent  is not a minima of gðaÞ or equivalent to the trust-region constraint. The equivalence lies in the following fact: if a T    kra gðaÞk ¼ 2kS ðaÞWdRðaÞk – 0, there would hold,

lim kdak k ¼ 0;

k!þ1

 Þ ¼ 1 > qcr lim qðdak ; a

k!þ1

ð13Þ

where dak is obtained via the Tikhonov regularization (10). Proof of the equalities in (13) can be found in reference [3]. It is implicated in Eq. (13) that there always exists a large regularization parameter k such that the update dak is small and the agreement condition (12) is satisfied. To this end, one can follow a recursive procedure to get a proper regularization parameter k and a proper update dak thereof; that is, 1. Get the initial regularization parameter k. Usually, the initial parameter is obtained from the L-curve method [21,22]  Þ. In prinregarding only the linearized problem (10) and to be distinct, such a regularization parameter is denoted by kL ða  Þ  Sða  Þdak kW and the update term kdak k, and the ciple, the L-curve method tries to balance the residual term kdRða  Þ is attained at the corner of the curve plotted by the residual term and the update term. For specific calresulting kL ða  Þ as the initial value is that such culation tools on the L-curve method, refer to the work in [22]. The reason for choose kL ða a choice as the final regularization parameter has already been shown to give good identification results for small damages [10,11]. 2. Compute the update dak as in (10).  Þ as in (11). 3. Calculate the agreement indicator qðdak ; a

6

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

4. If the agreement condition (12) is satisfied, terminate the recursion; otherwise, increase the regularization parameter k up to a factor c > 1, i.e., k :¼ ck and return to step 2. Eventually, with the proper update dak obtained, an iterative algorithm to solve the nonlinear least-squares optimization problem (5) can be established and details are presented in Table 1. A remarkable property of the enhanced response sensitivity approach is that with the trust-region constraint, the algorithm has been proved to be weakly convergent [3]; that is,

lim kra gðaðkÞ Þk ¼ 0:

ð14Þ

k!þ1

4. Optimal weight matrix under measurement and model errors 4.1. Optimal weight matrix under measurement errors In the previous section, an enhanced response sensitivity approach has been presented to identify the damage with userdefined weight matrix W and from the measured data irrespective of the measurement errors. Notwithstanding, a good weight matrix W remains to be determined and this is to be fulfilled in this section. To begin with, assume that the measured data is of the following form,

^ ¼ Rðaex Þ þ R R

ð15Þ

where a represent the exact damage parameters and R is the measured error vector pertaining to the Gaussian distribution subject to zero means and positive definite covariance E½R TR  ¼ Q R with E½ denoting the expectation operator. To have an impression on how to calibrate the measurement errors, one can refer to the work [23] for calibration of accelerometers for the acceleration response data and the literature [24] for assessment of displacement measurement errors. Without measurement errors R ¼ 0, the algorithm in Table 1 is expected to give almost exact identification of the damage parameters (see numerical examples in reference [3]); however, errors that are often small would always exist in practice R – 0 and on considering this, the identified results a obtained by the enhanced response sensitivity approach would always admit some small perturbations, i.e., a – aex or a  aex are often of small errors. Different weight matrices W would lead to different identification results a and thereof different identification accuracy a  aex . An optimal weight matrix shall render the expectation of squares of the identification errors a  aex minimized; that is to say, optimal choice of the weight matrix shall be modeled as: ex

Find : positive definite W;

ð16Þ

Minimize : E½ða  aex Þ ða  aex Þ: T

To solve (16), the relationship between W and E½ða  aex ÞT ða  aex Þ is essentially derived. To go further, the goal function gðaÞ (5) is further analyzed by linearization,

gðaÞ ¼ kRðaex Þ þ R  RðaÞk2W ¼ kR  S  ða  aex Þ þ h:o:tk2W

ð17Þ

where S :¼ Sðaex Þ and h:o:t represents higher order terms that are negligible with respect to R when a is in the vicinity of aex . Note that a is in the vicinity of aex and will also be the minima for (17) among all a in the vicinity of aex ; this means, 1

a  aex  ½ST WS ST WR :¼ LW R

ð18Þ

1 T

where LW ¼ ½ST WS S W. Then, with Eq. (18), the following would hold

E½ða  aex Þ ða  aex Þ  E½TR LTW LW R  ¼ E½trðTR LTW LW R Þ ¼ E½trðLW R TR LTW Þ ¼ trðLW Q R LTW Þ T

where trðÞ is the trace operator and the last equality arises due to E½  (16) approximately with Eq. (19) is presented.

T R R

ð19Þ

¼ Q R . Next, how to solve the optimization problem

Theorem 1. The optimal choice of W to minimize FðWÞ :¼ trðLW Q R LTW Þ is

Wopt ¼ cQ 1 R ; c > 0 a constant

ð20Þ

and then, 1

FðWopt Þ ¼ trð½ST Q 1 R S Þ:

ð21Þ

Proof. Herein, for convenience, manipulations will be put on the auxiliary matrix LW (see Eq. (18)) rather than directly on W. A remarkable feature of LW is

LW S ¼ I

ð22Þ

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

7

Table 1 Algorithmic details for enhanced response sensitivity approach. - set initial damage parameters að1Þ 2 A as those of the intact structure, define error tolerance RelErr (e.g., ¼ 106 ) for convergence criterion, fix the maximum number of iterations Nmax (e.g., =1000), fix trust-region parameters qcr 2 ½0:25; 075 (e.g., =0.5) and c > 1 (e.g., =2), set the maximum number of steps for trust-region procedure Ntr, (e.g.,=20), ^ and user-defined weight matrix W, - load the measured data R - for k ¼ 1 : Nmax

-

- solve (1) and (2) to get RðaðkÞ Þ,

^  RðaðkÞ Þ, - compute the residual dR ¼ R

- do sensitivity analysis (7) and (8) to get the sensitivity matrix SðaðkÞ Þ,

- use L-curve method to get the initial regularization parameter k ¼ kL ðaðkÞ Þ, - for i ¼ 1 : Ntr % Trust-region procedure - k ¼ ci1 kL ðaðkÞ Þ, 1 T

- compute the update da ¼ ½ST ðaðkÞ ÞWSðaðkÞ Þ þ kI - if a

ðkÞ

þ da R A

S ðaðkÞ ÞWdR,

continue,

- solve (1) and (2) to get RðaðkÞ þ daÞ, ^  RðaðkÞ þ daÞ, - compute new residual dRnew ¼ R - calculate the agreement indicator q ¼ - if q P qcr - end for

kdRk2 kdRnew k2 , kdRk2 kdRSðaðkÞ Þdak2

break,

- update damage parameters aðkþ1Þ ¼ aðkÞ þ da,

- if kdak=kaðkþ1Þ k 6 RelErr - end for

break.

and this forms a critical constraint in minimizing trðLW Q R LTW Þ. To this end, the Lagrangian multipliers H are introduced to easily enforce the constraints (22),

LðLW ; HÞ ¼ trðLW Q R LTW Þ þ trððI  LW SÞHÞ:

ð23Þ

where LðLW ; HÞ is also called the Lagrange function. Optimization/minimization of LðLW ; HÞ over LW leads to,

2Q R LTW ¼ SH ) LW ¼

1 T T 1 H S QR : 2

ð24Þ

Then, combination of (24) and (22) yields, 1

H ¼ 2ðST Q 1 R SÞ :

ð25Þ

Along with (24), optimal LW is found to be 1

T 1 T 1 Lopt W ¼ ðS Q R SÞ S Q R :

ð26Þ

On considering the definitions of LW and FðWÞ, optimal choice of the weight matrix in (20) as well as the minimum of FðWÞ in (21) follows naturally. h Theorem 1 has presented an optimal choice of the weight matrix in the sense that (approximate) expectation of squares error of identification results is minimized. The optimal weight matrix Wopt ¼ cQ 1 R as proportional to the reciprocal of the measurement error covariance is independent of the model problem as well as the exact damage parameters, rather, it is merely related to the measurement errors; such an optimal choice of the weight matrix also coincides with the results in linear estimation problems [2,18] as well as those undesignedly used in the literature [16,17]. For illustrative purpose, assume in the measured data (3) that the errors in displacement (resp. acceleration) responses at different time points and different positions of a structure are independent and have the same variance r2u (resp. r2u€ ) while the errors in different eigenvalues are also independent and pertain to the same variance r2K . Under this circumstance, the optimal weight matrix is 2 2 a diagonal matrix with diagonals being r2 € K in correspondence to the eigenvalue data, ru to the displacement data and ru to the acceleration data. In this way, the goal function gðaÞ under the optimal weight matrix (c ¼ 1) becomes

2  2 2     ^€  K ^  d ðaÞ € X Xd    ^  KðaÞ d i i j  dj ðaÞ 2 ^  gðaÞ ¼ kR  RðaÞkWopt ¼   þ  þ   :    rK   r r € u u  i2S j2S  u

ð27Þ

€ u

On the one hand, note that rK ; ru ; ru€ are of the same dimensions with eigenvalues, displacements and accelerations, respec^

^

aÞ di di ðaÞ tively and therefore, KKð rK ; ru ;

^ € ðaÞ € d d j j

ru€

are dimensionless; in other words, a scaling effect has been reached by the optimal

8

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

weight matrix. On the other hand, the optimal weight is proportional to the reciprocal of the covariance of measurement errors, implying that larger measurement errors correspond to less weights; this is intuitively reasonable. Actually, in addition to minimal expectation of squares error of the whole identification results, the optimal weight matrix (20) can also exhibit some other attractive features; these are to be discussed below. Theorem 2. Introduce an order operator  (or inversely ) for two symmetric matrix A; B 2 RN N of the same size as follows,

A  Bðor B AÞ () v T Av 6 v T Bv ;

8v 2 RN :

ð28Þ

Then, covariance of the identification errors E½ða  aex Þða  aex ÞT  ¼ E½LW R TR LTW  ¼ LW Q R LTW would also arrive at a minimum under the optimal weight matrix Wopt ¼ cQ 1 R in the sense of the order operator , i.e., T

T 1 opt Lopt W Q R ðLW Þ ¼ ½S Q R S

1

 LW Q R LTW ; 8 positive definite W:

ð29Þ

Proof. For an arbitrary positive definite weight matrix W, one can have

LW ¼ Lopt W þ DL

ð30Þ

with DL a matrix satisfying DL  S ¼ 0. Then, recalling the expression of Lopt W in (26), it is also found that 1

1

T T 1 T 1 T 1 T T Lopt W Q R DL ¼ ðS Q R SÞ S Q R Q R DL ¼ ðS Q R SÞ ðDL  SÞ ¼ 0:

ð31Þ

With Eqs. (30) and (31), there holds T

T

opt opt opt opt opt T T T LW Q R LTW ¼ ðLopt W þ DLÞQ R ðLW þ DLÞ ¼ LW Q R ðLW Þ þ DLQ R DL þ LW Q R DL þ ðLW Q R DL Þ T

T

T 1 opt opt opt T ¼ Lopt W Q R ðLW Þ þ DLQ R DL LW Q R ðLW Þ ¼ ½S Q R S

1

T

ð32Þ

T

where the inequality arises due to the fact that DLQ R DL is nonnegative definite. h Remark 2. From the estimation theory [18], the minimum covariance of the damage parameter estimates is given by the Cramér-Rao lower bound which is the inverse of the Fisher information matrix ST Q 1 R S. Then, it is found from Theorem 2 that the Cramér-Rao lower bound is reached if and only if the optimal weight matrix Wopt is used in the damage identification. In fact, one can find that the results in Theorem 1 are special cases of the results in Theorem 2 because the following holds 1

½ST Q 1 R S

1

T  LW Q R LTW ) trð½ST Q 1 R S Þ 6 trðLW Q R LW Þ:

ð33Þ

Moreover, let ej 2 R ; j ¼ 1; 2; . . . ; m denote a column vector with the jth element being 1 and others zeros and then, considering (28) and (29), there is m

T

T opt T  ex  ex T T opt E½ðaj  aex j Þ  ¼ ej E½ða  a Þða  a Þ ej ¼ ej LW Q R LW ej P ej LW Q R ðLW Þ ej 2

T

ð34Þ

implying that expectation of the identification squares error in an arbitrary/ individual damage parameter is also minimized under the optimal weight matrix (20). Furthermore, replacement of ej by ja1ex j ej in (34) would also show that the optimal j   ex 2 jaj aj j in every individual damweight matrix (20) would again lead to minimal expectation of the squares relative error jaex j j

age parameter. Above all, the optimal choice of the weight matrix Wopt ¼ cQ 1 R under measurement errors not only makes the expectation of squares error of the whole identification results minimized, but also gives rise to minimal expectation of the identification (absolute or relative) error in an arbitrary damage parameter. A noteworthy thing is that such an optimal weight matrix (20) is attained when only measurement errors are considered, however, model errors may also exist which may change the optimal choice of the weight matrix; this is to be specifically investigated in what follows. 4.2. Additional considerations on model errors In addition to the damage parameters a that are to be identified, the model parameters that govern the two models (1) and (2) are also taken in account herein; they are designated as b ¼ ½b1 ; b2 ; . . . ; bm0 T . Throughout the paper, the model parameters b refer to the possible system parameters excluding the damage parameters. For structures, the model parameters b may include the mass density in the mass matrix M, the two Rayleigh damping coefficients a1 ; a2 in defining the Rayleigh damping matrix C, stiffness of the support in the stiffness matrix K and the amplitude of an impact load in the load vector f , etc. Generally, the already known model parameters ^ b are deemed to be accurate or at least to provide reasonably good approximations of the exact values bex so that they are directly applicable for damage identification with the algorithm in

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

9

Table 1; otherwise, they would be included into the damage parameters to be identified. For further analysis, the relationship b pertains to between the unknown exact model parameters bex and the known values ^

^b ¼ bex þ b

ð35Þ

where b represent the small perturbation errors subject to the Gaussian distribution with zero means and covariance E½b Tb  ¼ Q b . The model errors b are assumed to be independent from the measurement errors R . Following the sensitivity-based damage identification in Section 3, the derived data under the known (or estimated) model parameters ^ ¼ Rðaex ; bex Þ þ R with the measurement errors ^ bÞ, while for measured data, there holds R b, is extensively expressed as Rða; ^

R defined in Section 4.1; in this way, the goal function for damage identification becomes ^  Rða; ^bÞk2 gðaÞ ¼ kR W

ð36Þ

for which the minima or identified results are again denoted by a . To further estimate the covariance of the identification errors a  aex , similar derivations to those in Eqs. (17)–(19) are invoked by resorting to the linear Taylor expansions for a in the vicinity of aex ,

gðaÞ ¼ kRðaex ; bex Þ þ R  Rða; ^bÞk2W ¼ kRðaex ; bex Þ  Rðaex ; ^bÞ þ R þ Rðaex ; ^bÞ  Rða; ^bÞk2W  krb Rðaex ; ^bÞðbex  ^bÞ þ R þ Sðaex  aÞk2W ¼ kR  Sb b  S  ða  aex Þk2W where S ¼ ra Rðaex ; ^ bÞ; Sb ¼ rb Rðaex ; ^ bÞ :¼ ½

@Rðaex ;^ bÞ @b1

;

@Rðaex ;^ bÞ @b2

;...;

@bm0

. Minimization of (37) yields,

a  a  LW ðR  Sb b Þ 

ð37Þ

@Rðaex ;^ bÞ

ð38Þ

ex

where LW is of the same form with that in (18). Note that R  Sb b are also pertaining to Gaussian distribution with zero means and the following covariance,

E½ðR  Sb b ÞðR  Sb b ÞT  ¼ Q R þ Sb Q b STb :

ð39Þ

Analogous to the proof of the optimal weight matrix in Section 4.1, the optimal weight matrix herein under both measurement and model errors is easily found to be 1

Wopt ¼ cðQ R þ Sb Q b STb Þ ;

c > 0:

ð40Þ

Remark 3. Unlike the optimal weight matrix (20) where only measurement errors are concerned, the optimal weight matrix (40) under both measurement and model errors is no longer independent from the model problem; rather, due to the term Sb , the optimal weight matrix (40) should be obtained with additional sensitivity analysis to the model parameters. ^ because the required aex is unknown and could not However, difficulty would arise for sensitivity analysis of Sb ¼ rb Rðaex ; bÞ be known exactly in practice. To circumvent this difficulty, there are two ways. One easy but rough way is to approximately use the undamaged parameters að1Þ as aex in sensitivity analysis; this is often effective when damage is small. In another ~  , that are obtained through damage identification (36) with the weight matrix way, the approximate damage parameters a

W ¼ cQ 1 and are often good approximations of aex , are adopted for computation of the sensitivity matrix Sb and the R optimal weight matrix (40) thereof. It is also noteworthy that in case when the response data is insensitive to the model

parameters, i.e., all singular values of Sb are small enough, one can approximately set Wopt  cQ 1 R in practice and for this reason, Wquasi ¼ cQ 1 R is termed the quasi-optimal weight matrix when model errors exist. bÞ requires additional sensitivity analysis of the model Eqs. (1) and (2) to the model Remark 4. Computation of Sb ¼ rb Rðaex ; ^ parameters b. To illustrate, three typical kinds of model parameters are taken into account for example,
Support stiffness (see Fig. 1(a)). Fixed boundary condition is often encountered in analysis of beam structures, however, this kind of boundary condition is rather idealized and is practically treated as the elastic boundary condition with large support stiffness including translational stiffness kt ; kv and rotational stiffness kh . For simplicity, three normalized model ^t ð1 þ et Þ; kv ¼ k ^v ð1 þ ev Þ; kh ¼ k ^h ð1 þ eh Þ parameters et ; ev ; eh regarding the support stiffness are introduced such that kt ¼ k ^t ; k ^v ; k ^h represent the known values of support stiffness with possible errors and et ; ev ; eh indicate possible relative where k errors. For the normalized model parameters, the known values are simply found to be ^ et ¼ 0; ^ ev ¼ 0; ^eh ¼ 0. As with the sensitivity analysis, only sensitivity to et is presented and sensitivity to other two parameters ev ; eh is analyzed likewise. To proceed further, assume that the model parameter et corresponds to the kthe DOF (or the horizontal displacement of ^t and zeros elsewhere. Then, to acquire the support) and then, there is @K ¼ Kt , a matrix with K t ðk; kÞ ¼ k @et

@R @et

€ @d

@K @di ¼ ½@e ; ; i 2 S u ; @etj ; j 2 S u€ , the sensitivity analysis for the dynamic Eqs. (1) and the modal Eqs. (2) is invoked in the t @et

same way with those in Eqs. (7) and (8) except that Kk now is replaced by Kt .

10

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 1. Schematic illustrations for support stiffness and impact load.


Damping parameters. In the dynamic Eqs. (1), the Rayleigh damping assumption with two coefficients a1 ; a2 is used. Usually, a1 ; a2 can be determined using two different eigenfrequencies and corresponding modal damping factors through modal tests; they may admit small perturbations. For sensitivity analysis is

@K @ak

@R @ak

€ @d

@K @di ¼ ½@a ; @a ; i 2 S u ; @a j ; j 2 S u€ ; k ¼ 1; 2, there k

k

k

¼ 0 since the damping coefficients are not involved in the eigenvalue Eq. (2) of an undamped system and

@di @ak

€ @d

; @a j

k

@u are obtained through the similar sensitivity analysis to that in (7) except that the main variable herein becomes @a k and the right-hand end should be replaced by zk ¼ Mu_ for k ¼ 1 and zk ¼ Ku_ for k ¼ 2. That is to say, the following € @u @ u_ @u @u @ u_ ðtÞ þ C @a ðtÞ þ K @a ðtÞ ¼ zk ðtÞ; k ¼ 1; 2 with @a ð0Þ ¼ @a ð0Þ ¼ 0. dynamic equation is invoked as M @a k

k

k

k

k


Load parameters (refer to Fig. 1(b)). The impact load is taken into account for example because it can invoke the highfrequency information of the structure [12,15] and therefore, is beneficial for damage identification. Assume that the peak ^ þ rÞ where A ^ is the known peak value and in value A concerning the total energy of the impact load is of the form A ¼ Að1 this way, the load vector shall be parameterized as f ¼ ^ f ð1 þ rÞ with the model parameter r verifying ^r ¼ 0. Usually, for impact tests, the initial displacements and velocities are zero, i.e., u0 ¼ 0; u_ 0 ¼ 0. Then, sensitivity analysis is given as fol¼ 0; for the dynamic Eqs. (1), sensitivity to r yields lows: for the modal Eqs. (2), there is obviously @K @r

€ @u @ u_ @u @f ^ ðtÞ þ C ðtÞ þ K ðtÞ ¼ ¼ f; @r @r @r @r @u @ u_ ð0Þ ¼ 0; ð0Þ ¼ 0 @r @r

M

meaning that

0

@u : ^ ¼u @r

t > 0;

after comparison between (1) and (41). To conclude, the sensitivity

1

0 @R : B d ^ C ¼@ i ; i 2 S u A @r ^€ dj ; j 2 S u€

ð41Þ

@R @r

is simply obtained as

ð42Þ

without additional computations.

Remark 5. With the optimal weight matrix (40), the enhanced response sensitivity approach in Table 1 can directly result in the identified damage parameters a , nevertheless, the model parameters can also be updated or corrected. The procedure is given as follows,

^  Rða ; bÞk2 1 þ kb  ^bk2 1 b ¼ arg minkR Q Q b

R

b

ð43Þ

for which the optimal weight matrix as the reciprocal of covariance of the measured data (refer to the derivation in Section 4.1) is applied. Such kind of update procedure can be also found in [2,17]. As a result of (43), the corrected model parameters are

: 1 1 T 1 ^  ^ b ¼^b þ ½STb Q 1 R Sb þ Q b  Sb Q R ðR  Rða ; bÞÞ where Sb has already been obtained in Remark 3.

ð44Þ

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

11

5. Numerical tests Two numerical examples concerning a simple plane frame and a rectangular Mindlin-Reissner plate are studied to verify the proposed optimal weight matrix under measurement and model errors and to see the specific effect on the choice of the weight matrix for damage identification. The algorithm in Table 1 is adopted for damage identification and the algorithmic parameters RelErr; Nmax; qcr ; c; Ntr are taken as the default values in the following brackets. For model parameters, the known parameters ^ b are given and the exact values bex are obtained as

^ bex k ¼ bk ð1 þ eb Rand Þ;

k ¼ 1; 2; . . . ; m0

ð45Þ

where Rand is a random number pertaining to standard normal distribution and eb represents the model error level for which eb ¼ 0 means no model error. For later simulation, ^ b would keep unchanged for an example while from (45), bex may change. As with measurement noises, there are

^ 2k ¼ ðxk Þ2 ð1 þ eK Rand Þ; k ¼ 1; 2; . . . ; q x ^ ¼ d þ eu  R d i i andu  stdðdi Þ; ^ € € € Þ; dj ¼ dj þ ea  Randa  stdðd j

i 2 Su

ð46Þ

j 2 S u€

€ are simulated data taken from Rðaex ; bex Þ; eK ; eu ; ea are levels of measurement errors, stdðÞ denotes where ðxk Þ2 ; di ; d j the standard deviation of a vector and Rand ; Randu ; Randa are random number and vectors pertaining to standard normal distribution. As with the constant c in the weight matrices (20) and (40), it is set to c ¼ 1 for later damage identification. id In the examples, the optimal weight matrix Wopt ¼ cQ 1 R and the identity weight matrix W ¼ I will be considered when there are only measurement errors. While when additional model errors exist, the optimal weight matrix 1

id Wopt ¼ cðQ R þ Sb Q b STb Þ , the quasi-optimal weight matrix Wquasi ¼ cQ 1 R and the identity weight matrix W ¼ I are taken into account for comparison. Furthermore, to measure how much a weight matrix A is similar/approximate to another weight matrix B, a similarity indicator cosðA; BÞ is defined as follows,

trðAT BÞ cosðA; BÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ½0; 1: trðAT AÞtrðBT BÞ

ð47Þ

Usually, greater cosðA; BÞ means that A is more approximate to B. Moreover, cosðA; BÞ ! 1 indicates that A and B are almost the same, while cosðA; BÞ ! 0 implies that A and B are completely different. 5.1. A plane frame A single-span plane frame structure [3,10] as shown in Fig. 2 is studied herein to see the effect of the weight matrix in damage identification with measurement and/or model errors. Geometric parameters are: the height and width of the frame H ¼ 1:2 m and L ¼ 0:6 m, and the rectangular cross-section of dimensions b ¼ 0:01 m and h ¼ 0:02 m with the inertial moment and area being I ¼ bh =12 ¼ 6:67 109 m4 and A ¼ bh ¼ 2 104 m2. The material properties of the intact frame 3

are: the mass density q ¼ 2:7 103 kg/m3 and the elastic modulus E ¼ 6:9 1010 N/m2. As for the finite element modeling, the frame is uniformly divided into 11 Euler-Bernoulli beam elements with 12 nodes each with three DOFs where x-, ytranslational and rotational displacements at the kth node are denoted by uk ; v k ; hk . Specific expressions of the elementary

mass and stiffness matrices MðkÞ ; KðkÞ for the kth element can be found in [10,25]. Practically, the fixed nodes 1 and 12 of ^v ¼ 1:5 1010 kN/m and k ^h ¼ 1:5 109 ^t ¼ k the frame are modeled with the same large support stiffness (see Fig. 2(c)): k kNm/rad. Eventually, the global stiffness and mass matrices are ðkÞ K ¼ A11 k¼1 K ;

Kð1Þ

ðkÞ M ¼ A11 k¼1 M ; ^t ; k ^v ; k ^h ; 0; 0; 0Þ; ¼ Kð1Þ þ diagð½k

^ ;k ^ ;k ^ Kð11Þ ¼ Kð11Þ þ diagð½0; 0; 0; k t v h Þ where A is the assembly operator and diagðÞ denotes diagonalization of a vector. To get the Rayleigh damping matrix, two different natural frequencies xi – xj in a reasonable range as well as corresponding modal damping ratios fi ; fj are often tested and measured, and then, the two damping coefficients are calculated [26] as



a1 a2



¼

2xi xj x2j  x2i



xj xi 1=xj 1=xi



fi fj



:

ð48Þ

12

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 2. A plane frame model: mesh, load and boundary conditions.

In this example, the first two frequencies of the intact frame are found to be 13.095 and 57.308 Hz and for simplicity, the first ^1 ¼ 1:3395s1 and two modal damping ratios are both assumed to be 0.01. Then, by (48), the Rayleigh coefficients are a ^1 M þ a ^2 K is obtained. As for the external load, the horizontal impact load ^2 ¼ 4:52 105 s so that the damping matrix C ¼ a a ^ FðtÞ lasting for 0.1s duration as shown in Fig. 2(b) is enforced on node 2. ^v ; k ^h ; a ^t ; k ^ are directly used to get the simulated ^1 ; a ^2 ; FðtÞ In the absence of model errors, the measured model parameters k data; while if model errors are considered, the exact (or simulated) model parameters are obtained as in Eq. (45) and then, numerical simulation is performed under the exact model parameters. For this example, the simulated data—displacements or accelerations is obtained from numerical simulation during the time period [0, 4s] with a sampling rate of 500 Hz and the measured data which is directly used for damage parameter identification is obtained by addition of measurement errors at certain levels (e.g.,=5%) as in (46). As for the choice of the weight matrix, if no model error exists, the optimal weight matrix id Wopt ¼ cQ 1 R in (20) and the conventional identity weight matrix (W ¼ I) [13] are both used for comparison; while if model

errors are considered, three weight matrices—the identity weight matrix, the quasi-optimal weight matrix Wquasi ¼ cQ 1 R 1

only accounting for measurement errors and the optimal weight matrix Wopt ¼ cðQ R þ Sb Q b STb Þ in (40) shall be adopted for comparative purpose. Consequently, with the measured data obtained and the weight matrix defined, the enhanced response sensitivity approach in Table 1 can give rise to the identified damage parameters. It is noteworthy that due to the randomness of the (measurement or model) errors and the requirement to compute the expectation of squares of the identification errors E½ða  aex Þða  aex ÞT , 40 sets of the measurement data (and exact model parameters data when considering model errors) will be generated by Monte Carlo sampling. For damage identification, the damage is assumed to be certain reduction of elastic modulus and therefore, the elastic modulus in the 11 elements are treated as the damage parameters. In what follows, five damage scenarios with measurement and/or model errors and using different measured quantities that are listed in Table 2 with ‘—’ meaning no model error are studied in detail along the above lines. Scenarios 1–3 use accelerations and/or displacements as the measurement data and consider different measurement error levels for accelerations/displacements at different positions (see Table 2); the damages are moderate (around 20%) and even large (40%). The three scenarios are assumed to be free from model errors. Based on the Monte Carlo sampling and the damage identification algorithm in Table 1, damage parameters can be identified under the optimal weight matrix Wopt ¼ cQ 1 R and the identity weight matrix Wid ¼ I and detailed results on means and standard deviations of individual damage parameters are schematically exhibited in Figs. 3–5 where ‘scaled identified results’ refer to

ak aintact k

with aintact denoting damage k

parameters of the intact frame. The benefits brought by the optimal weight matrix are obvious since for all the three scenarios, more accurate averages and less standard deviations are obtained under the optimal weight matrix. Moreover, to quantify the accuracy of the identified results, the ‘relative error’ for the kth parameter is introduced as

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u "  2 # u ak  aex k t

100% ‘relative error’ ¼ E ex

ak

ð49Þ

13

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24 Table 2 Damage identification scenarios for plane frame structure. Damage scenario

Damage location (element No.)

Reduction in elastic modulus

Measurement (with noise)

Model error

1 2 3 4 5

3 4, 4, 3, 3,

40% 20%, 20%, 20%, 20%,

€ 2 ð5%Þ; u € 9 ð10%Þ u € 2 ð2%Þ; u € 4 ð10%Þ; u € 9 ð10%Þ u € 2 ð2%Þ; u4 ð10%Þ; u € 9 ð10%Þ u € 9 ð5%Þ € 2 ð2%Þ; u u € 2 ð2%Þ; u € 9 ð5%Þ u

– – – FðtÞ (4%) kt ; kv ; kh (10%)

5 5 6, 9 6, 9

25% 25% 15%, 20% 15%, 20%

Fig. 3. Bar graph with means and standard deviations for damage identification of Scenario 1.

Fig. 4. Bar graph with means and standard deviations for damage identification of Scenario 2.

where ak denotes identified parameter. Then, the ‘relative errors’ for all identified damage parameters of the three scenarios are also computed and presented in Table 3 with the largest ‘relative error’ underlined.

14

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 5. Bar graph with means and standard deviations for damage identification of Scenario 3.

Table 3 ‘Relative errors’ for identified elastic modulus in all elements (#k) of Scenarios 1–3. Scena.

Weight

#1

#2

#3

#4

#5

#6

#7

#8

#9

#10

#11

1

Wid

0.45%

0.96%

0.67%

1.86%

1.17%

1.08%

4.49%

2.94%

3.60%

2.06%

1.53%

Wopt

0.20%

0.57%

0.44%

1.20%

0.72%

0.80%

3.42%

1.61%

2.33%

0.98%

0.62%

Wid

0.64%

1.72%

2.30%

2.58%

4.70%

1.31%

3.42%

2.95%

2.98%

2.12%

1.49%

Wopt

0.19%

0.57%

0.75%

0.50%

1.67%

0.28%

0.93%

1.28%

0.53%

1.13%

0.30%

Wid

0.80%

2.32%

3.67%

2.43%

8.89%

1.52%

4.70%

3.69%

3.53%

2.91%

1.79%

Wopt

0.14%

0.47%

0.67%

0.38%

1.54%

0.24%

0.86%

1.01%

0.42%

0.93%

0.26%

2

3

As is seen, the maximal ‘relative errors’ reach respective 4.49%, 4.70% and 8.89% for the three scenarios 1–3 under the identity weight matrix and decrease to 3.42%, 1.67% and 1.54% when the optimal weight matrix is used; both weight matrices shall lead to reasonably good identification of damage parameters in the absence of model errors. Nevertheless, careful comparisons show that the ‘relative errors’ by the optimal weight matrix of all damage parameters are less than those corresponding to the identity weight matrix and this is reasonable by referring to the analysis in Remark 2; the optimal weight matrix has led to evidently superior damage identification results to the identity matrix in all three scenarios. The advantage of the optimal weight matrix over the identity matrix is particularly substantial for scenario 3 where both accelerations and displacements are measured since usage of the optimal weight matrix renders the ‘relative errors’ almost all reduced to 1/4 of those by the identity weight matrix. Indeed, the optimal choice of the weight matrix can work as expected to improve the accuracy of the damage identification results when different levels of measurement errors or different types of response data are invoked. Scenarios 4–5 adopt accelerations as measurements and additionally take the model errors into account; in Scenario 4, the noise of level 4% is added to the impact load while the support stiffness at node 1 is assumed to admit the model error of level 10% in Scenario 5. Practically, how to analyze such model parameters and incorporate model error covariance into the optimal weight matrix (40) has already been illustrated in Remarks 3 and 4. Again, along with the Monte Carlo sampling and the enhanced response sensitivity approach, identification results for 40 sampling data sets are obtained for both scenarios with the three specific weight matrices—the identity matrix Wid ¼ I, the quasi-optimal weight matrix Wquasi ¼ cQ 1 R 1

and the optimal weight matrix Wopt ¼ cðQ R þ Sb Q b STb Þ . Based on the 40 sets of identification results, statistical bar graphs for the two scenarios are depicted in Figs. 6 and 7 and the specific ‘relative errors’ of the identified parameters are also calculated and presented in Table 4. For Scenario 4, the maximal ‘relative error’ for the identity weight matrix arrives at 22.20%, approaching the actual damage levels and therefore, implying that damage may become unidentifiable even though the means in Fig. 6 are satisfactory, while the maximal ‘relative error’ is reduced to 6.48% for the quasi-optimal weight matrix and 1.59% for the optimal weight

15

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 6. Bar graph with means and standard deviations for damage identification of Scenario 4.

Fig. 7. Bar graph with means and standard deviations for damage identification of Scenario 5.

Table 4 ‘Relative errors’ for identified elastic modulus in all elements (#k) of Scenarios 4–5. Scena. 4

5

#1

#2

#3

#4

#5

#6

id

W

9.36%

6.98%

8.31%

4.82%

5.98%

3.80%

Wquasi

6.45%

1.19%

1.24%

2.78%

Wopt

1.21%

0.41%

6.48% 1.34%

0.52%

Wid

0.22%

0.52%

0.61%

Wquasi

0.12%

0.31%

0.41%

Wopt

0.12%

0.31%

0.41%

Weight

#7

#8

#9

#10

#11

18.69%

1.75%

11.96%

0.92%

22.20% 4.07%

8.41% 2.04%

5.85%

3.13%

6.03%

0.98%

0.38%

1.59%

1.01%

1.57%

0.89%

1.18%

0.63%

1.43%

0.46%

2.12%

1.40%

1.29%

0.92%

0.39%

0.43%

0.92%

0.41%

1.49%

1.17%

1.17%

0.84%

0.38%

0.43%

0.92%

0.41%

1.48%

1.17%

1.17%

0.84%

0.38%

16

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

matrix; these two levels of identification errors are well acceptable for damage identification. This indicates that with the load error, the identity weight matrix may lead to unacceptable identification errors, while the quasi-optimal weight concerning merely the measurement errors can give rise to reasonable good identification results and nearly perfect identification results are gained using the optimal weight matrix where both the measurement and model errors are considered. More specifically, the ‘relative error’ of an arbitrary damage parameter under the optimal weight matrix is substantially less than that under the other two weight matrices; this to some extent verifies the optimality of the choice (40). For Scenario 5, it is concluded from Fig. 7 and Table 4 that the ‘relative errors’ under the quasi-optimal weight matrix and the optimal weight matrix are almost the same, but are absolutely less than those under the identity weight matrix; nevertheless, damages are all well identified using the three weight matrices. The model errors of the large support stiffness are found to have little influence on the identification accuracy. To get more insights into the difference between the identification results for Sce^ and the stiffness €2 ; u € 9 to the load parameter r with FðtÞ ¼ ð1 þ rÞ  FðtÞ narios 4 and 5, sensitivity of the measured responses u ^t ð1 þ et Þ; kv ¼ k ^v ð1 þ ev Þ; kh ¼ k ^h ð1 þ eh Þ is additionally analyzed and shown in Fig. 8. As is parameters et ; ev ; eh with kt ¼ k seen, the measured responses are insensitive to the large stiffness parameters while are quite sensitive to the load parameter; actually, similar results have already been reported in [10]. This to some extent explains why the identification results are sensitive to the load error but are insensitive to errors in large stiffness and why the load error should be much more counted than the large stiffness errors in the choice of optimal weight matrix. Identification results for Scenario 4 show that the load error and the corresponding optimal weight matrix thereof play an important role in damage identification. In fact, the 40 sets of sampling load parameters involved above can be corrected along the lines of (44) in Remark 5 and on the basis of the identified results under the optimal weight matrix. Detailed results on the peak values of the corrected and the exact (refer to (45)) impact loads are displayed in Fig. 9. It is seen that the relative jbk bex j k — jbex j k

errors between the corrected parameters and the exact ones—

are all less than 1%, while may exceed 5% between the

measured and exact load parameters; the correction in (44) evidently results in better approximations of the actual load/model parameters. The above results have shown that the choice of the weight matrix has played a significant role in the accuracy and robustness of the damage identification results. The optimal weight matrix derived in Section 4 has been shown to work as well as it was expected. To further see the similarity of the optimal weight matrix to the identity weight matrix or the quasi-optimal weight matrix, the similarity indicator (47) is calculated for all the five scenarios and the detailed results are listed in Table 5. For Scenarios 1–3 where only measurement errors are enforced, the largest cosðWopt ; Wid Þ occurring for Scenario 1 means that the identity weight matrix Wid is most approximate to the optimal weight matrix Wopt ; this to some extent explains in Table 3 why the improvement of the ‘relative error’ by the optimal weight matrix for Scenario 1 is less obvious than those for the other two scenarios. For Scenarios 4–5, the indicator cosðWopt ; Wquasi Þ tending to 1 is obviously greater than cosðWopt ; Wid Þ, indicating that the quasi-optimal weight matrix Wquasi is indeed a good approximation of the optimal weight matrix Wopt when model errors exist.

€ 2 and (b) u € 9 to stiffness and load parameters in the frame. Fig. 8. Sensitivity of measurements—(a) u

17

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 9. Corrected peak loads versus Monte Carlo sampled exact peak loads of the frame.

Table 5 Values of similarity indicators for the weight matrices used for Scenarios 1–5. Damage scenario

1

2

3

4

cosðWopt ; Wid Þ

0.7300

0.5822

0.5774

0.7222

5 0.7217

cosðWopt ; Wquasi Þ

–

–

–

0.9998

0.9996

Fig. 10. A simply-supported rectangular plate.

5.2. A rectangular plate The second example concerns damage identification in a simply supported rectangular plate for which the geometry, finite element mesh and external impact loads are shown in Fig. 10. Material properties of the intact plate are: Young’s modulus E ¼ 25 GPa, mass density q ¼ 2:8 103 kg/m3 and Poisson ratio m ¼ 0:2. Herein, the Young’s modulus is assumed to be

18

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

constant in each element and is again treated as damage/stiffness parameters. For practical analysis, the plate pertains to the Mindlin-Reissner plate theory for which the mass and stiffness matrices arising from the finite element analysis [3,25] are obtained as 3

K¼

h 12 Z

M¼

X

Z X

eb ð/ÞT Cb eb ð/ÞdX þ jh

Z X

es ð/ÞT Cs es ð/ÞdX;

/T Cm /dX;

0 1 m E B Cb ¼ @m 1 1  m2 0 0

0

1

0 C A; 1m 2

E Cs ¼ 2ð1 þ mÞ



1 0 0

1

 ;

0 B Cm ¼ q@

3

h =12

0

0

0

h

1

C 3 h =12 0 A

0 0

where h ¼ 0:06 m is thickness of the plate, j ¼ 5=6 is the shear correction factor, / collects the bilinear finite element shape P functions /i such that ½hx ; hy ; w ¼ /u :¼ ni¼1 /i ½ðhx Þi ; ðhy Þi ; wi  with ðhx Þi ; ðhy Þi ; wi being rotations and deflection DOFs at node i and

eb ð½hx ; hy ; wÞ ¼ ½@h@xx ; @h@yy ; @h@yx þ @h@xy ; es ð½hx ; hy ; wÞ ¼ ½@w  hx ; @w  hy  are the respective in-plane strain and out-plane shear @x @y

^1 ¼ 1:6118 s1 and strain. As for the Rayleigh damping matrix, the two Rayleigh coefficients are taken as a 5 ^2 ¼ 6:08 10 s. a Four damage scenarios using acceleration and/or eigenfrequency data with measurement and/or model errors are considered in this example; the details are given in Table 6. The measured data is obtained with the sampling rate of 1000 Hz as demonstrated at the beginning of Section 5 and 20 sets of sampling data are generated by Monte Carlo simulation for each case. Analogous to the first frame example, the definitions of Wid ; Wopt for parameter identification with only measurement errors and Wid ; Wquasi ; Wopt under model errors are also adopted for convenience purpose. Then, the algorithm in Table 1 is adopted for damage identification. Scenarios 6–8 try to study the effect of the weight matrix in hybrid sensitivity-based damage identification [13] where the response data and/or the eigenfrequency data are used. The three scenarios are free from model errors. In Scenarios 6–7, the optimal weight matrix and possible identity weight matrix are used while for Scenario 8 with hybrid measured data—the response and eigenfrequency data, both the identity weight matrix and the optimal weight matrix are considered; then, detailed identification results are displayed in Figs. 11–13. For Scenario 6 with the response data, though average identification results are almost of the same good quality for both the identity weight matrix and the optimal weight matrix, much smaller individual standard deviations of identified damage parameters are obtained under the optimal weight matrix, indicating that superiority of the optimal weight matrix is well exhibited. The eigenfrequency data alone in Scenario 7 even with the optimal weight matrix is unable to locate and quantify the damage (see Fig. 12), which is reasonable since the amount of the measured eigenfrequency data is too few with respect to the complexity of the damage parameters to identify the damage. A remarkable thing is that the hybrid data with the identity weight matrix in Scenario 8 gives unsatisfactory results, which are much worse than the lonely usage of the response data in Scenario 6; only when the optimal weight matrix is incorporated with the hybrid data in Scenario 8, the accuracy would become slightly better than that of Scenario 6. In experience, more measured data would lead to better identification results, however, it is concluded herein that this is achieved only when the optimal weight matrix is used. In this way, the importance of the choice of the weight matrix is evidently seen. To see the quantitative benefits brought by the optimal weight matrix, the ‘relative errors’, defined as in (49), are further computed for Scenarios 6–8; the results are graphically shown in Figs. 14–16. Obviously, the maximal ‘relative errors’ for Scenario 7 under the optimal weight matrix and Scenario 8 under the identity weight matrix reach 19.33% and 18.91% which are almost at the same level of damage intensity; this quantitatively means failure of the damage identification procedure. For Scenario 6, the maximal ‘relative error’ approaches 7.02% under the identity weight matrix and is reduced to 3.90% with the optimal weight matrix; the superior performance of the optimal weight matrix as illustrated in Remark 2 is assured. Moreover, the hybrid data with the optimal weight matrix leads to the best identification results among Scenarios

Table 6 Damage identification scenarios for rectangular plate. Damage scenario

Damage location (element No.)

Reduction in elastic modulus

Measurement (with noise)

Model error

6 7

1, 16, 18, 23

10%, 15%, 20%, 10%

€ 20 ð1%Þ; w € 37 ð5%Þ w

– –

8 9

5

fx2i gi¼1 ð0:15%Þ € 37 ð5%Þ; fx2i g5i¼1 ð0:15%Þ € 20 ð1%Þ; w w € 37 ð5%Þ € 20 ð1%Þ; w w

– a1 ; a2 ð5%Þ

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

19

Fig. 11. Bar graph with means and standard deviations for damage identification of Scenario 6.

Fig. 12. Bar graph with means and standard deviations for damage identification of Scenario 7.

6–8. Indeed, optimal identification accuracy is achieved by the optimal weight matrix, and what’s more, the optimal weight shall be particularly important in hybrid sensitivity-based damage identification. Scenario 9 attempts to pose the effects of the damping coefficient model errors and the corresponding optimal weight ^1 ; a ^2 are fixed and then, actual values of the 20 Monte matrix. To get the simulated data, the measured damping parameters a

20

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 13. Bar graph with means and standard deviations for damage identification of Scenario 8.

Fig. 14. ‘Relative errors’ for identified elastic modulus in all elements of Scenario 6 under (a) identity weight and (b) optimal weight.

Fig. 15. ‘Relative errors’ for identified elastic modulus in all elements of Scenario 7 under optimal weight.

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

21

Fig. 16. ‘Relative errors’ for identified elastic modulus in all elements of Scenario 8 under (a) identity weight and (b) optimal weight.

Fig. 17. Bar graph with means and standard deviations for damage identification of Scenario 9.

Carlo sets are obtained as in (45). Measurement noises are enforced as in (46). Then, damages are identifiable under the three weight matrices—Wid ; Wquasi ; Wopt . Detailed identification results are shown in Fig. 17 and corresponding ‘relative errors’ are displayed in Fig. 18. As is seen, the average identification results under the three weight matrices can well reflect the damage locations and intensities, however, the standard deviations are quite different. A first sight finds that the quasi-optimal weight matrix considering only the measurement error can give rise to better results than the identity weight matrix, but besides the two weight matrices, the optimal weight matrix has led to much more accurate averages and much less standard deviations. This sight is again perfectly testified by the results of ‘relative errors’ in Fig. 18 where the maximal ‘relative errors’ are respective 35.01%, 16.70% and 5.04% for the three weight matrices Wid ; Wquasi ; Wopt . All results not only show the optimal accuracy brought by the optimal weight matrix, but also pose the necessity of incorporating the damping coefficient model errors into the optimal weight matrix. To see why the damping coefficient errors have so huge influence on optimal weight selection and damage identification, ^k ð1 þ ek Þ; k ¼ 1; 2 is graph€ 20 ; w € 37 to the two model parameters e1 ; e2 with ak ¼ a the sensitivity of the measured quantities w ically plotted in Fig. 19. Obviously, the response data is sensitive to the errors in the Rayleigh damping coefficients and this just explains why damage identification results are sensitive to the Rayleigh damping coefficient errors in the example. In ^1 ; a ^2 addition, with satisfactory identified damage parameters under the optimal weight matrix, the damping coefficients a

22

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 18. ‘Relative errors’ for identified elastic modulus in all elements of Scenario 9 under (a) identity weight, (b) quasi-optimal weight and (c) optimal weight.

€ 20 and (b) w € 37 to damping parameters in the plate. Fig. 19. Sensitivity of measurements—(a) w

can be further corrected by resorting to Remark 5. Details on corrected results corresponding to the 20 sets of Monte Carlo sampled damping coefficients are depicted in Fig. 20. Again, the correction in (44) gives nearly perfect approximations to the exact/simulated damping coefficients.

23

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24

Fig. 20. Corrected damping coefficients—(a) a1 and (b) a2 versus Monte Carlo sampled exact values of the plate.

Table 7 Values of similarity indicators for the weight matrices used for Scenarios 6–9. Damage scenario opt

cosðW

id

;W Þ

cosðWopt ; Wquasi Þ

6

7

8

0.7193

0.4711

0.7189

9 0.7190

–

–

–

0.9995

Finally, the similarity indicator (47) is again calculated to see the similarity of the optimal weight matrix to the identity weight matrix and the quasi-optimal weight matrix for the four scenarios; the results are given in Table 7. For Scenarios 6–8 where only measurement errors are considered, the similarity indicators cosðWopt ; Wid Þ are all less than 0.72, meaning that the optimal weight matrix and the identity weight matrix are quite different. For Scenario 9 with both measurement and model errors, the similarity indicator cosðWopt ; Wquasi Þ ¼ 0:995 is quite near to 1, implying that the quasi-optimal weight matrix is a reasonably good approximation of the optimal weight matrix. 6. Conclusions In this paper, how to choose the weight matrix for response sensitivity-based damage identification has been throughout discussed. By resorting to the criterion that the optimal weight matrix should minimize the expectation of squares error in identification results (see Eq. (16)), the optimal weight matrix is found to be inverse proportional to the covariance of measurement errors when only measurement errors are considered. In case of model errors, the same criterion shows that the optimal weight matrix depends not only on measurement errors, but also on model errors along with sensitivity to model parameters. Further analysis finds that the optimal weight matrix can also render the ‘relative error’ (see Eq. (49)) of every individual damage parameter minimal. Then, numerical tests on a plane frame and a simply-supported plate, using response and/or eigenfrequency data with measurement and/or model errors, are conducted and results show that.
The optimal performance that the ‘relative error’ of every individual damage parameters is minimized is indeed reached by the proposed optimal weight matrices under measurement and/or model errors;
The benefits brought by the optimal weight matrix are great especially when hybrid types of data—acceleration response, displacement response and/or eigenfrequencies are used;
The response data is sensitive to load and damping parameters, but is insensitive to large support stiffness;

24

Z.-R. Lu et al. / Mechanical Systems and Signal Processing 114 (2019) 1–24


For model errors that the response data is sensitive to, the optimal weight matrix involving measurement error covariance, model error covariance and response sensitivity to model parameters is preferred in order to get satisfactory identification results;
For model errors that the response data is insensitive to, the quasi-optimal weight matrix is suggested for its convenient accessibility and nearly optimal performance;
For model parameters with errors, the correction procedure in Remark 5 can substantially correct the measured model parameters. Acknowledgement The present investigation was performed under the support of National Natural Science Foundation of China (No. 11572356 and No. 11702336), Guangdong Province Natural Science Foundation (No. 2017A030313007) and the Fundamental Research Funds of the Central Universities (No. 17lgjc42 and No. 17lgpy54). References [1] W. Fan, P. Qiao, Vibration-based damage identification methods: a review and comparative study, Struct. Health Monitor. 10 (2011) 83–111. [2] M.I. Friswell, J.E. Mottershead, Finite Element Model Updating in Structural Dynamics, Kluwer Academic Publishers, Dordrecht, 1995. [3] Z.R. Lu, L. Wang, An enhanced response sensitivity approach for structural damage identification: convergence and performance, Int. J. Numer. Methods Eng. (2017), doi:https://doi.org/10.1002/nme.5502. [4] J.E. Mottershead, M. Link, M.I. Friswell, The sensitivity method in finite element model updating: a tutorial, Mech. Syst. Signal Process. 25 (2011) 2275– 2296. [5] G.R. Liu, S.C. Chen, Flaw detection in sandwich plates based on time-harmonic response using genetic algorithm, Comput. Methods Appl. Mech. Eng. 190 (2001) 5505–5514. [6] Z.H. Ding, M. Huang, Z.R. Lu, Structural damage detection using artificial bee colony algorithm with hybrid search strategy, Swarm Evol. Comput. 28 (2016) 1–13. [7] C. Farhat, F.M. Hemez, Updating finite element dynamic models using an element-by-element sensitivity methodology, AIAA J. 31 (1993) 1702–1711. [8] H.P. Chen, T.S. Maung, Regularised finite element model updating using measured incomplete modal data, J. Sound Vib. 333 (2014) 5566–5582. [9] P.G. Bakir, E. Reynders, G. De Roeck, Sensitivity-based finite element model updating using constrained optimization with a trust region algorithm, J. Sound Vib. 305 (2007) 211–225. [10] Z.R. Lu, S.S. Law, Features of dynamic response sensitivity and its application in damage detection, J. Sound Vib. ;303 (2007) 305–329. [11] Y.Z. Fu, Z.R. Lu, J.K. Liu, Damage identification in plates using finite element model updating in time domain, J. Sound Vib. ;332 (2013) 7018–7032. [12] M. Link, M. Weiland, Damage identification by multi-model updating in the modal and in the time domain, Mech. Syst. Signal Process. 23 (2009) 1734– 1746. [13] Z.R. Lu, X.X. Lin, Y.M. Chen, M. Huang, Hybrid sensitivity matrix for damage identification in axially functionally graded beams, Appl. Math. Model. 41 (2017) 604–617. [14] B. Titurus, M.I. Friswell, Regularization in model updating, Int. J. Numer. Methods Eng. 75 (2008) 440–478. [15] J. Cattarius, D.J. Inman, Time domain analysis for damage detection in smart structures, Mech. Syst. Signal Process. 1 (3) (1997) 409–423. [16] T.A.N. Silva, N.M.M. Maia, M. Link, J. Mottershead, Parameter selection and covariance updating, Mech. Syst. Signal Process. 7–71 (2016) 269–283. [17] Y. Govers, M. Link, Stochastic model updating–covariance matrix adjustment from uncertain experiment modal data, Mech. Syst. Signal Process. 24 (2010) 696–706. [18] K.A. Brownlee, Statistical Theory and Methodology in Science and Engineering, John Wiley & Sons, 1965. [19] J.J. More, The Levenberg-Marquardt algorithm: implementation and theory, Chapter Numerical analysis, volume 630 of the series Lecture Notes in Mathematics, 1978, pp: 105–116. [20] A.M. Tikhonov, On the solution of ill-posed problems and the method of regularization, Soviet Math. 4 (1963) 1035–1038. [21] P.C. Hansen, Analysis of discrete ill-posed problems by means of the L-curves, SIAM Rev. 34 (4) (1992) 561–580. [22] P.C. Hansen, Regularization tools—a matlab package for analysis and solution of discrete ill-posed problem, Numer. Algorithms 6 (1) (1994) 1–35. [23] S.P. Won, F. Golnaraghi, A triaxial accelerometer calibration method using a mathematical model, IEEE Trans. Instrum. Meas. 59 (8) (2010) 2144–2153. [24] M. Bornert, F. Bremand, P. Doumalin, et al, Assessment of digital image correlation measurement errors: methodology and results, Exp. Mech. 49 (2009) 353–370. [25] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method for Solid and Structural Mechanics, sixth ed., Elsevier, 2005. [26] A.K. Chopra, Dynamics of Structures, Prentice Hall, 2000.