- Email: [email protected]
Sensor validation using minimum mean square error estimation
ARTICLE IN PRESS Mechanical Systems and Signal Processing 24 (2010) 1444–1457
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp
Sensor validation using minimum mean square error estimation Jyrki Kullaa a,b, a b
Department of Applied Mechanics, Helsinki University of Technology, Finland Department of Mechanical and Production Engineering, Helsinki Metropolia University of Applied Sciences, P.O. Box 4021, FIN-00079 Metropolia, Finland
a r t i c l e i n f o
abstract
Article history: Received 14 February 2009 Received in revised form 26 November 2009 Accepted 8 December 2009 Available online 21 December 2009
Sensor fault can be detected and corrected in a multichannel measurement system with enough redundancy using solely the measurement data. A single or multiple sensors can be estimated from the remaining sensors if training data from the functioning sensor network are available. The method is based on the minimum mean square error (MMSE) estimation, which is applied to the time history data, e.g. accelerations. The faulty sensor can be identified and replaced with the estimated sensor. Both spatial and temporal correlation of the sensors can be utilized. Using the temporal correlation is justified if the number of active structural modes is larger than the number of sensors. The disadvantages of the temporal model are discussed. Experimental multichannel vibration measurements are used to verify the proposed method. Different, and also simultaneous, sensor faults are studied. The effects of environmental variability and structural damage are discussed. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Sensor fault Minimum mean square error estimation Structural health monitoring
1. Introduction Structures equipped with sensors are becoming common due to the development of the sensor and sensor network technology as well as different applications exploiting sensor information. For example, vibration-based structural health monitoring is heavily based on the measurement data recorded during a long period. It is therefore most important that the sensors are working properly and possible sensor faults are detected in order to maintain the reliability of the system. Otherwise a sensor fault can be misinterpreted as a structural damage. Monitoring systems typically include several sensors at different locations of the structure in order to extract features for damage detection or damage localization. Such sensor network constitutes a redundant system, which can be used to detect sensor malfunction or failure, identify, and even correct the faulty sensor. These topics are discussed in this study. Dunia et al. [1] studied detection, isolation, and reconstruction of a faulty sensor using principal component analysis (PCA). They reconstructed each sensor in turn using the remaining sensors, and calculated a sensor validity index from residuals before and after reconstruction for fault detection and isolation. They proposed two approaches to reconstruct a single sensor, an iterative method and optimization, both resulting in the same closed-form solution. Different types of fault, bias, complete failure, drifting, and precision degradation were studied. Different residuals were investigated. They proposed a sensor validity index (SVI) to identify the faulty sensor by replacing one sensor in turn with the estimated sensor. The minimum SVI was found in the case where the faulty sensor was replaced with the estimated one. Dunia and
Corresponding author at: Department of Mechanical and Production Engineering, Helsinki Metropolia University of Applied Sciences, P.O. Box 4021, FIN-00079 Metropolia, Finland, Tel.: + 358 207836176; fax: +358 207836102. E-mail address: jyrki.kullaa@metropolia.fi
0888-3270/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2009.12.001
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1445
Qin [2] also derived conditions for fault detectability, reconstructability, and identifiability. Both sensor and process faults could be identified, but the fault direction must be known, which is easy for a sensor fault, but difficult for process faults. Kerschen et al. [3,4] also applied PCA to sensor validation and used the angle between the principal subspaces as a feature for fault detection. The principal component analysis was made separately for each time record. They also derived a closed-form equation to estimate a single sensor. The faulty sensor was identified by removing one sensor in turn. The faulty sensor was the removed sensor in the case with the minimum angle. Friswell and Inman [5,6] introduced a sensor validation method based on the modal model of the structure. The advantage is that the sensors are not assumed to be functioning initially but instead a modal model has to be available. Two approaches were studied, modal filtering and PCA. They concluded that the modal filtering approach performs better if an accurate modal model is available. They also found that a multiplicative sensor fault is more difficult to locate than an additive fault. Abdelghani and Friswell [7,8] studied model-based methods to validate sensors with additive or multiplicative faults. Two residual generation techniques were studied with an additive fault, the parity space and the modal filtering approaches. For the multiplicative fault, they introduced a correlation index to isolate the faulty sensor. The index was computed from proposed new residuals. Hernandez-Garcia and Masri [9] compared PCA, independent component analysis (ICA), and modified ICA (MICA) for sensor fault detection. They used two features for fault detection: the Hotelling’s T2 statistic and the squared prediction error (SPE). They concluded using a numerical model that the detection performance of ICA and MICA was higher than that of PCA. A disadvantage of PCA is that the number of principal components must be determined, which is not always straightforward with noisy data. It has also been reported that the estimation of several variables simultaneously is complex. This is needed if more than one sensor is faulty or if temporal correlation is also applied. Kullaa [10] introduced factor analysis for sensor validation and compared it with PCA. Similar to PCA, factor analysis also has one parameter, the number of factors, which has to be chosen. However, the sensitivity to the parameter was seen to be lower than in PCA. Also, a closed-form solution to reconstruct several sensors using the remaining sensors was derived. Kullaa [11] introduced the minimum mean square error (MMSE) estimation to sensor validation. It works directly in the data space without a need to determine any model parameters, but with a cost of a slower performance. Also, a spatiotemporal extension was introduced [12], which can be applied if the number of sensors is lower than the number of active modes in the structure. The estimation error was also derived but not applied in the study. The faulty sensor was identified by removing one sensor in turn and performing MMSE estimation to the remaining sensors. The faulty sensor was the missing sensor in the analysis with the lowest mean-square residual compared to that of the training data. The main disadvantage is that the analysis has to be made for each sensor resulting in a slow performance. In this paper, a considerably faster algorithm for sensor validation is derived based on the MMSE estimation. The performance of both the sensor estimation and faulty sensor identification are increased. Also, the estimation error is utilized in the faulty sensor identification, which results in a higher reliability compared to the previously proposed method [11]. In most of the aforementioned studies it was assumed that the number of sensors is higher than the number of excited modes. If the number of active modes is larger than the number of sensors, there is a lack of redundancy between sensors. The time histories could then be filtered to a limited bandwidth. Another solution is to use temporal correlation together with spatial correlation. Spatial correlation contains information between the sensors at the same time instant (a snapshot of the structural motion), whereas temporal correlation contains frequency information. Bearing this in mind, it would be advantageous to use spatial correlation only to distinguish between structural damage and sensor fault. Natural frequencies have been observed to be more sensitive to structural damage or environmental variations than the mode shapes. These effects are also studied in this paper. The linear MMSE estimator model is derived in the next section followed by the spatiotemporal correlation model. Algorithms for sensor fault detection, identification, and correction are proposed. The method uses solely the measurement data, no numerical model of the structure is needed. The model is then applied to an experimental study of a wooden bridge. It is first assumed that one sensor only is faulty but later simultaneous sensor faults are also investigated. Different fault types, bias, drifting, precision degradation, and gain are studied. Additive and multiplicative faults result in a change of different statistics. Also, the effects of the environment and minor damage are studied. The spatiotemporal correlation model is compared to the purely spatial correlation model. Finally, conclusions and practical suggestions from the research are given. 2. Minimum mean square error (MMSE) estimator Sensor validation is approached by estimating one or more sensors using the remaining ones. The likelihood ratio test is then applied to detect sensor fault. It is assumed that training data have been acquired from the sensor network, in which all sensors are functioning. It is also assumed that the number of active structural modes is less than the number of sensors. This is a necessary condition for enough redundancy. Later, this condition can be relaxed if temporal correlation is also taken into account. In multivariable measurements with enough redundancy, a subset of observation can be estimated using the remaining variables. Each observation is divided into observed variables v and missing variables u: u x¼ ð1Þ v
ARTICLE IN PRESS 1446
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
with a partitioned covariance matrix (estimated using the training data): " #
Ruu Rvu
R¼
Ruv Rvv
ð2Þ
The objective is to estimate the missing variables u using the observed variables v. There are two cases in which the estimator of u is linear [13]: 1. the data x are multinormally distributed, 2. the estimator u^ is constrained to be a linear function of the observed data: u^ ¼ Kv
ð3Þ
where K is an unknown matrix, and u and v are assumed to be random vectors with zero means and known covariances, but the joint distribution need not be Gaussian. Both cases result in identical estimators [13]. Here we consider case 2 only, as the assumption of a Gaussian distribution is not generally valid, e.g. for harmonic vibration. The optimal linear estimator u^ is obtained by minimizing the mean-square error (MSE) eMS ^ T ðuuÞÞ ^ eMS ¼ EððuuÞ
ð4Þ
where E( ) is an expectation operator. The mean-square error is minimized by choosing [13]: K ¼ Ruv R1 vv
ð5Þ
and finally EðujvÞ u^ ¼ Kv ¼ Ruv R1 vv v
ð6Þ
where u|v denotes u given v. The error covariance matrix is [13]: ^ ^ T Þ ¼ Ruu Ruv R1 covðujvÞ EððuuÞðu uÞ vv Rvu
ð7Þ
The main disadvantage of using Eq. (6) for sensor fault detection and identification is that when each variable is estimated in turn, a large matrix Rvv has to be repeatedly inverted. The dimension of the matrix is (p 1) (p 1), where p is the number of sensors. This makes the algorithm slow with a large sensor network. In the following, a faster algorithm is derived, which involves inverting a matrix with a low dimension (usually one!) for each sensor. The inverse of the full covariance matrix is nevertheless needed but it is computed only once. 2.1. Linear estimator Modified formulas for the linear estimator (Eqs. (6) and (7)) can be obtained using the precision matrix C, which is the inverse of the covariance matrix R (C = R 1). Using the matrix inversion in block form, the partitioned covariance matrix in terms of the partitioned precision matrix is " # " #1
Ruu Rvu "
¼
Ruv Cuu ¼ Rvv Cvu
Cuv Cvv
1 ðCuu Cuv C1 vv Cvu Þ
1 ðCuu Cuv C1 Cuv C1 vv Cvu Þ vv
1 1 C1 vv Cvu ðCuu Cuv Cvv Cvu Þ
1 1 1 C1 Cuv C1 vv þ Cvv Cvu ðCuu Cuv Cvv Cvu Þ vv
#
ð8Þ
Applying the matrix inversion lemma, and after some manipulation, the necessary submatrices can be shown to be 1 1 1 Ruu ¼ C1 Cvu C1 uu þ Cuu Cuv ðCvv Cvu Cuu Cuv Þ uu
ð9Þ
1 R1 vv ¼ Cvv Cvu Cuu Cuv
ð10Þ
1 1 Ruv ¼ C1 uu Cuv ðCvv Cvu Cuu Cuv Þ
ð11Þ
1 Rvu ¼ RTuv ¼ ðCvv Cvu C1 Cvu C1 uu Cuv Þ uu
ð12Þ
Substituting Eqs. (9)–(12) into Eqs. (6) and (7) results in 1 EðujvÞ ¼ Ruv R1 vv v ¼ Cuu Cuv v
ð13Þ
1 covðujvÞ ¼ Ruu Ruv R1 vv Rvu ¼ Cuu
ð14Þ
and
The dimension of Cuu is typically much lower than that of Rvv resulting in a faster algorithm. The precision matrix C is first computed by inverting the full covariance matrix R, but it is made only once.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1447
With enough redundancy in the sensor network, ujv can be assumed normally distributed pðujvÞ ¼ NðKv; WÞ
ð15Þ
where K ¼ C1 uu Cuv
W ¼ C1 uu
ð16Þ
The application of the MMSE estimation to sensor fault detection, identification, and reconstruction is discussed in the following sections. 3. Spatiotemporal correlation The variables in the sensor fault identification are the sequential sensor values from each sensor. If the process can be assumed stationary with zero mean, temporal correlation can also be utilized. The time-shifted covariance matrix estimate, with a time shift i, is computed by Ri ¼
N i 1 X xðkÞT xðkþ iÞ Ni k ¼ 1
ð17Þ
where x(k)= x(kDt) is the kth sample at a time instant kDt, where Dt is the sampling period and N is the number of samples. If the sensor variables are iT h ð18Þ x ¼ xðkmÞT ; . . . ; xðk1ÞT ; xðkÞT ; xðkþ 1ÞT ; . . . ; xðkþ nÞT then the covariance matrix 2 R0 R1 6 R0 6 R1 R¼6 6 ^ ^ 4 Rmn Rmn þ 1
estimate R is 3 2 R0 Rm þ n 7 6 T Rm þ n1 7 6 R1 7¼6 7 6 ^ & ^ 5 4 RTm þ n R0
R1
R0
^
&
RTm þ n1
Rm þ n
3
Rm þ n1 7 7 7 7 ^ 5 R0
ð19Þ
where m+n is the model order. Hence, instead of the time shifts m and n, one parameter only, the model order, has to be defined. The last equality in Eq. (19) results from the fact that Ri ¼ RTi If the sensor to be investigated is sensor j, then the missing variables are h iT u ¼ xj ðkmÞT ; . . . ; xj ðk1ÞT ; xj ðkÞT ; xj ðk þ 1ÞT ; . . . ; xj ðkþ nÞT
ð20Þ
ð21Þ
corresponding to row j and column j of each Ri. All these missing values are estimated simultaneously. As a result, the number of each estimated variable is equal to m+n + 1 except for the beginning and end of the time histories. A more accurate estimate is obtained by averaging. It should be noted that appending two time histories from different measurements and using a model order larger than zero, an erroneous estimation will occur around the time discontinuity. There are several disadvantages in using the spatiotemporal model instead of the spatial model: (1) The process has to be stationary; (2) the model order has to be determined; (3) the time discontinuities introduce estimation error; (4) an increased number of variables results in a slower performance; and (5) separation of structural fault from a sensor fault may be difficult. Therefore, spatiotemporal correlation should only be used if the number of sensors is inadequate for required redundancy. In that case, filtering can be used as an alternative to decrease the active modes in the data. 4. Sensor fault detection, identification, and correction This section discusses the application of the presented theory to sensor validation: fault detection, faulty sensor identification, and sensor reconstruction. 4.1. Fault detection The first step in the sensor validation is fault detection. The MMSE model is built using the training data from the functioning sensor network. This model is the null hypothesis H0, and the hypothesis test is H0 :
EðujvÞ ¼ K0 v; covðujvÞ ¼ W0
H1 :
EðujvÞaK0 v; covðujvÞaW0
ð22Þ
ARTICLE IN PRESS 1448
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
where K0 and W0 are estimated from the training data. The generalized likelihood ratio test (GLRT) [14,15] decides H1, if LG ðujvÞ ¼
X
ln
pðujv; H1 Þ 4g pðujv; H0 Þ
ð23Þ
where g is some threshold and assuming u|v is normally distributed, pðujv; H1 Þ ¼ ln pðujv; H1 Þln pðujv; H0 Þ pðujv; H0 Þ 1 jW1 j 1 1 ðuK1 vÞT W1 ðuK0 vÞT W1 ¼ ln 1 ðuK1 vÞ þ 0 ðuK0 vÞ 2 jW0 j 2 2
ln
ð24Þ
where K1 and W1 are estimated separately for each measurement. The feature used for fault detection is the log-ratio given by Eq. (24). The number of variables is equal to the number of sensors. This number is usually too large for a statistical analysis and a dimensionality reduction is made using the principal component analysis (PCA) [16]. The largest principal components point at the directions of the largest variance in the multidimensional space and they are finally used to design a control chart [17] for sensor fault detection. For different types of sensor faults, a suitable fault indicator depends on whether the fault affects the mean or the variance. Shewhart chart: The Shewhart chart, or xbar chart [17], monitors the change of the mean value, and the plotted variable is the subgroup mean x. The upper and lower control limits, UCL and LCL, respectively, are: ) UCL ð25Þ ¼ x 7 A3 S LCL where x ¼ meanðxÞ is the average of the subgroup means, S is the average of the subgroup standard deviations, and A3 is a constant depending on the subgroup size. Here, the subgroup size is 20, and A3 is 0.6797. S Chart: The S chart monitors the variability [17] and the plotted variable is the sample standard deviation S. The upper and lower control limits, UCL and LCL, respectively, are UCL ¼ B4 S LCL ¼ B3 S
ð26Þ
where B3 and B4 depend on the subgroup size. With a subgroup size of 20, B3 =0.5102 and B4 = 1.4898. Shewhart T chart: The multivariate counterpart of the Shewhart chart is the Shewhart T control chart [17], where the plotted characteristic is: T 2 ¼ nðxx ÞT S1 ðxx Þ
ð27Þ
where x is the subgroup average, x is the process average, which is the mean of the subgroup averages when the process is in control, and S is the matrix consisting of the grand average of the subgroup variances and covariances. The upper control limit is UCL ¼
pðmþ 1Þðn1Þ Fa;p;mnmp þ 1 mnmp þ 1
ð28Þ
where p is the dimension of the variable, n is the subgroup size, m is the number of subgroups, when the process is assumed to be in control, and Fa,p,mn m p + 1 denotes the a percentage point of the F distribution with p and mnmp+1 degrees-offreedom. The MMSE estimation is applied to sensor fault detection using the following procedure. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.
Define the training data with no sensor faults for hypothesis H0. Estimate mean l and covariance matrix R of the data. Compute the precision matrix: C = R 1. Select sensor and set it as a missing variable u. Form matrices lu, lv, Cuu, Cvv, Cuv, and Cvu by partitioning l and C. Compute K and C using Eq. (16). Return to 4 until all sensors have been modelled. Define the test data (typically one time series) for hypothesis H1. Repeat steps 2–7 for the test data. Compute the log-ratio (24) for each sample. Define the in-control samples with no sensor faults. Perform a dimensionality reduction using the principal component analysis. Perform a statistical analysis using control charts. If the control chart signals, it is an indication of sensor fault. Move on to fault isolation and reconstruction.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1449
4.2. Faulty sensor identification Once a sensor fault is detected, the faulty sensor must be identified. It is proposed that the largest log-likelihood ratio reveals the faulty sensor. The MMSE estimation is applied to faulty sensor identification using the following procedure. 1. Compute the log-likelihood ratio for each sensor using Eq. (23). 2. Find the largest log-likelihood ratio, which indicates the faulty sensor. 4.3. Sensor correction The sensor correction can be made once the faulty sensor has been isolated. The model is built using the training data from all sensors without fault. The faulty sensor is reconstructed by removing the sensor from the network and estimating it using the remaining sensors. The procedure is as follows. 1. 2. 3. 4.
Set the faulty sensor as the missing variable u. Compute the mean and variance of u|v for each sample using Eqs. (13) and (14), respectively. EðujvÞ is the reconstructed sensor. covðujvÞ can be used to plot the confidence limits.
5. Experimental results Sensor faults were investigated with a monitoring system built in the laboratory. The structure was a wooden bridge shown in Fig. 1 with a total mass of 36 kg. A random excitation was generated with an electrodynamic shaker to activate the vertical, transverse, and torsional modes. Fifteen Kistler 8712A5M1 piezo-electric accelerometers were used for the response measurements. They were located at three different longitudinal positions, each containing five accelerometers: vertical and transverse accelerations were measured at both top flanges but transverse acceleration at one of the bottom flanges only. The sampling frequency was 256 Hz and the measurement period was 32 s. The modal parameters were identified from the output-only data, which resulted in 20 natural modes below 128 Hz. As this number was larger than the number of sensors, there was lack of redundancy. Therefore, the signals were low-pass filtered below 64 Hz and re-sampled. In this frequency range, the structure contained 14 natural modes. The data used for sensor validation were three subsequent measurements and a fourth measurement acquired later on the same day. These four records were appended after each other. Different fault types in different sensors were studied: bias, drifting, precision degradation, and gain. The sensor fault was introduced to the fourth measurement. The first two measurements were used to build the model and to design the control charts.
Fig. 1. Wooden bridge.
ARTICLE IN PRESS 1450
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
In sensor fault detection, the MMSE estimation resulted in log-ratios (24) for each sensor. The dimensionality of the feature vector was thus 15, and further dimensionality reduction was made using PCA. The largest principal component was then used in sensor fault detection using control charts. Once a fault was detected, further procedures were performed to identify and correct the faulty sensor. Also, the effects of the environmental variability and structural damage (point masses) on the sensor validation were studied. 5.1. Gain fault A gain fault was studied multiplying sensor 14 by a constant 1.3. Because the mean acceleration is zero, the gain fault affects the change in the variance and the S chart is a proper fault indicator. The S control chart for the largest principal component of the log-ratios (24) is shown in Fig. 2a to detect a sensor fault. Notice the logarithmic scale. The subgroup size was 20. The sensor fault was clearly detected. The faulty sensor identification based on the maximum likelihood ratio is shown in Fig. 3, plotting the log-likelihood ratio for each sensor. The faulty sensor was correctly identified by the proposed method. After identifying the faulty sensor, it was removed and the remaining sensors were re-analyzed. The resulting control chart in Fig. 2b is considerably different indicating no clear fault. It can be concluded that the fault was localized to sensor 14. 5.2. Spatiotemporal analysis Spatiotemporal correlation model was studied with sensor data in which there was not enough redundancy. The sampling frequency was 256 Hz and no low-pass filtering was applied. The number of active modes was thus higher than the number of sensors.
S
101
100 100
200
300
400
500
600
700
800
600
700
800
Sample Number
S
101
100 100
200
300
400
500
Sample Number Fig. 2. S control charts for sensor fault detection. Gain fault in sensor 14: (a) all sensors included and (b) sensor 14 removed.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1451
Faulty Sensor Identified: 14 1000
Log−Likelihood Ratio
800
600
400
200
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number Fig. 3. Faulty sensor identification. Gain fault in sensor 14.
102 102
Tave
Tave
100
100
10−2
10−4 10−2 500
1000 Sample Number
500
1500
1000 Sample Number
1500
Tave
102
100
10−2
500
1000 Sample Number
1500
Fig. 4. Shewhart T control charts for sensor fault detection. Bias fault in sensor 1: (a) model order 0 (spatial correlation); (b) model order 5 and (c) model order 5 with sensor 1 removed.
ARTICLE IN PRESS 1452
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
Faulty Sensor Identified: 1
Log−Likelihood Ratio
250
200
150
100
50
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number x 105
Faulty Sensor Identified: 1
8
Log−Likelihood Ratio
7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number Fig. 5. Faulty sensor identification. Bias fault in sensor 1: (a) model order 0 (spatial correlation) and (b) model order 5.
Bias fault was simulated by adding a constant value 0.01 to sensor 1. This resulted in a change in the mean and the Shewhart T control chart for the largest principal component of the log-ratios was applied to detect the sensor failure. With spatial correlation (model order= 0), the sensor fault was detected (Fig. 4a), but using also the temporal correlation with a model order 5, the fault was more clearly detected (Fig. 4b, notice the logarithmic scale). In either case, the faulty sensor was correctly identified from the maximum log-likelihood (Fig. 5). Also, the reconstruction was more accurate using the spatiotemporal model. Fig. 6 shows the true and estimated sensor (top, a detail), their difference (bottom), and the difference between the measured (faulty) and estimated sensor (middle) from both the spatial and spatiotemporal analyses. The sensor estimated using the spatiotemporal model agrees well with the true sensor, whereas the estimation error is relatively high using the spatial model. After identifying the faulty sensor, it was removed and the remaining data were re-analyzed. The resulting control chart is shown in Fig. 4c, which now indicates no fault. It can be thus concluded that the fault was successfully localized to sensor 1.
5.3. Effects of environment and damage A relevant question in structural health monitoring is if sensor fault can be distinguished from structural damage or environmental or operational variability. In the cases above, different sensor faults were detected from the data recorded at similar environmental and structural conditions as the training data. In this section, the fourth measurement is taken from a different environmental condition than the training data. Also, minor structural change was present with an added mass of size 0.197 kg on the structure with a total mass of 36 kg. The
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1453
Reconstructed sensor: 1 0.1 0 −0.1 0.1 0 −0.1 1
2
3
4
1
2
3 Measurement
4
0.1 0 −0.1
Reconstructed sensor: 1 0.1 0 −0.1 0.1 0 −0.1 1
2
3
4
1
2
3 Measurement
4
0.1 0 −0.1
Fig. 6. Sensor validation with bias fault of sensor 1. Model order is (a) 0 (spatial correlation) and (b) 5. Detail of the true (blue) and estimated (green) sensor (top); residual: measured – estimated (middle); residual: true – estimated (bottom). The red lines are the estimated 3s error bounds. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
effects of these changes on the modal properties have been found to be significant [18]. In the present study, the objective is to separate a sensor fault from other changes in the system. In this study, sensor 10 suffered from drifting (a linear change in the mean from 0 to 0.01 within measurement). The fault was detected and the faulty sensor correctly identified using a spatial correlation model. The control chart is shown in Fig. 7a showing a clear change in the system. Sensor 10 was then removed resulting in the control chart shown in Fig. 7b and suggesting sensor 9 is faulty. Also, removing sensor 9 resulted in the control chart shown in Fig. 7c. Now, looking at all three charts it can be seen that removing sensor 10 had a large influence causing the drifting effect to vanish, whereas a further removal of sensor 9 did not make any enhancement. Therefore, it can be concluded that sensor 10 was faulty but there were also other influences present.
5.4. Simultaneous sensor faults Multiple sensor faults were studied introducing different faults to three sensors simultaneously: precision degradation (s = 0.01) to sensor 2, drifting (linear change of mean from 0 to 0.1) to sensor 5, and gain (multiplication by 1.3) to sensor 14. The data were filtered below 64 Hz. The sensor fault was easily detected, but the main interest here was the identification of the faulty sensors. Increasing the estimated number of faulty sensors, the possible combinations to be investigated increased. For example, with 15 sensors, estimating 1, 2, 3, or 4 faulty sensors, the number of combinations was 15, 105, 455, and 1365, respectively. The results of faulty sensor identification with different number of estimated faulty sensors are shown in Fig. 8, each bar representing a different sensor combination. If one sensor only was assumed faulty, sensor 5 was found (Fig. 8a). With two faulty sensors estimated, sensors 2 and 5 were found (Fig. 8b). With the correct number of three faulty sensors, all faulty sensors were correctly identified as faulty (Fig. 8c). If four faulty sensors were searched for, sensors 2, 5, 13, and 14 were
ARTICLE IN PRESS 1454
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
0
0
−20
−20
−40
−40 xave
xave
−60 −80 −100
−60 −80
−120
−100
−140 −120
−160 200
400 600 Sample Number
200
800
400 Sample Number
600
800
0
xave
−20
−40
−60
−80 100
200
300 400 500 Sample Number
600
700
800
Fig. 7. Shewhart control charts for sensor validation with measurement 4 at different environment and with a minor structural damage. Drifting of sensor 10: (a) all sensors present; (b) sensor 10 removed and (c) sensors 9 and 10 removed.
identified as faulty (not shown). Using spatiotemporal correlation with model order 4 also resulted in a correct identification with a higher sensitivity, but at a cost of a higher computational time (Fig. 8d). The reconstructed and true sensors 2, 5, and 14 are shown in Fig. 9 using spatial correlation. It can be seen that sensors 2 and 5 were accurately reconstructed, whereas the estimation error of sensor 14 was higher. This was also anticipated as can be seen from the estimated 3s error bounds in Fig. 9b showing a detail of the sensor estimation. It is interesting that if one sensor only was assumed faulty, the residuals of the other two faulty sensors 2 or 14 showed no remarkably high values (Fig. 8a) and thus could not be identified. Therefore, the correct number of faulty sensors has to be known for a reliable identification of all faulty sensors. Estimating the correct number of faulty sensors remains a subject for a further study.
6. Conclusion Sensor fault detection, identification, and correction were studied using the minimum mean square error (MMSE) estimation with spatial and spatiotemporal correlation between the variables. The proposed method was validated using experimental vibration data. The following observations were made. Different statistics were used to detect additive and multiplicative sensor faults. Additive faults (drifting and bias) can be detected from the change in the mean of the log-ratios, whereas multiplicative faults (gain and precision degradation) are detected from the change in the variance of the log-ratios. Therefore, two types of control charts were used: the univariate or multivariate Shewhart chart to monitor the mean, and the univariate S chart to monitor the variance. In practice, it is not known in advance what type of sensor fault is present, and therefore both charts should be constructed.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
x 107
Faulty Sensor Identified: 5
3
3
2.5
2.5
Log−Likelihood Ratio
Log−Likelihood Ratio
x 107
2 1.5 1 0.5
1455
Faulty Sensors Identified: 2 5
2 1.5 1 0.5
0
0 0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number x 10
7
20
40
60
80
100
Sensor Combination
Faulty Sensors Identified: 2 5 14
x
107
Faulty Sensors Identified: 2 5 14
3 8 Log−Likelihood Ratio
Log−Likelihood Ratio
2.5 2 1.5 1
6
4
2 0.5 0
0 0
100
200
300
Sensor Combination
400
0
100
200 300 Sensor Combination
400
Fig. 8. Faulty sensor identification with different numbers of estimated faulty sensors: (a) 1, (b) 2 (c) 3 using spatial correlation and (d) 3 using spatiotemporal correlation with model order 4. Each bar represents one sensor combination. Sensors 2, 5, and 14 are faulty.
A faster algorithm for the sensor estimation was developed compared to the previously introduced method. The results are identical, but the proposed method requires that the full covariance matrix be inverted only once and the computational performance for each sensor is high. In the previous method, each sensor estimation required an inverse of a large matrix. The estimation resulted in the expected value of each sensor as well as the error covariance. The expected value can be used as the reconstruction of a faulty sensor. The error covariance can be used to assess the confidence limits. Because the origin of the sensor data in this study was structural vibration, the sensor data also had temporal correlation. This additional correlation can be utilized if the number of active modes is larger than the number of sensors, for which the spatial correlation model may result in poor estimation. On the other hand, if there is enough redundancy, the spatial correlation model should be preferred, because the spatiotemporal model includes several disadvantages or assumptions that may be violated. The proposed method was relatively easy to expand to simultaneous sensor faults. The method could find all faulty sensors if the number of faulty sensors was known. With an estimated lower number of faulty sensors, all sensors that were identified as faulty were indeed faulty. The computational effort increases considerably with an increasing number of faulty sensors. Therefore, there is a practical upper limit of faulty sensors that can be found in a large sensor network. Also, if several sensors are assumed faulty, the remaining number of sensors used for estimation decreases, resulting in a higher estimation error.
ARTICLE IN PRESS 1456
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
Reconstructed sensor: 2 5 14 0.05 0 −0.05 1
2
3
4
1
2
3
4
1
2
3 Measurement
4
0.05 0 −0.05 0.1 0 −0.1
Reconstructed sensor: 2 5 14 0.05 0 −0.05 0.05 0 −0.05 0.1 0 −0.1 Time Fig. 9. Reconstructed (red) and true (blue) sensors 2 (top), 5 (middle), and 14 (bottom): (a) all measurements; (b) detail with estimated 3s error bounds (black dotted line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
A method was proposed to estimate if the detected fault was due to a sensor fault or another effect. The identified faulty sensor was removed from the set and the remaining sensors were re-analyzed. If the control chart showed a significant change, the removed sensor was probably faulty. On the other hand, if the control chart did not change much, the primary source of change was elsewhere (global effect). It should be noted, however, that the effects of environment or minor damage could not be removed. This problem remains a subject for a future study. It is anticipated that in the future more structures are equipped with a sensor network having numerous sensors. This fact makes the proposed method valuable and available for practical applications.
Acknowledgements This research was performed in the Multidisciplinary Institute of Digitalisation and Energy (MIDE) research programme launched at Helsinki University of Technology to celebrate its centennial as a university. References [1] R. Dunia, S.J. Qin, T.F. Edgar, T.J. McAvoy, Identification of faulty sensors using principal component analysis, AiChe Journal 42 (10) (1996) 2797–2812. [2] R. Dunia, S.J. Qin, Joint diagnosis of process and sensor faults using principal component analysis, Control Engineering Practice 6 (1998) 457–469. [3] G. Kerschen, P. De Boe, J.-C. Golinval, K. Worden, Sensor validation for on-line vibration monitoring, in: Proceedings of the Second European Workshop on Structural Health Monitoring, Munich, Germany, 2004, pp. 819–827. [4] G. Kerschen, P. De Boe, J.-C. Golinval, K. Worden, Sensor validation using principal component analysis, Smart Materials and Structures 14 (2005) 36–42.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1457
[5] M.I. Friswell, D.J. Inman, Sensor validation of smart structures, in: Proceedings of IMAC-XVIII: A Conference on Structural Dynamics, San Antonio, Texas, 2000, pp. 483–489. [6] M.I. Friswell, D.J. Inman, Sensor validation of smart structures, Journal of Intelligent Material Systems and Structures 10 (1999) 973–982. [7] M. Abdelghani, M.I. Friswell, Sensor validation for structural systems with additive sensor faults, Structural Health Monitoring 3 (3) (2004) 265–275. [8] M. Abdelghani, M.I. Friswell, Sensor validation for structural systems with multiplicative sensor faults, Mechanical Systems and Signal Processing 21 (2007) 270–279. [9] M.R. Hernandez-Garcia, S.F. Masri, Multivariate statistical analysis for detection and identification of faulty sensors using latent variable methods, Advances in Science and Technology 56 (2008) 501–507. [10] J. Kullaa, Three models for sensor validation, in: Proceedings of the 7th International Workshop on Structural Health Monitoring, Stanford, CA, 2009, pp. 529–536. [11] J. Kullaa, Sensor fault identification and correction in structural health monitoring, in: Proceedings of ISMA2006, International Conference on Noise and Vibration Engineering, Leuven, Belgium, 2006, pp. 873–884. [12] J. Kullaa, Sensor fault identification and correction in vibration-based multichannel structural health monitoring, in: Proceedings of the 6th International Workshop on Structural Health Monitoring, Stanford, CA, 2007, pp. 606–613. [13] H.W. Sorenson, Parameter Estimation. Principles and Problems, Marcel Dekker, New York, 1980. [14] M. Basseville, I.V. Nikiforov, Detection of Abrupt Changes—Theory and Application, Prentice-Hall, Englewood Cliffs, NJ, 1993. [15] S.M. Kay, Fundamentals of Statistical Signal Processing, Detection theory, Prentice-Hall, Upper Saddle River, NJ, 1988. [16] S. Sharma, Applied Multivariate Techniques, John Wiley & Sons, New York, 1996. [17] D.C. Montgomery, Introduction to Statistical Quality Control, third ed., John Wiley & Sons, New York, 1997. [18] J. Kullaa, Elimination of environmental influences from damage-sensitive features in a structural health monitoring system, in: Proceedings of the First European Workshop on Structural Health Monitoring, Paris, France, 2002, pp. 742–749.
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp
Sensor validation using minimum mean square error estimation Jyrki Kullaa a,b, a b
Department of Applied Mechanics, Helsinki University of Technology, Finland Department of Mechanical and Production Engineering, Helsinki Metropolia University of Applied Sciences, P.O. Box 4021, FIN-00079 Metropolia, Finland
a r t i c l e i n f o
abstract
Article history: Received 14 February 2009 Received in revised form 26 November 2009 Accepted 8 December 2009 Available online 21 December 2009
Sensor fault can be detected and corrected in a multichannel measurement system with enough redundancy using solely the measurement data. A single or multiple sensors can be estimated from the remaining sensors if training data from the functioning sensor network are available. The method is based on the minimum mean square error (MMSE) estimation, which is applied to the time history data, e.g. accelerations. The faulty sensor can be identified and replaced with the estimated sensor. Both spatial and temporal correlation of the sensors can be utilized. Using the temporal correlation is justified if the number of active structural modes is larger than the number of sensors. The disadvantages of the temporal model are discussed. Experimental multichannel vibration measurements are used to verify the proposed method. Different, and also simultaneous, sensor faults are studied. The effects of environmental variability and structural damage are discussed. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Sensor fault Minimum mean square error estimation Structural health monitoring
1. Introduction Structures equipped with sensors are becoming common due to the development of the sensor and sensor network technology as well as different applications exploiting sensor information. For example, vibration-based structural health monitoring is heavily based on the measurement data recorded during a long period. It is therefore most important that the sensors are working properly and possible sensor faults are detected in order to maintain the reliability of the system. Otherwise a sensor fault can be misinterpreted as a structural damage. Monitoring systems typically include several sensors at different locations of the structure in order to extract features for damage detection or damage localization. Such sensor network constitutes a redundant system, which can be used to detect sensor malfunction or failure, identify, and even correct the faulty sensor. These topics are discussed in this study. Dunia et al. [1] studied detection, isolation, and reconstruction of a faulty sensor using principal component analysis (PCA). They reconstructed each sensor in turn using the remaining sensors, and calculated a sensor validity index from residuals before and after reconstruction for fault detection and isolation. They proposed two approaches to reconstruct a single sensor, an iterative method and optimization, both resulting in the same closed-form solution. Different types of fault, bias, complete failure, drifting, and precision degradation were studied. Different residuals were investigated. They proposed a sensor validity index (SVI) to identify the faulty sensor by replacing one sensor in turn with the estimated sensor. The minimum SVI was found in the case where the faulty sensor was replaced with the estimated one. Dunia and
Corresponding author at: Department of Mechanical and Production Engineering, Helsinki Metropolia University of Applied Sciences, P.O. Box 4021, FIN-00079 Metropolia, Finland, Tel.: + 358 207836176; fax: +358 207836102. E-mail address: jyrki.kullaa@metropolia.fi
0888-3270/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2009.12.001
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1445
Qin [2] also derived conditions for fault detectability, reconstructability, and identifiability. Both sensor and process faults could be identified, but the fault direction must be known, which is easy for a sensor fault, but difficult for process faults. Kerschen et al. [3,4] also applied PCA to sensor validation and used the angle between the principal subspaces as a feature for fault detection. The principal component analysis was made separately for each time record. They also derived a closed-form equation to estimate a single sensor. The faulty sensor was identified by removing one sensor in turn. The faulty sensor was the removed sensor in the case with the minimum angle. Friswell and Inman [5,6] introduced a sensor validation method based on the modal model of the structure. The advantage is that the sensors are not assumed to be functioning initially but instead a modal model has to be available. Two approaches were studied, modal filtering and PCA. They concluded that the modal filtering approach performs better if an accurate modal model is available. They also found that a multiplicative sensor fault is more difficult to locate than an additive fault. Abdelghani and Friswell [7,8] studied model-based methods to validate sensors with additive or multiplicative faults. Two residual generation techniques were studied with an additive fault, the parity space and the modal filtering approaches. For the multiplicative fault, they introduced a correlation index to isolate the faulty sensor. The index was computed from proposed new residuals. Hernandez-Garcia and Masri [9] compared PCA, independent component analysis (ICA), and modified ICA (MICA) for sensor fault detection. They used two features for fault detection: the Hotelling’s T2 statistic and the squared prediction error (SPE). They concluded using a numerical model that the detection performance of ICA and MICA was higher than that of PCA. A disadvantage of PCA is that the number of principal components must be determined, which is not always straightforward with noisy data. It has also been reported that the estimation of several variables simultaneously is complex. This is needed if more than one sensor is faulty or if temporal correlation is also applied. Kullaa [10] introduced factor analysis for sensor validation and compared it with PCA. Similar to PCA, factor analysis also has one parameter, the number of factors, which has to be chosen. However, the sensitivity to the parameter was seen to be lower than in PCA. Also, a closed-form solution to reconstruct several sensors using the remaining sensors was derived. Kullaa [11] introduced the minimum mean square error (MMSE) estimation to sensor validation. It works directly in the data space without a need to determine any model parameters, but with a cost of a slower performance. Also, a spatiotemporal extension was introduced [12], which can be applied if the number of sensors is lower than the number of active modes in the structure. The estimation error was also derived but not applied in the study. The faulty sensor was identified by removing one sensor in turn and performing MMSE estimation to the remaining sensors. The faulty sensor was the missing sensor in the analysis with the lowest mean-square residual compared to that of the training data. The main disadvantage is that the analysis has to be made for each sensor resulting in a slow performance. In this paper, a considerably faster algorithm for sensor validation is derived based on the MMSE estimation. The performance of both the sensor estimation and faulty sensor identification are increased. Also, the estimation error is utilized in the faulty sensor identification, which results in a higher reliability compared to the previously proposed method [11]. In most of the aforementioned studies it was assumed that the number of sensors is higher than the number of excited modes. If the number of active modes is larger than the number of sensors, there is a lack of redundancy between sensors. The time histories could then be filtered to a limited bandwidth. Another solution is to use temporal correlation together with spatial correlation. Spatial correlation contains information between the sensors at the same time instant (a snapshot of the structural motion), whereas temporal correlation contains frequency information. Bearing this in mind, it would be advantageous to use spatial correlation only to distinguish between structural damage and sensor fault. Natural frequencies have been observed to be more sensitive to structural damage or environmental variations than the mode shapes. These effects are also studied in this paper. The linear MMSE estimator model is derived in the next section followed by the spatiotemporal correlation model. Algorithms for sensor fault detection, identification, and correction are proposed. The method uses solely the measurement data, no numerical model of the structure is needed. The model is then applied to an experimental study of a wooden bridge. It is first assumed that one sensor only is faulty but later simultaneous sensor faults are also investigated. Different fault types, bias, drifting, precision degradation, and gain are studied. Additive and multiplicative faults result in a change of different statistics. Also, the effects of the environment and minor damage are studied. The spatiotemporal correlation model is compared to the purely spatial correlation model. Finally, conclusions and practical suggestions from the research are given. 2. Minimum mean square error (MMSE) estimator Sensor validation is approached by estimating one or more sensors using the remaining ones. The likelihood ratio test is then applied to detect sensor fault. It is assumed that training data have been acquired from the sensor network, in which all sensors are functioning. It is also assumed that the number of active structural modes is less than the number of sensors. This is a necessary condition for enough redundancy. Later, this condition can be relaxed if temporal correlation is also taken into account. In multivariable measurements with enough redundancy, a subset of observation can be estimated using the remaining variables. Each observation is divided into observed variables v and missing variables u: u x¼ ð1Þ v
ARTICLE IN PRESS 1446
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
with a partitioned covariance matrix (estimated using the training data): " #
Ruu Rvu
R¼
Ruv Rvv
ð2Þ
The objective is to estimate the missing variables u using the observed variables v. There are two cases in which the estimator of u is linear [13]: 1. the data x are multinormally distributed, 2. the estimator u^ is constrained to be a linear function of the observed data: u^ ¼ Kv
ð3Þ
where K is an unknown matrix, and u and v are assumed to be random vectors with zero means and known covariances, but the joint distribution need not be Gaussian. Both cases result in identical estimators [13]. Here we consider case 2 only, as the assumption of a Gaussian distribution is not generally valid, e.g. for harmonic vibration. The optimal linear estimator u^ is obtained by minimizing the mean-square error (MSE) eMS ^ T ðuuÞÞ ^ eMS ¼ EððuuÞ
ð4Þ
where E( ) is an expectation operator. The mean-square error is minimized by choosing [13]: K ¼ Ruv R1 vv
ð5Þ
and finally EðujvÞ u^ ¼ Kv ¼ Ruv R1 vv v
ð6Þ
where u|v denotes u given v. The error covariance matrix is [13]: ^ ^ T Þ ¼ Ruu Ruv R1 covðujvÞ EððuuÞðu uÞ vv Rvu
ð7Þ
The main disadvantage of using Eq. (6) for sensor fault detection and identification is that when each variable is estimated in turn, a large matrix Rvv has to be repeatedly inverted. The dimension of the matrix is (p 1) (p 1), where p is the number of sensors. This makes the algorithm slow with a large sensor network. In the following, a faster algorithm is derived, which involves inverting a matrix with a low dimension (usually one!) for each sensor. The inverse of the full covariance matrix is nevertheless needed but it is computed only once. 2.1. Linear estimator Modified formulas for the linear estimator (Eqs. (6) and (7)) can be obtained using the precision matrix C, which is the inverse of the covariance matrix R (C = R 1). Using the matrix inversion in block form, the partitioned covariance matrix in terms of the partitioned precision matrix is " # " #1
Ruu Rvu "
¼
Ruv Cuu ¼ Rvv Cvu
Cuv Cvv
1 ðCuu Cuv C1 vv Cvu Þ
1 ðCuu Cuv C1 Cuv C1 vv Cvu Þ vv
1 1 C1 vv Cvu ðCuu Cuv Cvv Cvu Þ
1 1 1 C1 Cuv C1 vv þ Cvv Cvu ðCuu Cuv Cvv Cvu Þ vv
#
ð8Þ
Applying the matrix inversion lemma, and after some manipulation, the necessary submatrices can be shown to be 1 1 1 Ruu ¼ C1 Cvu C1 uu þ Cuu Cuv ðCvv Cvu Cuu Cuv Þ uu
ð9Þ
1 R1 vv ¼ Cvv Cvu Cuu Cuv
ð10Þ
1 1 Ruv ¼ C1 uu Cuv ðCvv Cvu Cuu Cuv Þ
ð11Þ
1 Rvu ¼ RTuv ¼ ðCvv Cvu C1 Cvu C1 uu Cuv Þ uu
ð12Þ
Substituting Eqs. (9)–(12) into Eqs. (6) and (7) results in 1 EðujvÞ ¼ Ruv R1 vv v ¼ Cuu Cuv v
ð13Þ
1 covðujvÞ ¼ Ruu Ruv R1 vv Rvu ¼ Cuu
ð14Þ
and
The dimension of Cuu is typically much lower than that of Rvv resulting in a faster algorithm. The precision matrix C is first computed by inverting the full covariance matrix R, but it is made only once.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1447
With enough redundancy in the sensor network, ujv can be assumed normally distributed pðujvÞ ¼ NðKv; WÞ
ð15Þ
where K ¼ C1 uu Cuv
W ¼ C1 uu
ð16Þ
The application of the MMSE estimation to sensor fault detection, identification, and reconstruction is discussed in the following sections. 3. Spatiotemporal correlation The variables in the sensor fault identification are the sequential sensor values from each sensor. If the process can be assumed stationary with zero mean, temporal correlation can also be utilized. The time-shifted covariance matrix estimate, with a time shift i, is computed by Ri ¼
N i 1 X xðkÞT xðkþ iÞ Ni k ¼ 1
ð17Þ
where x(k)= x(kDt) is the kth sample at a time instant kDt, where Dt is the sampling period and N is the number of samples. If the sensor variables are iT h ð18Þ x ¼ xðkmÞT ; . . . ; xðk1ÞT ; xðkÞT ; xðkþ 1ÞT ; . . . ; xðkþ nÞT then the covariance matrix 2 R0 R1 6 R0 6 R1 R¼6 6 ^ ^ 4 Rmn Rmn þ 1
estimate R is 3 2 R0 Rm þ n 7 6 T Rm þ n1 7 6 R1 7¼6 7 6 ^ & ^ 5 4 RTm þ n R0
R1
R0
^
&
RTm þ n1
Rm þ n
3
Rm þ n1 7 7 7 7 ^ 5 R0
ð19Þ
where m+n is the model order. Hence, instead of the time shifts m and n, one parameter only, the model order, has to be defined. The last equality in Eq. (19) results from the fact that Ri ¼ RTi If the sensor to be investigated is sensor j, then the missing variables are h iT u ¼ xj ðkmÞT ; . . . ; xj ðk1ÞT ; xj ðkÞT ; xj ðk þ 1ÞT ; . . . ; xj ðkþ nÞT
ð20Þ
ð21Þ
corresponding to row j and column j of each Ri. All these missing values are estimated simultaneously. As a result, the number of each estimated variable is equal to m+n + 1 except for the beginning and end of the time histories. A more accurate estimate is obtained by averaging. It should be noted that appending two time histories from different measurements and using a model order larger than zero, an erroneous estimation will occur around the time discontinuity. There are several disadvantages in using the spatiotemporal model instead of the spatial model: (1) The process has to be stationary; (2) the model order has to be determined; (3) the time discontinuities introduce estimation error; (4) an increased number of variables results in a slower performance; and (5) separation of structural fault from a sensor fault may be difficult. Therefore, spatiotemporal correlation should only be used if the number of sensors is inadequate for required redundancy. In that case, filtering can be used as an alternative to decrease the active modes in the data. 4. Sensor fault detection, identification, and correction This section discusses the application of the presented theory to sensor validation: fault detection, faulty sensor identification, and sensor reconstruction. 4.1. Fault detection The first step in the sensor validation is fault detection. The MMSE model is built using the training data from the functioning sensor network. This model is the null hypothesis H0, and the hypothesis test is H0 :
EðujvÞ ¼ K0 v; covðujvÞ ¼ W0
H1 :
EðujvÞaK0 v; covðujvÞaW0
ð22Þ
ARTICLE IN PRESS 1448
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
where K0 and W0 are estimated from the training data. The generalized likelihood ratio test (GLRT) [14,15] decides H1, if LG ðujvÞ ¼
X
ln
pðujv; H1 Þ 4g pðujv; H0 Þ
ð23Þ
where g is some threshold and assuming u|v is normally distributed, pðujv; H1 Þ ¼ ln pðujv; H1 Þln pðujv; H0 Þ pðujv; H0 Þ 1 jW1 j 1 1 ðuK1 vÞT W1 ðuK0 vÞT W1 ¼ ln 1 ðuK1 vÞ þ 0 ðuK0 vÞ 2 jW0 j 2 2
ln
ð24Þ
where K1 and W1 are estimated separately for each measurement. The feature used for fault detection is the log-ratio given by Eq. (24). The number of variables is equal to the number of sensors. This number is usually too large for a statistical analysis and a dimensionality reduction is made using the principal component analysis (PCA) [16]. The largest principal components point at the directions of the largest variance in the multidimensional space and they are finally used to design a control chart [17] for sensor fault detection. For different types of sensor faults, a suitable fault indicator depends on whether the fault affects the mean or the variance. Shewhart chart: The Shewhart chart, or xbar chart [17], monitors the change of the mean value, and the plotted variable is the subgroup mean x. The upper and lower control limits, UCL and LCL, respectively, are: ) UCL ð25Þ ¼ x 7 A3 S LCL where x ¼ meanðxÞ is the average of the subgroup means, S is the average of the subgroup standard deviations, and A3 is a constant depending on the subgroup size. Here, the subgroup size is 20, and A3 is 0.6797. S Chart: The S chart monitors the variability [17] and the plotted variable is the sample standard deviation S. The upper and lower control limits, UCL and LCL, respectively, are UCL ¼ B4 S LCL ¼ B3 S
ð26Þ
where B3 and B4 depend on the subgroup size. With a subgroup size of 20, B3 =0.5102 and B4 = 1.4898. Shewhart T chart: The multivariate counterpart of the Shewhart chart is the Shewhart T control chart [17], where the plotted characteristic is: T 2 ¼ nðxx ÞT S1 ðxx Þ
ð27Þ
where x is the subgroup average, x is the process average, which is the mean of the subgroup averages when the process is in control, and S is the matrix consisting of the grand average of the subgroup variances and covariances. The upper control limit is UCL ¼
pðmþ 1Þðn1Þ Fa;p;mnmp þ 1 mnmp þ 1
ð28Þ
where p is the dimension of the variable, n is the subgroup size, m is the number of subgroups, when the process is assumed to be in control, and Fa,p,mn m p + 1 denotes the a percentage point of the F distribution with p and mnmp+1 degrees-offreedom. The MMSE estimation is applied to sensor fault detection using the following procedure. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.
Define the training data with no sensor faults for hypothesis H0. Estimate mean l and covariance matrix R of the data. Compute the precision matrix: C = R 1. Select sensor and set it as a missing variable u. Form matrices lu, lv, Cuu, Cvv, Cuv, and Cvu by partitioning l and C. Compute K and C using Eq. (16). Return to 4 until all sensors have been modelled. Define the test data (typically one time series) for hypothesis H1. Repeat steps 2–7 for the test data. Compute the log-ratio (24) for each sample. Define the in-control samples with no sensor faults. Perform a dimensionality reduction using the principal component analysis. Perform a statistical analysis using control charts. If the control chart signals, it is an indication of sensor fault. Move on to fault isolation and reconstruction.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1449
4.2. Faulty sensor identification Once a sensor fault is detected, the faulty sensor must be identified. It is proposed that the largest log-likelihood ratio reveals the faulty sensor. The MMSE estimation is applied to faulty sensor identification using the following procedure. 1. Compute the log-likelihood ratio for each sensor using Eq. (23). 2. Find the largest log-likelihood ratio, which indicates the faulty sensor. 4.3. Sensor correction The sensor correction can be made once the faulty sensor has been isolated. The model is built using the training data from all sensors without fault. The faulty sensor is reconstructed by removing the sensor from the network and estimating it using the remaining sensors. The procedure is as follows. 1. 2. 3. 4.
Set the faulty sensor as the missing variable u. Compute the mean and variance of u|v for each sample using Eqs. (13) and (14), respectively. EðujvÞ is the reconstructed sensor. covðujvÞ can be used to plot the confidence limits.
5. Experimental results Sensor faults were investigated with a monitoring system built in the laboratory. The structure was a wooden bridge shown in Fig. 1 with a total mass of 36 kg. A random excitation was generated with an electrodynamic shaker to activate the vertical, transverse, and torsional modes. Fifteen Kistler 8712A5M1 piezo-electric accelerometers were used for the response measurements. They were located at three different longitudinal positions, each containing five accelerometers: vertical and transverse accelerations were measured at both top flanges but transverse acceleration at one of the bottom flanges only. The sampling frequency was 256 Hz and the measurement period was 32 s. The modal parameters were identified from the output-only data, which resulted in 20 natural modes below 128 Hz. As this number was larger than the number of sensors, there was lack of redundancy. Therefore, the signals were low-pass filtered below 64 Hz and re-sampled. In this frequency range, the structure contained 14 natural modes. The data used for sensor validation were three subsequent measurements and a fourth measurement acquired later on the same day. These four records were appended after each other. Different fault types in different sensors were studied: bias, drifting, precision degradation, and gain. The sensor fault was introduced to the fourth measurement. The first two measurements were used to build the model and to design the control charts.
Fig. 1. Wooden bridge.
ARTICLE IN PRESS 1450
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
In sensor fault detection, the MMSE estimation resulted in log-ratios (24) for each sensor. The dimensionality of the feature vector was thus 15, and further dimensionality reduction was made using PCA. The largest principal component was then used in sensor fault detection using control charts. Once a fault was detected, further procedures were performed to identify and correct the faulty sensor. Also, the effects of the environmental variability and structural damage (point masses) on the sensor validation were studied. 5.1. Gain fault A gain fault was studied multiplying sensor 14 by a constant 1.3. Because the mean acceleration is zero, the gain fault affects the change in the variance and the S chart is a proper fault indicator. The S control chart for the largest principal component of the log-ratios (24) is shown in Fig. 2a to detect a sensor fault. Notice the logarithmic scale. The subgroup size was 20. The sensor fault was clearly detected. The faulty sensor identification based on the maximum likelihood ratio is shown in Fig. 3, plotting the log-likelihood ratio for each sensor. The faulty sensor was correctly identified by the proposed method. After identifying the faulty sensor, it was removed and the remaining sensors were re-analyzed. The resulting control chart in Fig. 2b is considerably different indicating no clear fault. It can be concluded that the fault was localized to sensor 14. 5.2. Spatiotemporal analysis Spatiotemporal correlation model was studied with sensor data in which there was not enough redundancy. The sampling frequency was 256 Hz and no low-pass filtering was applied. The number of active modes was thus higher than the number of sensors.
S
101
100 100
200
300
400
500
600
700
800
600
700
800
Sample Number
S
101
100 100
200
300
400
500
Sample Number Fig. 2. S control charts for sensor fault detection. Gain fault in sensor 14: (a) all sensors included and (b) sensor 14 removed.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1451
Faulty Sensor Identified: 14 1000
Log−Likelihood Ratio
800
600
400
200
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number Fig. 3. Faulty sensor identification. Gain fault in sensor 14.
102 102
Tave
Tave
100
100
10−2
10−4 10−2 500
1000 Sample Number
500
1500
1000 Sample Number
1500
Tave
102
100
10−2
500
1000 Sample Number
1500
Fig. 4. Shewhart T control charts for sensor fault detection. Bias fault in sensor 1: (a) model order 0 (spatial correlation); (b) model order 5 and (c) model order 5 with sensor 1 removed.
ARTICLE IN PRESS 1452
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
Faulty Sensor Identified: 1
Log−Likelihood Ratio
250
200
150
100
50
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number x 105
Faulty Sensor Identified: 1
8
Log−Likelihood Ratio
7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number Fig. 5. Faulty sensor identification. Bias fault in sensor 1: (a) model order 0 (spatial correlation) and (b) model order 5.
Bias fault was simulated by adding a constant value 0.01 to sensor 1. This resulted in a change in the mean and the Shewhart T control chart for the largest principal component of the log-ratios was applied to detect the sensor failure. With spatial correlation (model order= 0), the sensor fault was detected (Fig. 4a), but using also the temporal correlation with a model order 5, the fault was more clearly detected (Fig. 4b, notice the logarithmic scale). In either case, the faulty sensor was correctly identified from the maximum log-likelihood (Fig. 5). Also, the reconstruction was more accurate using the spatiotemporal model. Fig. 6 shows the true and estimated sensor (top, a detail), their difference (bottom), and the difference between the measured (faulty) and estimated sensor (middle) from both the spatial and spatiotemporal analyses. The sensor estimated using the spatiotemporal model agrees well with the true sensor, whereas the estimation error is relatively high using the spatial model. After identifying the faulty sensor, it was removed and the remaining data were re-analyzed. The resulting control chart is shown in Fig. 4c, which now indicates no fault. It can be thus concluded that the fault was successfully localized to sensor 1.
5.3. Effects of environment and damage A relevant question in structural health monitoring is if sensor fault can be distinguished from structural damage or environmental or operational variability. In the cases above, different sensor faults were detected from the data recorded at similar environmental and structural conditions as the training data. In this section, the fourth measurement is taken from a different environmental condition than the training data. Also, minor structural change was present with an added mass of size 0.197 kg on the structure with a total mass of 36 kg. The
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1453
Reconstructed sensor: 1 0.1 0 −0.1 0.1 0 −0.1 1
2
3
4
1
2
3 Measurement
4
0.1 0 −0.1
Reconstructed sensor: 1 0.1 0 −0.1 0.1 0 −0.1 1
2
3
4
1
2
3 Measurement
4
0.1 0 −0.1
Fig. 6. Sensor validation with bias fault of sensor 1. Model order is (a) 0 (spatial correlation) and (b) 5. Detail of the true (blue) and estimated (green) sensor (top); residual: measured – estimated (middle); residual: true – estimated (bottom). The red lines are the estimated 3s error bounds. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
effects of these changes on the modal properties have been found to be significant [18]. In the present study, the objective is to separate a sensor fault from other changes in the system. In this study, sensor 10 suffered from drifting (a linear change in the mean from 0 to 0.01 within measurement). The fault was detected and the faulty sensor correctly identified using a spatial correlation model. The control chart is shown in Fig. 7a showing a clear change in the system. Sensor 10 was then removed resulting in the control chart shown in Fig. 7b and suggesting sensor 9 is faulty. Also, removing sensor 9 resulted in the control chart shown in Fig. 7c. Now, looking at all three charts it can be seen that removing sensor 10 had a large influence causing the drifting effect to vanish, whereas a further removal of sensor 9 did not make any enhancement. Therefore, it can be concluded that sensor 10 was faulty but there were also other influences present.
5.4. Simultaneous sensor faults Multiple sensor faults were studied introducing different faults to three sensors simultaneously: precision degradation (s = 0.01) to sensor 2, drifting (linear change of mean from 0 to 0.1) to sensor 5, and gain (multiplication by 1.3) to sensor 14. The data were filtered below 64 Hz. The sensor fault was easily detected, but the main interest here was the identification of the faulty sensors. Increasing the estimated number of faulty sensors, the possible combinations to be investigated increased. For example, with 15 sensors, estimating 1, 2, 3, or 4 faulty sensors, the number of combinations was 15, 105, 455, and 1365, respectively. The results of faulty sensor identification with different number of estimated faulty sensors are shown in Fig. 8, each bar representing a different sensor combination. If one sensor only was assumed faulty, sensor 5 was found (Fig. 8a). With two faulty sensors estimated, sensors 2 and 5 were found (Fig. 8b). With the correct number of three faulty sensors, all faulty sensors were correctly identified as faulty (Fig. 8c). If four faulty sensors were searched for, sensors 2, 5, 13, and 14 were
ARTICLE IN PRESS 1454
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
0
0
−20
−20
−40
−40 xave
xave
−60 −80 −100
−60 −80
−120
−100
−140 −120
−160 200
400 600 Sample Number
200
800
400 Sample Number
600
800
0
xave
−20
−40
−60
−80 100
200
300 400 500 Sample Number
600
700
800
Fig. 7. Shewhart control charts for sensor validation with measurement 4 at different environment and with a minor structural damage. Drifting of sensor 10: (a) all sensors present; (b) sensor 10 removed and (c) sensors 9 and 10 removed.
identified as faulty (not shown). Using spatiotemporal correlation with model order 4 also resulted in a correct identification with a higher sensitivity, but at a cost of a higher computational time (Fig. 8d). The reconstructed and true sensors 2, 5, and 14 are shown in Fig. 9 using spatial correlation. It can be seen that sensors 2 and 5 were accurately reconstructed, whereas the estimation error of sensor 14 was higher. This was also anticipated as can be seen from the estimated 3s error bounds in Fig. 9b showing a detail of the sensor estimation. It is interesting that if one sensor only was assumed faulty, the residuals of the other two faulty sensors 2 or 14 showed no remarkably high values (Fig. 8a) and thus could not be identified. Therefore, the correct number of faulty sensors has to be known for a reliable identification of all faulty sensors. Estimating the correct number of faulty sensors remains a subject for a further study.
6. Conclusion Sensor fault detection, identification, and correction were studied using the minimum mean square error (MMSE) estimation with spatial and spatiotemporal correlation between the variables. The proposed method was validated using experimental vibration data. The following observations were made. Different statistics were used to detect additive and multiplicative sensor faults. Additive faults (drifting and bias) can be detected from the change in the mean of the log-ratios, whereas multiplicative faults (gain and precision degradation) are detected from the change in the variance of the log-ratios. Therefore, two types of control charts were used: the univariate or multivariate Shewhart chart to monitor the mean, and the univariate S chart to monitor the variance. In practice, it is not known in advance what type of sensor fault is present, and therefore both charts should be constructed.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
x 107
Faulty Sensor Identified: 5
3
3
2.5
2.5
Log−Likelihood Ratio
Log−Likelihood Ratio
x 107
2 1.5 1 0.5
1455
Faulty Sensors Identified: 2 5
2 1.5 1 0.5
0
0 0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sensor number x 10
7
20
40
60
80
100
Sensor Combination
Faulty Sensors Identified: 2 5 14
x
107
Faulty Sensors Identified: 2 5 14
3 8 Log−Likelihood Ratio
Log−Likelihood Ratio
2.5 2 1.5 1
6
4
2 0.5 0
0 0
100
200
300
Sensor Combination
400
0
100
200 300 Sensor Combination
400
Fig. 8. Faulty sensor identification with different numbers of estimated faulty sensors: (a) 1, (b) 2 (c) 3 using spatial correlation and (d) 3 using spatiotemporal correlation with model order 4. Each bar represents one sensor combination. Sensors 2, 5, and 14 are faulty.
A faster algorithm for the sensor estimation was developed compared to the previously introduced method. The results are identical, but the proposed method requires that the full covariance matrix be inverted only once and the computational performance for each sensor is high. In the previous method, each sensor estimation required an inverse of a large matrix. The estimation resulted in the expected value of each sensor as well as the error covariance. The expected value can be used as the reconstruction of a faulty sensor. The error covariance can be used to assess the confidence limits. Because the origin of the sensor data in this study was structural vibration, the sensor data also had temporal correlation. This additional correlation can be utilized if the number of active modes is larger than the number of sensors, for which the spatial correlation model may result in poor estimation. On the other hand, if there is enough redundancy, the spatial correlation model should be preferred, because the spatiotemporal model includes several disadvantages or assumptions that may be violated. The proposed method was relatively easy to expand to simultaneous sensor faults. The method could find all faulty sensors if the number of faulty sensors was known. With an estimated lower number of faulty sensors, all sensors that were identified as faulty were indeed faulty. The computational effort increases considerably with an increasing number of faulty sensors. Therefore, there is a practical upper limit of faulty sensors that can be found in a large sensor network. Also, if several sensors are assumed faulty, the remaining number of sensors used for estimation decreases, resulting in a higher estimation error.
ARTICLE IN PRESS 1456
J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
Reconstructed sensor: 2 5 14 0.05 0 −0.05 1
2
3
4
1
2
3
4
1
2
3 Measurement
4
0.05 0 −0.05 0.1 0 −0.1
Reconstructed sensor: 2 5 14 0.05 0 −0.05 0.05 0 −0.05 0.1 0 −0.1 Time Fig. 9. Reconstructed (red) and true (blue) sensors 2 (top), 5 (middle), and 14 (bottom): (a) all measurements; (b) detail with estimated 3s error bounds (black dotted line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
A method was proposed to estimate if the detected fault was due to a sensor fault or another effect. The identified faulty sensor was removed from the set and the remaining sensors were re-analyzed. If the control chart showed a significant change, the removed sensor was probably faulty. On the other hand, if the control chart did not change much, the primary source of change was elsewhere (global effect). It should be noted, however, that the effects of environment or minor damage could not be removed. This problem remains a subject for a future study. It is anticipated that in the future more structures are equipped with a sensor network having numerous sensors. This fact makes the proposed method valuable and available for practical applications.
Acknowledgements This research was performed in the Multidisciplinary Institute of Digitalisation and Energy (MIDE) research programme launched at Helsinki University of Technology to celebrate its centennial as a university. References [1] R. Dunia, S.J. Qin, T.F. Edgar, T.J. McAvoy, Identification of faulty sensors using principal component analysis, AiChe Journal 42 (10) (1996) 2797–2812. [2] R. Dunia, S.J. Qin, Joint diagnosis of process and sensor faults using principal component analysis, Control Engineering Practice 6 (1998) 457–469. [3] G. Kerschen, P. De Boe, J.-C. Golinval, K. Worden, Sensor validation for on-line vibration monitoring, in: Proceedings of the Second European Workshop on Structural Health Monitoring, Munich, Germany, 2004, pp. 819–827. [4] G. Kerschen, P. De Boe, J.-C. Golinval, K. Worden, Sensor validation using principal component analysis, Smart Materials and Structures 14 (2005) 36–42.
ARTICLE IN PRESS J. Kullaa / Mechanical Systems and Signal Processing 24 (2010) 1444–1457
1457
[5] M.I. Friswell, D.J. Inman, Sensor validation of smart structures, in: Proceedings of IMAC-XVIII: A Conference on Structural Dynamics, San Antonio, Texas, 2000, pp. 483–489. [6] M.I. Friswell, D.J. Inman, Sensor validation of smart structures, Journal of Intelligent Material Systems and Structures 10 (1999) 973–982. [7] M. Abdelghani, M.I. Friswell, Sensor validation for structural systems with additive sensor faults, Structural Health Monitoring 3 (3) (2004) 265–275. [8] M. Abdelghani, M.I. Friswell, Sensor validation for structural systems with multiplicative sensor faults, Mechanical Systems and Signal Processing 21 (2007) 270–279. [9] M.R. Hernandez-Garcia, S.F. Masri, Multivariate statistical analysis for detection and identification of faulty sensors using latent variable methods, Advances in Science and Technology 56 (2008) 501–507. [10] J. Kullaa, Three models for sensor validation, in: Proceedings of the 7th International Workshop on Structural Health Monitoring, Stanford, CA, 2009, pp. 529–536. [11] J. Kullaa, Sensor fault identification and correction in structural health monitoring, in: Proceedings of ISMA2006, International Conference on Noise and Vibration Engineering, Leuven, Belgium, 2006, pp. 873–884. [12] J. Kullaa, Sensor fault identification and correction in vibration-based multichannel structural health monitoring, in: Proceedings of the 6th International Workshop on Structural Health Monitoring, Stanford, CA, 2007, pp. 606–613. [13] H.W. Sorenson, Parameter Estimation. Principles and Problems, Marcel Dekker, New York, 1980. [14] M. Basseville, I.V. Nikiforov, Detection of Abrupt Changes—Theory and Application, Prentice-Hall, Englewood Cliffs, NJ, 1993. [15] S.M. Kay, Fundamentals of Statistical Signal Processing, Detection theory, Prentice-Hall, Upper Saddle River, NJ, 1988. [16] S. Sharma, Applied Multivariate Techniques, John Wiley & Sons, New York, 1996. [17] D.C. Montgomery, Introduction to Statistical Quality Control, third ed., John Wiley & Sons, New York, 1997. [18] J. Kullaa, Elimination of environmental influences from damage-sensitive features in a structural health monitoring system, in: Proceedings of the First European Workshop on Structural Health Monitoring, Paris, France, 2002, pp. 742–749.