Journal of Natural Disaster Science, Volume 24, Number 2, 2002, pp101-107

Health Monitoring Algorithm by the Monte Carlo Filter Basedon Non-Gaussian Noise
Ikumasa YOSHIDATokyo Electric Power Services Co.,Ldt., Tokyo, Japan
Tadanobu SATO
Disaster Prevention Research Institute, Kyoto University, Kyoto, Japan(Received for 4 Mar., 2002)

ABSTRACT
A basic study of the application of the Monte Carlo filter (MCF) to structural damage detection is reported. Inthat method, each probability distribution is expressed by many of its realizations, called particles or samples.
The advantage of the MCF is that it deals with non-linear and non-Gaussian problems. In terms of damage detection, non-Gaussian noise may be preferable because the damage tends to be concentrated on a specific part of astructure. Two kinds of numerical examples are shown. First, the stiffness of the experimental model for theshaking table test is identified by the MCF and EKF (Extended Kalman Filter). Based on hypothetical data,
numerical simulations of damage detection with non-Gaussian process noise then are performed and discussed.
Because the MCF results are given by many particles (samples), the detailed probabilistic nature of the identified
parameters also can be discussed.

1. INTRODUCTION
After a large earthquake, rapid detection of damage to structures,
especially important ones such as hospitals, bridges, fire stations, isneeded to prevent secondary disasters. Extensive time and costs are
involved in estimating damage to such structures visually, if
inspected in detail, the conventional way. Recently, many important structures are being monitored with sensors such as seismometers. A practical method for rapid detection of damage to the structures that uses monitored data is desirable from the standpoint ofdamage detection or health monitoring.
Linear/Nonlinear identification methodology has been wellstudied. Of the various reported methods, the Kalman filter is oneof the best known and widely used algorithm owing to its beautiful,sophisticated theory which requires only the first and second
moment of probabilistic nature because the filter is well establishedon the basis of linear Gaussian assumptions. In the field of structural identification, the Kalman filter has been one of the most widelyused tools (Yun & Shinozuka 1980, Hoshiya & Saito 1984).
Although various nonlinear identification methods have been proposed, most have the same limitation as the Kalman filter as they
are based on Gaussian noise, leading to the quadratic form of theobjective function.
Structural damage must be non-stationary phenomenon, therefore process noise has an important function because it works as aforgetting factor (Koh & See 1994, Hoshiya & Yoshida 1998, Sato& kaji 2001). Although Gaussian process noise is used in the
Kalman filter, other types of noise may be preferable for damage
detection when the nature of the damaged structure is considered.
Minor damage caused by seismic forces tends to be localized; e.g.,
it tends to be concentrated on a specific part of a structure. When alarge destructive earthquake hits, damage will occur in many partsof a structure simultaneously. Cautious inspections of structureswith such major damage is not required because the severity of thedamage is clear. In damage detection, the assessment of minor tomoderate damage is important. To detect such damage, Gaussian
process noise, which derives the quadratic form of objective function, may not be appropriate.
Non-Gaussian identification methods are being investigated(Kitagawa 1996). Basically, the non-Gaussian approach requires
much computation but because of remarkable advances in computer
performance, new approaches that use the Monte Carlo techniquehave become practical; e.g., the Genetic Algorithm. The Monte
Carlo Filter proposed by Kitagawa (1996) is a Monte Carlo methodthat deals with non-Gaussian noises. In it, probability distributions
are approximated by many of their realizations; i.e., the many particles or samples. The probabilistic nature of a state vector isexpressed by the many particles instead by the first and second
moment as in the Kalman filter. A method of damage detection thatuses the Monte Carlo filter is here proposed and demonstrated.
Formulation of identification by means of Monte Carlo filter first ispresented. It is the natural extension of the Kalman filter (linear
Gaussian) but essentially a different method in that it is not necessary to use Gaussian noise. Two kinds of numerical examples areshown; identification of the stiffness of the experimental model forthe shaking table test as reported by Loh (2000), performed by
MCF and ordinary EKF (Extended Kalman filter), and numericalsimulation of damage detection by the MCF with Gaussian or nonGaussian process noise based on hypothetical data.

KEY WORDS: damage detection, health monitoring, state space, non-Gaussian filter

102

I. YOSHIDA AND T. SATO

2.

METHODOLOGY

2.1 The Kalman and Monte Carlo Filters
The general form of the state space model, including the non-linear
state/observation equation and non-Gaussian noise, is described asfollows. The state equation is
xk = F (x k-1 , wk)(1)

where, xk is the state vector at the k-th step, and wk the process
noise. The observation equation isyk = H (xk , vk)(2)

where, yk is the observation vector at the k-th step, and vk the observation noise. When the state and observation equations are linear,and the noises Gaussian, the Kalman filter technique, simple, very
sophisticated filtering method, is applicable. The algorithm is presented in Fig. 1, where Q and R respectively represent the covariance matrix of wk and vk . The probability distribution of a statevector is described solely by the first and second moments. The
Kalman filter algorithm has two parts, a time-updating process andan observation-updating process. The time-updating process is a
one-step ahead prediction based on the information at the (k-1)-thstep. The predicted state vector and its covariance matrix aredenoted as xk / k-1 and Pk / k-1 . The observation-updating process consists of the estimation based on observation data yk at the k-th stepand information predicted by the time-updating process. The estimated state vector and its covariance matrix are given as xk / k andPk / k .
Several methods have been proposed for general, non-linear,and non-Gaussian problems, in which the probability function isapproximated by a group of step functions, linear segments, orGaussian distributions (Kitagawa 1996). These methods basicallyare limited to low dimensional problems because of the impracticality of extensive computation time when the dimension is large.

Fig. 1In the Monte Carlo filter technique, probability functions areapproximated by many of their realizations; i.e., particles or samples. The algorithm of MCF also is shown in Fig. 1. Its generalflow is similar to that of the Kalman filter, but in the MCF theprobabilistic nature of the state vector is described by many realizations instead of the first and second moments. The most important point of the MCF is the resampling in the observation updatingprocess. The likelihood of each sample realization with respect toobservation noise first is calculated. Next, the samples are resampled according to their likelihood ratios. Samples of large likelihood may be resampled several times, whereas samples of small
likelihood may not be resampled at all.
2.2

State Transfer and Observation Equations for Damage
Detection
The Monte Carlo filter algorithm is applied to structural damage
detection problems. The equation of motion is given by
Mur + Cur + Kur = -Mhz

(3)

where, M, C, and K represent the mass, damping, and stiffness
matrix respectively and ,and ur the relative acceleration,
velocity, and displacement to the input motion. The notation
represents input acceleration, and h is a vector whose components
are unity ( h = (1,1,1,... ,1)T ).
Two types of formulation are introduced. The first, proposedby Sato & Kaji (2000), is called Type-1 in this paper. The statevector, x, is composed simply of unknown parameters such as
stiffness and the damping ratio. The state equation is given by thesimple form.
xk = xk-1 + wk

(4)

The observation equation isyk = Hk xk + vk = -M -1Ck ur -M -1Kk ur -hz + vk

Time-updating and observation-updating processes of the filtering algorithms.
(a) Kalman filter; (b) Monte Carlo filter(5)

HEALTH MONITORING ALGORITHM BY THE MONTE CARLO FILTER BASED ON NON-GAUSSIAN NOISE

In this formulation, the observation vector, yk , represents the
response accelerations of all the nodes. Velocity and displacementof all the nodes also are necessary on the right side of Equation (5).
The equation of motion is considered in the observation equation,in which the state vector xk contains Ck and Kk . As both Equations(4)(5) are linear, the Kalman filter is applicable when the noises wk
and vk are Gaussian. The advantage of the Type-1 formulation isthat it is simple and practical, whereas the disadvantage is that itneeds motions of all degrees of freedom as observation data.
In the Type-2 formulation, the state vector x is composed ofthe motions of all the nodes and unknown parameters;ur( ) x = ( uu )

x
x= xa
b

a

(6)

r
r

where xb represents unknown parameters. The state equation isxa , k = F (xa, k-1, xb, k ) + wa, k(7a)xb, k = xb, k-1 + wb, k(7b)where F(x) is the function that represents one-step prediction. The
equation of motion is considered in the function, F(x). In the following numerical examples, the linear acceleration method wasused for F(x). The observation equation is written simply;

yk = Hxk + vk

(8)

where H is a constant simple matrix which indicates observation
point locations. In this formulation, identification with observed
data of part of the nodes can be performed theoretically. Only theType-2 formulation is discussed in the following numerical examples.
2.3 Type of Noises
In the state space model, process noise has an important role. The
probability density functions of Pearson’s family are shown in Fig.
2. The Cauchy, Student’s t, and Gausssian distributions are family
members. The Cauchy and Gausssian distributions have opposite
characteristics. The central part of the Gaussian distribution is the
thickest, whereas that of the Cauchy distribution is the thinnest,
indicative that the tail part of the Cauchy distribution is the heaviest. The distribution with heavier tail part than that of the Cauchy
distribution is proposed. It is a combination of the delta function
and uniform distribution. This distribution is illustrated in Fig.
2(b). This distribution here is called the U-D (Uniform-Delta) distribution.
Gaussian noise is interpreted as having L2-norm (Menke
1989) from its quadratic form of the objective function (penalty

Cauchyt ( f = 3)delta distributiont (f = 5)
Gaussuniform distribution

(a)

103

(b)

Fig. 2

Probability Density Functions. (a) Pearson’s family probability density functions, Cauchy, Student’s t(f : degree of freedom), Gaussian distribution; (b) Uniform-Delta (U-D) distribution

Fig. 3

Contour map shapes of probability density functions. (a) Norms of several types of probability density function; (b) Cauchy distribution; (c)Uniform and Delta distribution (U-D)

104

I. YOSHIDA AND T. SATO

function, cost functions). This means that distance is defined as inour physical world. Fig. 3 shows contour map shapes of three
types of probability density functions in two-dimensional space.
The Gaussian shape is a circle (ellipse), whereas that of the
Laplace distribution is a square, interpreted as having the L1-norm.
The shape of the Cauchy contour depends on the part of the probability density function. The shape of the contour map is almost acircle near the origin (mean), evidence that the norm is close to L2,
but at the tail part the norm seems to be lower. The shape of theU-D distribution is very simple, having only three levels of heightin two dimensional space. The highest part is at the mean point(origin), the second highest on the axes, and the other parts are thelowest. The shape of this probability density function also suggests a very low norm.
The type of probability density function determines the formof the objective function for identification and has a significant
effect on the estimated parameters. Identification made with low
norm noise tends to identify concentrated damage to a specific partof a structure as compared with identification made with higher
norm noises. In fact, however, minor to moderate damage causedby seismic forces tends to be localized, i.e., it tends to be concentrated. This suggests that lower norm noises are better for the purpose of damage detection, although Gaussian noises are widelyused. In the following numerical examples, Gaussian and U-Dnoises are compared.

steel structure is installed on the shaking table. The height of themodel is 650cm. Many sensors, acceleration meters, velocity
meters, and strain gauges are attached to the model. Shaking tests
were carried out under various conditions to obtain a set of earthquake response data. In this numerical study, data observed under
the following excitation conditions were used.
● Earthquake: Kobe (1995, meteorological observatory, verticalto fault component)● Scale for input motion: 48%
The type-2 formulation of the structural identification, given earlier was used for both the MCF and EKF methods. The formulation
for the EKF is essentially that proposed by Hoshiya & Saito(1984). The purpose of the identification is for estimation of thestiffness of the model, not for damage detection, therefore
Gaussian process noise also was used in the MCF. The particle

3. NUMERICAL EXAMPLES
3.1 Identification with experimental data
The experimental data of Loh (2000) was used to identify the stiffness of the model by two methods, Extended Kalman (EKF) andMonte Carlo (MCF) filters described earlier. Fig. 4 shows a sketchof the experimental model as well as the response accelerations on
the top floor and on the shaking table. A model of a five-story

Fig. 5

Fig. 4

Sketch of the experimental model and time histories of
response accelerations

Identified stiffnessess of the experimental model by the MCF (Monte Carlo filter) and EKF (Extended Kalman filter) methods

HEALTH MONITORING ALGORITHM BY THE MONTE CARLO FILTER BASED ON NON-GAUSSIAN NOISE(sample) number was 1000. In the EKF, no process noise wasconsidered. As no data on mass distribution was available, uniform mass distribution is assumed. Initial stiffness values were
determined based on the 1st predominant frequency of the Fourier
spectral ratio of response at the top floor to the motion of the shaking table.
The observed data for all the floors were used for the identification. The time histories of the identified stiffness are shown in
Fig. 5 as ratios to the initial value. Although the same stiffness
values are not identified by the two methods, the trend is the same.
The identified stiffness of the upper parts, especially the top floor,
is smaller than that of lower parts.
Identifications based on observed data for a specific floor also
were tried, but no stable solutions were obtained. The identified
stiffness varied depending on the observation point. When the
EKF was used, solutions sometimes were not obtained due tonumerical explosion. When the MCF was used, numerical explosions did not occur because the parameter ranges are defined, but
the solutions obtained also were not stable. Identifications madewith observation data for two specific floors also were tried. Morestable solutions were obtained, but solutions still sometimes
became unstable.
3.2 Damage detection with hypothetical data
In this numerical example, the hypothetical 5-DOF model, shownin Fig. 6, was used. All the elements (floors) have the same
weight, initial stiffness, and initial damping ratio; respectively
9.8tf/m3, 400 tf/m2, and 0.02. Damage is assumed to occur at element 3, where stiffness decreases from 400 to 360 tf/m2 and thedamping ratio increases from 0.02 to 0.04 during 6-7 seconds.
Due to this damage, the fundamental frequency of the model isreduced, from 1.0 to 0.99Hz.
The input and response motions for nodes 1 and 3 are shownin Fig. 7. The observation data obtained for damage detection add
3% (rms ratio) Gaussian white noises to the response data. The
number of particles is 1000. Identification for up to 10000 particles was tried, but the performance was almost the same. For
fewer than 1000 particles, however, the performance was bad.
When the SODF simulation was tried, good results were obtainedeven for a particle number of 100.

Damage detection by the MCF can be done under various
conditions. Results of the cases shown in Table 1 are here discussed. In all of them, the Cauchy distribution was used for theobservation noise. First, cases 1 and 2 results are discussed, inwhich the U-D distribution was used as the process noise. The
particle distributions of the estimated stiffness and damping ratioof element 3 are shown respectively in Figs. 8 and 9. These figures are considered to show the non-stationary probability density
function of the identified parameters. When data for all the nodes
was used, the stiffness distribution peak was clear as shown (Fig.
8). When the observation data was limited to nodes 1 and 3, thedistribution peak was not as clear (Fig. 9), but the trend for the
stiffness distribution peak to moves after about 6 seconds is stillclear. The trend of the damping ratio is observable, but generally
is not clear. Although many simulations under various conditions
were tried, identification of the change in the damping ratio wasvery difficult. Because the identified damping ratios are very
unstable, the discussion hereafter concentrates on stiffness.
Estimator needs to be determined from the particle distributions for practical damage detection. The most popular estimatoris the mean, which sometimes is inadequate. The histograms in
Fig. 10 represent the identified stiffness distributions of element 3
at 6 and 9 seconds in case 2, which corresponds to the time sections in Fig. 9. A single distribution peak is presented at 6 sec-

1
2
Fig. 7

Acceleration time histories of input motion and responses

3
4
5

Fig. 6

Hypothetical 5-DOF model with damage to element 3

105

Table 1. Numerical calculation cases

106

I. YOSHIDA AND T. SATO

onds, and the mean, median, and mode of distribution appear to bealmost the same. At 9 seconds, however, there are three peaks,
and the distribution shape is not symmetric. The mean value therefore may not be an appropriate estimator. Interestingly, the highest
peak is close to the true damage stiffness, and the second highest
one is close to the initial value. In this study, we used the highest
distribution peak, the mode of particles as the estimator.
The identified stiffnesses (highest distribution peak, mode) of
elements 2, 3, and 4 are shown in Fig. 11. When observation data
for all the nodes were used, damage was detected properly despite
the type of process noise. The identified stiffnesses of 2 and 4 arealmost constant, but the stiffness of 3 decreases after 6 seconds.
When the observed data were limited to nodes 1 and 3, U-D
process noise seemed to be better than Gaussian process noise inthis numerical simulation. These results suggest that the MCF
with non-Gaussian process noise is appropriate to use for damage

detection. The advantage using the U-D process noise could not
be definitively shown because the identified parameters are sensitive to numerical conditions. Further studies on the modeling/tuning of process noise and a more stable algorithm are needed.
4. CONCLUSIONSFormulation of damage detection by means of the Monte Carlo filter was introduced and the validity of that method demonstrated
for two types of numerical simulations. The conclusions are as
follows;
1) Formulation of damage detection by the Monte Carlo filter(MCF) is possible. The U-D distribution for which the norm
is very low, also can be used for damage detection.
2) Parameter identification in the experimental model can bemaded by the MCF and ordinary Extended Kalman filter

Stiffness

Stiffness
500.0

500.0

10.0

10.0
300.0

time

300.0

0.0

time
0.0

Damping Ratio

Damping Ratio

0.06

0.06

10.0
0.0
Fig. 8

10.0

time

0.0

0.0

Distribution of stiffness particles and the damping ratio of
element 3(case-1 observation data for all nodes, U-D Proc-ess
Noise)

Fig. 10

Fig. 9time
0.0

Distribution of stiffness particles and the damping ratio of
element 3(case-2 observation data for nodes 1 and 3, U-D process
noise)

Estimator of element 3 stiffness based on the particle distribution(case-2 observation data for nodes 1 and 3, U-D process noise)HEALTH MONITORING ALGORITHM BY THE MONTE CARLO FILTER BASED ON NON-GAUSSIAN NOISE

Fig. 11107

Identified stiffnesses of elements 2, 3, and 4. (a) U-D process noise; (b) Gaussian process noise(EKF) methods. Similar parameter trends are identified byboth methods.
3) Damage detection by the MCF is done with Gaussian or U-D
noise which has a very low norm in the objective function ofidentification. Numerical simulations findings suggest thatthe use of low norm process noise is suitable for damage
detection problems.
In the MCF, the parameters identified sometimes vary depending
on numerical conditions such as particle number, seed of random
number, and noise parameters. Further improvements such as
adaptive tuning, use of the proper type of process noise, and amore stable algorithm are needed for a practical tool. The MCF
appears to be a prospective tool to detect damage done to a structure, because it is applicable to a broad range of non-linear andnon-Gaussian problems.
REFERENCES
Kitagawa, G.., 1996. Monte Carlo Filter and smoother for non-Gaussianstate space models, Journal of Computational and Graphical Statistics,
5(1), 1-25.
Koh, C.G. and See, L.M., 1994. Identification and uncertainty estimation of

structural parameters, ASCE Journal of Engineering Mechanics,
120(6), 1219-1236.
Loh, C.-H., 2000. The report o f shaking table test results and numericalsimulation o f a five- story steel frame structure, NCREE Report, (inChinese).
Menke, W., 1989. Geophysical data analysis: discrete inverse theory,
Academic Press, Inc.
Hoshiya, M. and Saito, E., 1984. Structural identification by extended
Kalman Filter, ASCE Journal of Engineering Mechanics, 110(12),1757-1770.
Hoshiya, M. and Yoshida, I., 1998. Process noise and optimum observationin conditional stochastic fields, ASCE Journal of Engineering
Mechanics, 124(12), 1325-1330
Tadanobu SATO, Kai QI. 1998. Adaptive H∞ filter: Its application to
structural identification, ASCE Journal of Engineering Mechanics,
124(11), 1233-1240
Sato, T. and Kaji, K., 2001. Adaptive Monte Carlo filter and structure identification, Proceedings of the first International Conference on Monte
Carlo Simulation, pp.441〜pp.447
Yun, C.-B. and Shinozuka, M., 1980. Identification of nonlinear structuraldynamic systems, ASCE Journal of Engineering Mechanics, 8(2),187-203