Effective fault detection in structural health monitoring systems

Abstract

A new fault detection technique is considered in this article. It is based on kernel partial least squares, exponentially weighted moving average, and generalized likelihood ratio test. The developed approach aims to improve monitoring the structural systems. It consists of computing an optimal statistic that merges the current information and the previous one and gives more weight to the most recent information. To improve the performances of the developed kernel partial least squares model even further, multiscale representation of data will be used to develop a multiscale extension of this method. Multiscale representation is a powerful data analysis way that presents efficient separation of deterministic characteristics from random noise. Thus, multiscale kernel partial least squares method that combines the advantages of the kernel partial least squares method with those of multiscale representation will be developed to enhance the structural modeling performance. The effectiveness of the proposed approach is assessed using two examples: synthetic data and benchmark structure. The simulation study proves the efficiency of the developed technique over the classical detection approaches in terms of false alarm rate, missed detection rate, and detection speed.

Introduction

Structural health monitoring (SHM) is considered as an approach to monitor the structure systems. It consists of developing monitoring solutions based on the data collected by the sensors installed on the structure. Its purpose is to detect, identify, estimate, and/or locate the existing faults on the structure and to make decisions concerning its evolution and its state in order to improve its safety. In the current work, we will focus on the fault detection (FD) problem.

In the literature, many methods for detecting faults^1–3 have been considered. They can be defined as methods based on a model^4–6 or as data-driven methods.^7–12 Many physical systems of the actual civil structure are complex and difficult to construct their models.^13,14 Therefore, data-driven methods are widely used for monitoring and FD purposes. ¹⁵

Latent variable regression (LVR) is a well-known data-driven modeling technique that includes principal component analysis (PCA) ¹⁶ and partial least squares (PLS). ¹⁷ LVR modeling methods are a multivariate analysis that aims to reduce the dimensionality of the data and rely on the definition of a linear data transformation via an orthonormal matrix which is computed from the dataset itself.

However, most of the practical systems are input–output, nonlinear, and multivariate. To make the extension to input–output, nonlinear, and multivariate models, kernel partial least squares (KPLS) will be used for the modeling purposes. The benefits of using the KPLS model lie in the fact that it does not require nonlinear optimization and it only requires the solution of an eigenvalue problem.

To improve the performances of the KPLS model even further, multiscale representation of data will be used to develop a multiscale kernel partial least squares (MSKPLS) method. Multiscale representation is a powerful data analysis way that presents efficient separation of deterministic characteristics from random noise.

Once the model is built, the monitoring technique should be addressed in order to detect the faults in SHM systems. To do that, a novel detection chart based on generalized likelihood ratio test (GLRT) and exponentially weighted moving average (EWMA) chart statistics will be applied. The EWMA-GLRT chart depends upon evaluating the residuals, which are obtained from the MSKPLS to detect the fault in the network.

The benefits of the proposed EWMA-GLRT are to develop a new detection statistic that considers the current information and the previous one and strengthens the more recent information by providing more weight to it. Therefore, the developed FD technique consists of considering the MSKPLS for modeling purposes and EWMA-GLRT chart for detection and monitoring goals. The developed technique will be evaluated and compared to many other techniques using two examples: synthetic data and benchmark structure. The results show the efficiency of the proposed technique over the classical approaches in terms of false alarm rate (FAR), missed detection rate (MDR), and detection speed.

The rest of this article is arranged as follows: first, we describe the multiscale KPLS method in section “Multiscale KPLS description.” In section “OEWMA-GLRT detection chart description,” we present the new EWMA-GLRT chart. In section “Applications,” the effectiveness of the proposed approach is evaluated through synthetic data and a simulated SHM benchmarking data. Finally, we present the conclusions in section “Conclusion.”

Multiscale KPLS description

KPLS

Given an input matrix $X\in {\mathbb{R}}^{N\times m}$ consists of N samples with m process variables and an output matrix $Y\in {\mathbb{R}}^{N\times m}$ consists of N observations with p product quality variables. The PLS technique projects the available process measurements X into a low-dimensional space with $\ell$ latent variables, where $\ell$ can be determined by certain criteria, such as the cumulative percent variance (CPV) criterion. Then, Y can be constructed by these latent variables. The PLS model can be expressed ¹⁸ as follows

$$\begin{array}{l}X=T{P}^T+E=\sum _{i=1}^{\ell }{t}_i{p}_i^T+E\\ Y=U{Q}^T+F=\sum _{j=1}^{\ell }{u}_j{q}_j^T+F\end{array}$$

(1)

where $T=[\begin{array}{cc}{t}_1\hfill & t_2\cdots t_{\ell }\hfill \end{array}]\in {\mathbb{R}}^{N\times \ell }$ and $U=[\begin{array}{cc}{u}_1\hfill & u_2\cdots u_{\ell }\hfill \end{array}]\in {\mathbb{R}}^{N\times \ell }$ are the input and output score matrices, respectively. $P=[\begin{array}{cc}{p}_1\hfill & p_2\cdots p_{\ell }\hfill \end{array}]\in {\mathbb{R}}^{m\times \ell }$ and $Q=[\begin{array}{cc}{q}_1\hfill & q_2\cdots q_{\ell }\hfill \end{array}]\in {\mathbb{R}}^{N\times \ell }$ are the input and the output loading matrices, respectively. $E\in {\mathbb{R}}^{N\times m}$ and $F\in {\mathbb{R}}^{N\times p}$ denote the residual parts of X and Y, respectively, which are obtained from the PLS model.

Note that the matrices X and Y are first standardized to have zero mean and unit variance before constructing the PLS model. PLS assumes linear correlation among variables, which lead to prediction and modeling errors in cases of nonlinear processes. To address this issue and in order to extend this technique to deal with nonlinear input–output models, many extensions have been proposed to define the nonlinearities. The PLS has been extended to nonlinear regression through the use of kernel-based functions. The most well-known approach is called KPLS, which will be considered in the next section.

KPLS model

The KPLS method consists of using the classical PLS directly in the feature space. Its idea first is to map the original nonlinear data into a linear high-dimensional feature space and then to perform linear PLS in the feature space. The main idea of KPLS method is illustrated in Figure 1.

Figure 1.

Basic idea of KPLS.

Define $\varphi$ as a nonlinear map which projects the input vector from the original space into the feature space $\mathcal{F}$ in which they are related linearly approximately. After the nonlinear map, the input matrix X turns to the feature matrix $\mathrm{\Phi }$

$$\mathrm{\Phi }={\left[\begin{array}{cc}\varphi (x_1)\hfill & \varphi (x_2)\cdots \varphi (x_N)\hfill \end{array}\right]}^T\in {\mathbb{R}}^{N\times f}$$

(2)

where the dimension f can be arbitrarily large or even infinite.

In the feature space, $\varphi (x_i)$ $(i=1,2,...,N)$ should be scaled to zero mean

$$\overline{\varphi }(x_i)=\varphi (x_i)-\overline{\varphi }$$

(3)

where

$$\begin{array}{c}\multicolumn{1}{c}{\overline{\varphi }=\frac{1}{N}\sum _{i=1}^N\varphi (x_i)}\\ =\frac{1}{N}\left[\begin{array}{cc}\varphi (x(1))\hfill & \varphi (x(2))\cdots \varphi (x(N))\hfill \end{array}\right]I_N\hfill \\ =\frac{1}{N}{\mathrm{\Phi }}^T{I}_N\hfill \end{array}$$

(4)

with $I_N=[1\cdots 1{]}^T\in {\mathbb{R}}^N$ .

Thus, the zero mean of $\mathrm{\Phi }$ is

$$\begin{array}{c}\multicolumn{1}{c}{\overline{\mathrm{\Phi }}={\left[\begin{array}{cc}\overline{\varphi }(x_1)\hfill & \overline{\varphi }(x_2)\cdots \overline{\varphi }(x_N)\hfill \end{array}\right]}^T}\\ =\mathrm{\Phi }-{\left[\begin{array}{cc}\overline{\varphi }\hfill & \overline{\varphi }\cdots \overline{\varphi }\hfill \end{array}\right]}^T\hfill \\ =\mathrm{\Phi }-\frac{1}{N}{I}_N{I}_N^T\mathrm{\Phi }\hfill \end{array}$$

(5)

According to Sheriff et al., ¹⁹ the KPLS model of $(\overline{\mathrm{\Phi }},Y)$ is modeled as follows

$$\begin{array}{c}\overline{\mathrm{\Phi }}=T{P}^T+{\overline{\mathrm{\Phi }}}_r\hfill \\ \multicolumn{1}{c}{Y=U{Q}^T+Y_r}\end{array}$$

(6)

where $T\in {\mathbb{R}}^{N\times A}$ and $U\in {\mathbb{R}}^{N\times A}$ are the input and output score matrices, respectively. $P\in {\mathbb{R}}^{m\times A}$ and $Q\in {\mathbb{R}}^{N\times A}$ are the loading matrices of $\overline{\mathrm{\Phi }}$ and Y, respectively. ${\overline{\mathrm{\Phi }}}_r$ and $Y_r$ are residuals. A is the number of latent variables in the feature space determined by the CPV criterion. ¹⁹

Note that $\varphi (.)$ is not explicitly defined and the dimension f of $\mathrm{\Phi }$ is arbitrarily large or even infinite. To avoid the explicit using $\mathrm{\Phi }$ , a kernel matrix K is introduced and it is defined as

$$K=\mathrm{\Phi }{\mathrm{\Phi }}^T$$

(7)

where the entry $k_{\mathrm{ij}}$ of the matrix K is given by

$$k_{\mathrm{ij}}={\varphi }^T(x_i)\varphi (x_j)=k(x_i,x_j)$$

(8)

Thus, the kernel matrix K is given by

$$\begin{array}{c}\multicolumn{1}{c}{K=\mathrm{\Phi }{\mathrm{\Phi }}^T}\\ =\left[\begin{array}{ccc}{\varphi }^T(x_1)\varphi (x_1)\hfill & \cdots \hfill & {\varphi }^T(x_1)\varphi (x_N)\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\varphi }^T(x_N)\varphi (x_1)}& \cdots \hfill & {\varphi }^T(x_N)\varphi (x_N)\hfill \end{array}\right]\hfill \\ =\left[\begin{array}{ccc}k(x_1,x_1)\hfill & \cdots \hfill & k(x_1,x_N)\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{k(x_N,x_1)}& \cdots \hfill & k(x_N,x_N)\hfill \end{array}\right]\hfill \end{array}$$

(9)

The function $k(.,.)$ is called the kernel function. Note that there exist several types of these functions.^20,21 The Gaussian kernel function is considered as one of the most used functions which is defined as

$$k(x_i,x_j)=\exp (-\frac{{(x_i{-}_j)}^T(x_i-x_j)}{c})$$

(10)

where c is the width of a Gaussian function.

The kernel matrix K is centralized as

$$\overline{K}=\overline{\mathrm{\Phi }}{\overline{\mathrm{\Phi }}}^T$$

(11)

According to equation (5), $\overline{K}$ can be written as

$$\begin{array}{c}\multicolumn{1}{c}{\overline{K}=[\mathrm{\Phi }-\frac{1}{N}{I}_N{I}_N^T\mathrm{\Phi }]{[\mathrm{\Phi }-\frac{1}{N}{I}_N{I}_N^T\mathrm{\Phi }]}^T}\\ =[E-\frac{1}{N}{I}_N{I}_N^T]\mathrm{\Phi }{\mathrm{\Phi }}^T{[E-\frac{1}{N}{I}_N{I}_N^T]}^T\hfill \\ =[E-\frac{1}{N}{I}_N{I}_N^T]K{[E-\frac{1}{N}{I}_N{I}_N^T]}^T\hfill \end{array}$$

(12)

where I is the identity matrix.

From the centered data matrices K and Y, the regression coefficient matrix B can be obtained as^22,23

$$B={\mathrm{\Phi }}^{T}U(T^T\overline{K}U{)}^{-1}{T}^{T}Y$$

(13)

The prediction of the response variables is given by

$$\hat{Y}=\mathrm{\Phi }B=\mathrm{KU}(T^T\overline{K}U{)}^{-1}{T}^{T}Y$$

(14)

Equation (14) shows that the response variables (outputs) can be obtained from the inner products of the mapped vectors. Hence, for a new observation x of the predictor vector, the outputs are estimated by

$$\begin{array}{c}\multicolumn{1}{c}{\hat{y}=B^T\overline{\varphi }(x)=Y^{T}T{\left[U{(T^T\overline{K}U)}^{-1}\right]}^T\overline{\mathrm{\Phi }}\overline{\varphi }(x)}\\ =Y^{T}T{\left[U{(T^T\overline{K}U)}^{-1}\right]}^T\overline{k}(x)\hfill \end{array}$$

(15)

where $\overline{k}(x)$ is the vector of centered kernel functions. It is defined as

$$\begin{array}{c}\multicolumn{1}{c}{\overline{k}(x)=\overline{\mathrm{\Phi }}\overline{\varphi }(x)}\\ ={\left[\begin{array}{ccccc}\overline{\varphi }(x_1)\hfill & .\hfill & .\hfill & .\hfill & \overline{\varphi }(x_N)\hfill \end{array}\right]}^T\overline{\varphi }(x)\hfill \\ ={\left[\begin{array}{ccccc}{\overline{\varphi }}^T(x_1)\overline{\varphi }(x)\hfill & .\hfill & .\hfill & .\hfill & {\overline{\varphi }}^T(x_N)\overline{\varphi }(x)\hfill \end{array}\right]}^T\hfill \\ ={\left[\begin{array}{ccccc}\overline{k}(x_1,x)\hfill & .\hfill & .\hfill & .\hfill & \overline{k}(x_N,x)\hfill \end{array}\right]}^T\hfill \end{array}$$

(16)

From equations (4), (5), and (16), the following relationships can be defined as

$$\begin{array}{c}\multicolumn{1}{c}{\overline{k}(x)=\overline{\mathrm{\Phi }}\overline{\varphi }(x)}\\ =[\mathrm{\Phi }-\frac{1}{N}{I}_N{I}_N^T\mathrm{\Phi }][\varphi (x)-\frac{1}{N}{\mathrm{\Phi }}^T{I}_N]\hfill \\ =[I_N-\frac{1}{N}{I}_N{I}_N^T]\mathrm{\Phi }[\varphi (x)-\frac{1}{N}{\mathrm{\Phi }}^T{I}_N]\hfill \\ =[I_N-\frac{1}{N}{I}_N{I}_N^T][k(x)-\frac{1}{N}K{I}_N]\hfill \end{array}$$

(17)

where $k(x)$ is the vector of non-centered kernel functions and is defined analogously to $\overline{k}(x)$ . The main steps of KPLS modeling are illustrated in Algorithm 1.

Algorithm 1. KPLS algorithm.

1. Set $i=1$ , $K_1=K$ , and $Y_1=Y$ . 2. Initialize the score-vector ${\mathit{u}}_i(N\times 1)$ of the latent variable $u_i$ of $Y_i$ , as the maximum-variance column of $Y_i$ . 3. Compute the score-vector ${\mathit{t}}_i(N\times 1)$ of the latent variable $t_i$ of ${\mathrm{\Phi }}_i$ , as ${\mathit{t}}_i=K_i{\mathit{u}}_i/||K_i{\mathit{u}}_i||$ , $||{\mathit{t}}_i||=1$ . 4. Regress the columns of $Y_i$ on ${\mathit{t}}_i$ : ${\mathit{c}}_i=Y_i{\mathit{t}}_i$ , where ${\mathbf{c}}_i$ is a weighting vector. 5. Calculate the new score-vector: ${\mathit{u}}_i=Y_i{\mathit{c}}_i/||Y_i{\mathit{c}}_i||$ , $||{\mathit{u}}_i||=1$ . 6. Repeat steps (3) to (5) until the convergence of ${\mathit{t}}_i$ . 7. Deflate the matrices: $K_{i+1}=(I-{\mathit{t}}_i{\mathit{t}}_i^T)K_i(I-{\mathit{t}}_i{\mathit{t}}_i^T)$ , $Y_{i+1}=Y_i-{\mathit{t}}_i{\mathit{t}}_i^T{Y}_i$ . 8. Save the data in the matrices: $T\leftarrow t_i$ , $U\leftarrow u_i$ . 9. Set $i=i+1$ , and return to step (2). Stop when $i>\ell$ , with $\ell$ being the selected number of latent variables.

Multiscale KPLS technique

Teppola and Minkkinen ²⁴ are the first who have developed the linear multiscale PLS, where the multiscale representation was combined with linear PLS model. It has been applied essentially to remove the drift in the data set. Zhang and Hu ²⁵ have studied the nonlinear extension (MKPLS). The wavelet-based multiscale representation of data can ameliorate the KPLS modeling. Similar to the proposed method multiscale principal component analysis (MSPCA) by Bakshi and Top, ²⁶ we shall use in this work the multiscale KPLS algorithm. Similar to PLS, the modeling using KPLS can be proceeded in the mapped feature space.

Given a data set of training data, the variables are decomposed in the feature space through the discrete wavelet transform (DWT). At each individual scale, KPLS model is applied. In order to construct the data, important scale coefficients are selected according to a statistical threshold (see Figure 2). Then, KPLS model is applied on global scale for FD. At each scale, the statistical thresholding is proceeded as a data filtering stage which improves KPLS efficiency. The MSKPLS algorithm is shown in Algorithm 2 and its representation model is shown in Figure 2.

Figure 2.

Representation of MSKPLS model for FD.

Algorithm 2. MSKPLS algorithm.

The input and output data X and Y are each decomposed into coarse approximate scale and detail scales At each scale, KPLS algorithm is applied and the loading vectors and score vectors of X and Y matrices are computed For X and Y and at each scale, only important coefficients (higher then threshold values) are selected X and Y are reconstructed KPLS algorithm is performed on the new X and Y matrices.

Advantages of multiscale data representation

Multiscale representation has the ability to separate the noise from important features in the data. When data are decomposed at multiple scales by passing through low-pass and high-pass filters, noise is effectively separated from the important features. Random noise in a signal is normally present over all the coefficients, whereas deterministic features in the data are captured in a few but relatively large coefficients. The important features in the data are usually captured by the last scaled signals as well as any large wavelet coefficient (in the detail signals), whereas other small wavelet coefficients usually correspond to noise. ²⁶ Thus, multiscale representation provides an effective method for noise-feature separation as shown in Figure 3.

Figure 3.

A schematic diagram of data representation at multiple scales: (a) original signal, (b) first scaled signal, (c) first detail signal, (d) second scaled signal, (e) second detail signal, (f) third scaled signal, and (g) third detail signal. ²⁶

The benefits of multiscale representation present that they can help verify the assumptions of independence, normality, and noise level made by several univariate and multivariate monitoring approaches. ²⁷

The objective of the next section is to present MSKPLS-based optimized exponentially weighted moving average-generalized likelihood ratio test (OEWMA-GLRT) approach for FD purposes. The developed framework is addressed so that the modeling phase is addressed using the MSKPLS and the detection of the faults is achieved using the OEWMA-GLRT chart. The MSKPLS is used to compute the monitored residuals and the OEWMA-GLRT chart (which is presented next) is applied to the evaluated monitored residuals for FD purposes.

OEWMA-GLRT detection chart description

The new detection chart combines the advantages of the EWMA and GLRT charts. It also consists of optimizing the EWMA chart using a multi-objective function. It reduces the MDR, the FAR, and the detection speed $(\mathrm{AR}{L}_1)$ by selecting the optimal EWMA tuning parameters (L and $\lambda$ ).

Classical GLRT detection chart

The GLRT chart is a commonly used hypothesis testing technique for FD in model-based methods.^28,29 Let us denote for an observation vector $Y\in {\mathcal{R}}^m$ follows a Gaussian distribution by $\mathcal{N}(0;{\sigma }^2{I}_m)$ or $\mathcal{N}(\theta \ne 0;{\sigma }^2{I}_m)$ . Note that ${\sigma }^2>0$ represents the known variance and $\theta$ represents the mean vector. The hypothesis test problem can be defined as follows

$$\{\begin{array}{cc}{\mathcal{H}}_0=\{\mathcal{Y}~\mathcal{N}(0;{\sigma }^2{I}_m)\},\hfill & (\mathrm{null}\mathrm{hypothesis})\hfill \\ \multicolumn{1}{c}{{\mathcal{H}}_1=\{\mathcal{Y}~\mathcal{N}(\theta ;{\sigma }^2{I}_m)\},}& (\mathrm{alternative}\mathrm{hypothesis})\hfill \end{array}$$

(18)

The parameter $\theta$ may define the fault value in our study. The likelihood estimate of $\theta$ is calculated through maximizing the GLRT $\mathcal{T}(Y)$ as follows

$$\begin{array}{c}\mathcal{T}(Y)=2\log \frac{\underset{\theta \in {\mathcal{R}}^m}{\sup }{f}_{\theta }(Y)}{f_{\theta =0}(Y)}\hfill \\ =2\log (\underset{\theta }{\sup }\exp (-\frac{||Y-\theta |{|}_2^2}{2{\sigma }^2})/\exp (-\frac{||Y|{|}_2^2}{2{\sigma }^2}))\hfill \\ =\frac{1}{{\sigma }^2}(\underset{\theta }{min}||Y-\theta |{|}_2^2+||Y|{|}_2^2)\hfill \\ =\frac{1}{{\sigma }^2}(||Y-\hat{\theta }|{|}_2^2+||Y|{|}_2^2)=\frac{1}{{\sigma }^2}||Y|{|}_2^2\hfill \end{array}$$

(19)

where $\hat{\theta }$ represents the maximum likelihood estimate of $\theta$ . It is presented as $\hat{\theta }=\arg \underset{\theta }{min}||Y-\theta |{|}_2^2=Y$ , where $||.|{|}_2$ represents the Euclidean norm. The probability density function of Y is defined as

$$f_{\theta }(Y)=\frac{1}{{(2\pi )}^{\frac{m}{2}}{\sigma }^m}\exp \{-\frac{1}{2{\sigma }^2}||Y-\theta |{|}_2^2\}$$

Note that the maximization of the likelihood function corresponds to the maximization of its logarithm. The distribution of the GLRT decision function $\mathcal{T}(Y)$ under ${\mathcal{H}}_0$ permits establishing a statistical test with a preferred FAR, $\alpha$ , where the threshold ${\delta }_{\alpha }^2$ is selected to fulfill the probability of the false alarm as follows

$${\mathcal{P}}_0(\mathcal{T}(Y)\ge {\delta }_{\alpha }^2=\alpha$$

(20)

where ${\mathcal{P}}_0(A)$ represents the probability of an event A when Y is distributed according to the null hypothesis ${\mathcal{H}}_0$ . As Y follows a normal distribution (equation (18)), the statistics T follows the ${\chi }^2$ distribution law which is central under ${\mathcal{H}}_0$ and non-central under ${\mathcal{H}}_1$ . A parameter ${\kappa }_{\theta }$ of non-centrality is defined and is equal to ${\kappa }_{\theta }=1/{\sigma }^2||\theta |{|}_2^2$ . Thus, we need to find the distribution of the GLRT in order to choose a suitable threshold. The test statistic is distributed through a chi-square law^30,31

$$\mathcal{T}(Y)=\frac{1}{{\sigma }^2}||Y|{|}_2^2~{\chi }_m^2$$

(21)

The GLRT technique is derived to detect a shift in the mean (fault), and its explicit asymptotic statistic is computed. It is used to detect an additive fault $\theta$ through computing the maximum probability of detection for a given false alarm.

In this work, the FD problem is considered as a hypothesis test taking into account two possible hypotheses: a null hypothesis ${\mathcal{H}}_0$ of no change, where measurement vector $\mathbf{x}$ has no fault, and an alternative hypothesis ${\mathcal{H}}_1$ of the change point, where $\mathbf{x}$ has a fault. Thus, the statistic has to choose between ${\mathcal{H}}_0$ and ${\mathcal{H}}_1$ for the best detection. The GLRT will be applied to residuals generated from the MSKPLS model.

From equation (18), let $Y=e_j$ for $j=1,\dots ,m$ , then the problem of hypothesis test problem can be presented as

$$\{\begin{array}{cc}{\mathcal{H}}_0=e_j~\mathcal{N}(0,1)\},\hfill & (\mathrm{null}\mathrm{hypothesis})\hfill \\ \multicolumn{1}{c}{{\mathcal{H}}_1=e_j~\mathcal{N}(\theta ,1),}& (\mathrm{alternative}\mathrm{one})\hfill \end{array}$$

(22)

and the GLRT statistic is presented as follows

$$G(e_j)=e_j^2~{\chi }_{1,\alpha }^2,j=1,...,m$$

(23)

In order to choose the threshold $(G_{\alpha })$ , we need to determine its distribution. Since e is normally distributed, the GLRT statistic will follow a Chi-square distribution ${\chi }^2$

$$G~{\chi }^2$$

(24)

If the GLRT statistic is the upper threshold $(G_{\alpha })$ , there is a fault; otherwise, the process is considered to be operating under normal conditions.

OEWMA-GLRT detection chart

The EWMA chart was first developed by Roberts ³² in 1959. It can be computed as follows ³³

$$Z_i=\lambda X_i+(1-\lambda )Z_{i-1},i=1,\dots ,N$$

(25)

where $\lambda$ is a smoothing parameter $(0<\lambda \le 1)$ which represents the memory change of the detection statistic, $X_i$ represents the ith individual observation value, and $Z_0$ represents the initial value which is taken as the in-control mean ${\mu }_0$ .

The EWMA chart has two control limits: an upper control limit $(\mathrm{UCL})$ and a lower control limit $(\mathrm{LCL})$ . When $Z_i$ exceeds these control limits, the EWMA chart detects a fault.

These control limits are computed as follows ³⁴

$$\mathrm{UCL}={\mu }_0+L\sigma \sqrt{\frac{\lambda }{2-\lambda }[1-{(1-\lambda )}^{2i}]}$$

(26)

$$\mathrm{LCL}={\mu }_0-L\sigma \sqrt{\frac{\lambda }{2-\lambda }[1-{(1-\lambda )}^{2i}]}$$

(27)

where L represents the control width of the EWMA statistic Z and $\sigma$ represents the standard deviation of X. For better use of the EWMA method, the smoothing parameter $\lambda$ and control width L need to be carefully chosen. In this article, we propose to develop an optimized exponentially weighted moving average (OEWMA) statistic which computes the EWMA statistic and its control limits with optimized smoothing parameters of EWMA $\hat{\lambda }$ and $\hat{L}$ .

In order to compute the optimal values $\hat{\lambda }$ and $\hat{L}$ , we use a multi-objective optimization (MOO). ³⁵ The MOO is handled using three objective functions: MDR, FAR, and average run length $(\mathrm{AR}{L}_1)$ . ³⁵

The OEWMA-GLRT statistic $(T)$ can be computed as

$$T(Y)=\alpha \frac{G(Y)}{G_{\alpha }}+(1-\alpha )\frac{Z(Y)}{\mathrm{UCL}}$$

(28)

where $G_{\alpha }$ is the control limit of GLRT, UCL is the upper control limit of EWMA statistic, and $\alpha$ is the weight parameter $(0<\alpha <1)$ .

In order to detect the faults in the residual vector R that is obtained from MSKPLS model, MSKPLS-based OEWMA-GLRT is proposed. The filtering and detection are two independent phases implemented in the FD algorithm (see Figure 4). Nevertheless, when these two tasks are simultaneously implemented, this may improve the FD abilities. Thus, we develop MSKPLS-based OEWMA-GLRT where the KPLS model is established based on the wavelet coefficients and the OEWMA-GLRT is applied to detect faults.

Figure 4.

Schematic illustration of the developed MSKPLS-based OEWMA-GLRT approach.

The proposed MSKPLS-based OEWMA-GLRT technique for FD is described in Algorithm 3.

Algorithm 3. MSKPLS-based OEWMA-GLRT algorithm.

1. Split the data into training set and testing set. Training data:  2. Standardize the data to have zero mean and unit variance.  3. Decompose X and Y data into coarse approximate scale and detail scale after selecting the best decomposition depth.  4. Apply the KPLS at each scale.  5. Reconstruct the X and Y data with KPLS.  6. Construct the MSKPLS model.  7. Compute the OEWMA-GLRT statistic T.  8. Compute the OEWMA-GLRT threshold $T_t$ . Testing data:  9. Standardize the data to have zero mean and unit variance. 10. Decompose X and Y data into coarse approximate scale and detail scale after selecting the best decomposition depth. 11. Apply the KPLS at each scale. 12. Reconstruct the X and Y data with KPLS. 13. Construct the MSKPLS model. 14. Compute the OEWMA-GLRT statistic T. 15. If the statistic $T_t$ is under the threshold $T_t$ , the system is under normal conditions. Else, a fault is declared.

Applications

The validation of the developed approach is done through two simulated examples: synthetic data set ¹⁰ and benchmark model International Association for Structural Control-American Society of Civil Engineers (IASC-ASCE).^36–39

FD using simulated example

The simulated example proposed in Bakshi ⁴⁰ contains four variables where the first variable follows a normal distribution with zero mean and unit variance. The simulation is presented in the following equations

$$\{\begin{array}{c}{X}_{t_1}=[-2\times \mathrm{randn}(200,1);\mathrm{randn}(300,1);3\times \mathrm{randn}(500,1)]\hfill \\ \multicolumn{1}{c}{X_{t_2}=\sin {(0.1:0.1:100)}^T+\cos {(0.1:0.1:100)}^T}\\ \multicolumn{1}{c}{X_{t_3}=0.5\times X_{t_1}^2+0.5\times X_{t_2}}\\ \multicolumn{1}{c}{X_{t_4}=0.3\times X_{t_1}+0.7\times X_{t_2}^2}\end{array}$$

(29)

The input data matrix $\tilde{X}$ is contaminated with error that follows a Gaussian distribution with zero mean. This error consists of a measurement noise $v_{k-1}~\mathcal{N}(0,0.2)$ , so that

$$\{\begin{array}{c}{X}_1(t)=X_{t_1}+{\sigma }_v\times \mathrm{randn}(\mathrm{size}(X_{t_1}))\hfill \\ \multicolumn{1}{c}{X_2(t)=X_{t_2}+{\sigma }_v\times \mathrm{randn}(\mathrm{size}(X_{t_2}))}\\ \multicolumn{1}{c}{X_3(t)=X_{t_3}+{\sigma }_v\times \mathrm{randn}(\mathrm{size}(X_{t_3}))}\\ \multicolumn{1}{c}{X_4(t)=X_{t_4}+{\sigma }_v\times \mathrm{randn}(\mathrm{size}(X_{t_4}))}\end{array}$$

(30)

The four model variables generated in the previous equation will be arranged in a matrix of the form

$$X(t)=[X_1(t),X_2(t),X_3(t),X_4(t)]$$

This matrix consists of four columns presenting the variables and $1000$ rows presenting the samples that express the evolution of the model over time. The four input variables’ responses are shown in Figure 5.

Figure 5.

The time evolution of the input variables.

The output variable will be arranged in a data matrix Y and it is generated as the following equation

$$Y(t)=X^2(t)$$

(31)

In order to construct the KPLS model, the input and output data should be scaled first. Then, the data set is divided into training data which consist of $500$ samples and testing data which consist of $500$ samples, which will be used to validate the performance of developed FD technique.

To illustrate the performance of the proposed technique OEWMA-GLRT chart, we compare it with the other techniques such as the GLRT chart and the OEWMA chart. This comparison addresses the linear and the nonlinear cases processed by the PLS model and the KPLS model, respectively. In the KPLS model, the radial basis function (RBF) is used with a parameter $\eta =6$ . The number of the retained kernel principal components has opted for 2 using the CPV criterion. To test the data, a mean shift fault of size $4\sigma$ is added to the first variable $X_1$ of the testing data set at $[400\mathrm{to}500]$ . The detection comparison results between OEWMA-GLRT, OEWMA and GLRT methods for PLS and KPLS models are shown in Figures 6–11 and Table 1.

Figure 6.

Results using PLS-based GLRT method.

Figure 7.

Results using PLS-based OEWMA method.

Figure 8.

Results using PLS-based OEWMA-GLRT method.

Figure 9.

Results using KPLS-based GLRT method.

Figure 10.

Results using KPLS-based OEWMA method.

Figure 11.

Results using KPLS-based OEWMA-GLRT method.

Table 1.

Comparison in terms of MDR, FAR, and $\mathrm{AR}{L}_1$ .

Technique	GLRT			OEWMA			OEWMA-GLRT
	MDR	FAR	$\mathrm{AR}{L}_1$	MDR	FAR	$\mathrm{AR}{L}_1$	MDR	FAR	$\mathrm{AR}{L}_1$
PLS	52.277	51.253	37	32.574	55.03	1	30.28	49.002	24
KPLS	22.772	13.759	1	10.02	12.02	1	9.802	11.759	1

GLRT: generalized likelihood ratio test; OEWMA: optimized exponentially weighted moving average; OEWMA-GLRT: optimized exponentially weighted moving average-generalized likelihood ratio test; MDR: missed detection rate; FAR: false alarm rate; PLS: partial least squares; KPLS: kernel partial least squares.

Figures 6–8 show that the PLS-based OEWMA-GLRT chart provides better detection results comparing to both PLS-based OEWMA and PLS-based GLRT charts. Otherwise, the results show that the detection abilities based on linear PLS model present a lot of false alarm and MDRs. This is because the PLS method assumes a linear relationship between variables. However, the simulated synthetic example (29) is nonlinear, which means that the linear PLS method cannot tackle the issue of non-linearity. To deal with this issue, KPLS method is applied here for modeling purposes. Figures 9–11 and Table 1 show the detection performances based on KPLS model. The benefits of KPLS model on the detection performances over the linear PLS model are shown in Table 1.

Figures 9–11 illustrate that the developed KPLS-based OEWMA-GLRT method gives the best detection quality with respect to the classical KPLS-based GLRT and KPLS-based OEWMA charts. All of them outperform the detection charts based on PLS model (see Figures 6–9). Table 1 illustrates this in terms of the detection metrics: MDR, FAR, and $\mathrm{AR}{L}_1$ values.

Figures 12 and 13 and Table 2 show the effect of the multiscale representation on the detection abilities. In fact, the use of the multiscale representation of the data may enhance the monitoring quality by reducing the missed detection and the FARs. The wavelet-based multiscale representation of data may ameliorate the detection abilities of the KPLS- and PLS-based methods (Figures 8 and 11–13).

Figure 12.

Results using MSPLS-based OEWMA-GLRT method.

Figure 13.

Results using MSKPLS-based OEWMA-GLRT method.

Table 2.

Comparison in terms of MDR, FAR, and $\mathrm{AR}{L}_1$ .

Technique	MSPLS			MSKPLS
	MDR	FAR	$\mathrm{AR}{L}_1$	MDR	FAR	$\mathrm{AR}{L}_1$
OEWMA-GLRT	38.316	44.85	1	0	2.005	1

MSPLS: multiscale partial least squares; MSKPLS: multiscale kernel partial least squares; MDR: missed detection rate; FAR: false alarm rate; OEWMA-GLRT: optimized exponentially weighted moving average-generalized likelihood ratio test.

FD using simulated benchmark process

The validation of the proposed approach is also assessed using a nonlinear simulated benchmark structure.^36–39 The benchmark structures provided by the IASC-ASCE Structural Health Monitoring Group are studied to verify the effectiveness of damage detection.

The benchmark gives an analytical model based on the structure 12-degree-of-freedom (DOF) shear-building model. It consists of a $2.5\mathrm{m}\times 2.5\mathrm{m}$ plan with $3.6\mathrm{m}$ tall.^41,42 The sections are designed for a scale model. It contains four floors where the slabs and the beams move as rigid bodies with translation in the x and y directions and rotation $\theta$ about the center column. Hence, there are 3 degrees of freedom per floor. Figure 14 presents the analytical model diagram.

Figure 14.

Steel frame scale structure.

More details on the IASC-ASCE benchmark are presented in Chaabane et al. ⁴² This model was developed to generate simulated input–output response data. ⁴²

In order to construct the KPLS model, the input and output data should be scaled first. Then, the data set is divided into training data which consists of $1000$ samples and testing data which consists of $1000$ samples and that will be used to validate the performance of developed FD technique.

When we apply the multiscale KPLS to construct the model, the best decomposition depth needs to be chosen in order to get good detection. The best decomposition depth which gives the lower MDR, FAR, and $\mathrm{AR}{L}_1$ is equal to 3. To test the data, a mean shift fault of size $2\sigma$ is added to the third variable $X_3$ of the testing data set at $[800\mathrm{to}1000]$ .

The detection performances of GLRT-, OEWMA-, and OEWMA-GLRT-based linear and nonlinear PLS methods are shown in Figures 15–20 and Table 3. We can conclude from the results presented in Figures 15–17 and Table 3 that the PLS-based charts result in bad monitoring abilities. This is due to the fact that the PLS assumes that the relationships between variables are linear. Thus, it may not be appropriate for nonlinear analysis. Thus, in order to deal with nonlinear practical processes, a nonlinear PLS (KPLS) will be applied in the modeling phase. Figures 18–20 and Table 3 present the detection performances of KPLS-based charts. According to Figure 18 and Table 3, the detection efficiency of the KPLS-based OEWMA detection is lower than the KPLS-based OEWMA-GLRT (Figure 19 and Table 3). And all of them outperform the PLS-based approaches (Figures 15–17 and Table 3).

Figure 15.

Results using PLS-based OEWMA-GLRT method.

Figure 16.

Results using PLS-based OEWMA method.

Figure 17.

Results using PLS-based OEWMA-GLRT method.

Figure 18.

Results using KPLS-based OEWMA-GLRT method.

Figure 19.

Results using KPLS-based OEWMA method.

Figure 20.

Results using KPLS-based OEWMA-GLRT method.

Table 3.

MDR, FAR, and $\mathrm{AR}{L}_1$ evaluation using PLS and KPLS-based OEWMA-GLRT, OEWMA, and OEWMA-GLRT charts.

Technique	GLRT			OEWMA			OEWMA-GLRT
	MDR	FAR	$\mathrm{AR}{L}_1$	MDR	FAR	$\mathrm{AR}{L}_1$	MDR	FAR	$\mathrm{AR}{L}_1$
PLS	41.14	35.31	21	44.74	29.74	1	41.18	29.45	14
KPLS	15.25	6.29	1	4.52	91.82	1	13.72	5.47	1

GLRT: generalized likelihood ratio test; OEWMA: optimized exponentially weighted moving average; OEWMA-GLRT: optimized exponentially weighted moving average-generalized likelihood ratio test; MDR: missed detection rate; FAR: false alarm rate; PLS: partial least squares; KPLS: kernel partial least squares.

Figures 21 and 22 illustrate the detection performances and show the benefits of multiscale representation when using multiscale partial least squares (MSPLS)-based OEWMA-GLRT and MSKPLS-based OEWMA-GLRT methods. These figures show that MSKPLS-based OEWMA-GLRT method has superior effectiveness over linear MSPLS-based OEWMA-GLRT method (see Table 4) and both of them outperform the classical-based charts. This is because the multiscale representation of data might be effective to noises and has a great impact on the detection performances (Figures 21 and 22).

Figure 21.

Results using MSPLS-based OEWMA-GLRT method.

Figure 22.

Results using MSKPLS-based OEWMA-GLRT method.

Table 4.

MDR, FAR, and $\mathrm{AR}{L}_1$ evaluation using MSPLS and MSKPLS-based OEWMA-GLRT charts.

Technique	MSPLS			MSKPLS
	MDR	FAR	$\mathrm{AR}{L}_1$	MDR	FAR	$\mathrm{AR}{L}_1$
OEWMA-GLRT	9.116	1.45	1	0	1.35	1

MSPLS: multiscale partial least squares; MSKPLS: multiscale kernel partial least squares; MDR: missed detection rate; FAR: false alarm rate; OEWMA-GLRT: optimized exponentially weighted moving average-generalized likelihood ratio test.

Conclusion

In this work, a MSKPLS-based OEWMA and GLRT is proposed for FD. The MSKPLS method is applied for modeling and OEWMA-GLRT chart is used for detection. The developed OEWMA-GLRT consists of computing an optimal statistic that merges current and previous information providing more importance to the most recent data. It chooses the optimal tuning parameters that minimize the MDR, the FAR, and the out-of-control average run length $(\mathrm{AR}{L}_1)$ . This helps to provide a more accurate estimation of the GLRT statistic and provide a stronger memory that enables better decision making with respect to FD. The monitoring effectiveness of the proposed algorithm was compared to classical algorithms using two examples, synthetic data and benchmark structure. The results showed the detection efficiency of the proposed technique in terms of MDR, FAR, and $\mathrm{AR}{L}_1$ values.