A Novel Data-Driven Fault Detection Method Based on Stable Kernel Representation for Dynamic Systems

Abstract

With the steady improvement of advanced manufacturing processes and big data technologies, modern industrial systems have become large-scale. To enhance the sensitivity of fault detection (FD) and overcome the drawbacks of the centralized FD framework in dynamic systems, a new data-driven FD method based on Hellinger distance and subspace techniques is proposed for dynamic systems. Specifically, the proposed approach uses only system input/output data collected via sensor networks, and the distributed residual signals can be generated directly through the stable kernel representation of the process. Based on this, each sensor node can obtain the identical residual signal and test statistic through the average consensus algorithms. In addition, this paper integrates the Hellinger distance into the residual signal analysis for improving the FD performance. Finally, the effectiveness and accuracy of the proposed method have been verified in a real multiphase flow facility.

1. Introduction

As a result of intelligence and informatization, modern industrial processes have evolved to become more complicated. Any abnormal behavior of equipment components may affect productivity or even cause accidents. To guarantee the reliability and stability of industrial systems, fault detection (FD) plays a fundamental role and has received intensive attention from scholars and engineers [1,2,3].

Currently, the majority of FD methods are commonly classified as model-based techniques and data-driven techniques [4]. In the framework of model-based FD, it is necessary to obtain a precise mathematical representation of the systems. In addition, the accuracy of the FD results will depend on the accuracy of the modeling. According to the way of signal generation, existing model-based FD methods can be divided into three distinct groups, i.e., parameter estimation techniques, observer-based techniques, the subspace-based strategy [5]. In practical applications, such methods can successfully implement FD schemes when accurate mathematical models are available. Unfortunately, with the increasing size of modern engineering systems, the modeling of systems using first principles poses further challenges.

With the accelerated advancement of sensor networks and data processing methods, data-driven FD strategies have naturally become a critical research topic and are well developed [6,7,8,9,10]. Data-driven FD techniques are typically sorted into multivariate statistics, neural networks, and subspace-technique-aided schemes, etc. [11]. Specifically, traditional multivariate statistical methods [12,13,14], i.e., canonical correlation analysis (CCA), have been extensively studied and applied to modern industrial systems. The core of multivariate statistical methods is to analyze the correlation among process variables, followed by constructing appropriate test statistics for FD tasks. This group of methods can solve the FD problem well in static processes. However, such traditional FD strategies do not usually consider dynamic changes in the systems. Therefore, they are usually unable to perform FD tasks in dynamic systems. In the past few years, with increasing attention to statistical learning, neural network-based FD methods have been rapidly developed [15,16,17,18]. This group of methods usually uses historical data generated by the production process to train the neural network model. Then, the models classify the test data to determine whether faults have occurred. Due to the excellent fitting ability of neural networks, neural network-based methods have superior performance in dealing with FD problems for nonlinear systems. However, the training phase of neural networks requires the use of abundant labeled data, which has some limitations in practical applications [19]. In recent years, the subspace-technique-aided FD method has been widely studied because of its simple design and lower computational effort [20,21,22,23,24,25]. The core idea behind it is to identify the system parameters through the collected data. Because it takes into account the dynamic behavior in the systems, the subspace-aided FD method performs very well in dealing with process dynamics. In the subspace technique framework, Reference [23] proposes a stable kernel representation (SKR) of the systems. Remarkably, the proposed SKR scheme can directly construct residual generators using process data without identifying complex system models. Based on this, many data-driven FD algorithms are designed under the SKR framework.

At present, the majority of data-driven FD methods have been designed in a centralized framework, which involves collecting all process data in a central location to carry out necessary computations and FD tasks. With the growth of industrial process size, the centralized design becomes increasingly demanding in terms of memory and computing power, resulting in poor flexibility and high cost. Therefore, there is a strong research increase in distributed data-driven FD [26]. For example, ref. [27] evaluated the effectiveness of multiblock multivariate statistics-based approaches, such as PCA and PLS, for decentralized monitoring and assessed the individual contributions of each block. Although the above approaches can implement distributed FD for each sub-block, they do not take into account the connection among the sub-blocks. On this basis, ref. [28] designed a multiblock PCA strategy where information interaction among subblocks is considered. Similarly, considering the connection among neighboring nodes, ref. [29] designed a distributed CCA algorithm to achieve plant-wide process monitoring. The core idea behind it was to reduce uncertainty through information interaction among neighboring nodes. However, when there are few relevant variables in the historical data, the distributed FD results obtained based on this method are usually unreliable. In order to address this defect, ref. [30] designed a distributed, regularized CCA-based process monitoring algorithm. First, the traditional CCA algorithm is executed between the local and the neighboring nodes. In order to eliminate uncorrelated variables from the monitoring data, a GA-based regularization algorithm is then embedded in the traditional CCA technique. Finally, according to the local monitoring results, the corresponding residual signals and test statistics can be generated at the local node. In terms of technical systems, the multiphase flow facility used in this study is a device to achieve the separation of water, oil, and air at a given flow rate. The test zone, consisting of splitters and supply lines, can provide a mixture of water, oil, and air. The technical system is typically used to verify the effectiveness and accuracy of process monitoring and FD algorithms, such as the latent-variable-analysis-based FD method [31], the Kalman-filter-based FD method [32], and the multivariate-statistics-based process monitoring method [33]. In addition, the technical system used is described in detail in [34,35].

In general, the data-driven distributed FD methods above usually ignore the dynamic changes in the process. Therefore, these data-driven distributed methods have some limitations in dealing with FD problems in dynamic processes. In addition, there are abundant process data in dynamic systems, and the relationship among process variables is complicated. The abundant process data and strong coupling among variables bring new challenges to existing distributed data-driven FD solutions.

Motivated by the aforementioned points, a new data-driven distributed FD strategy is developed for dynamic systems. Compared with previous research, the key contributions of the developed FD solution are given as

• Compared with traditional SKR-based FD approaches, the proposed method is more sensitive to fault information by introducing the Hellinger distance (HD) in the residual signal.
• The consensus algorithm is embedded in the information interaction among sensor blocks. Therefore, each sensor block can obtain FD results without global fusion operations, thus remarkably improving FD efficiency.
• It has superior flexibility in the design of FD framework, particularly when the system models are not accurately obtained.

The structuring of this work is structured as follows. Section 2 provides information about the system descriptions, Hellinger distance, and the average consensus algorithms. In Section 3, a new distributed FD scheme for dynamic systems is presented. The effectiveness of the proposed FD algorithm is then demonstrated through a multiphase flow facility in Section 4. Finally, Section 5 presents the conclusion and prospects for future work.

2. Preliminaries

2.1. System Descriptions

Given a LTI system $H\left(z\right)$ with input factor $u\in R^u$ and output factor $y\in R^y$, the input–output (I/O) behavior is characterized as

$$\begin{array}{c}\hfill y\left(z\right)=H\left(z\right)u\left(z\right)\end{array}$$

(1)

where variable z represents the complex z-transform. In order to analyze the relationship among variables in dynamic systems, the state space model is used in this study. It not only reflects dynamic behavior in the process data but also provides a concise description, which is usually expressed in a standard form, as follows:

$$\begin{array}{c}\hfill x(k+1)=\left[\begin{array}{cc} A& B \end{array}\right]\left[\begin{array}{c}x\left(k\right)\\ u\left(k\right)\end{array}\right]\end{array}$$

(2)

$$\begin{array}{c}\hfill y\left(k\right)=\left[\begin{array}{cc} C& D \end{array}\right]\left[\begin{array}{c}x\left(k\right)\\ u\left(k\right)\end{array}\right]\end{array}$$

(3)

where A, B, C and D are system parameters; $u\left(k\right)\in R^u$, $y\left(k\right)\in R^y$, and $x\left(k\right)\in R^x$ refer to the system input, output, and state variables, respectively.

A sensor network consisting of $k_t$ blocks has been integrated into the considered system, as shown in Figure 1. In the sensor network, the network topology $\mathcal{G}$ can be represented using “node” $\mathcal{J}$ and “edge” $\mathcal{K}$ as

$$\begin{array}{c}\hfill \mathcal{G}=(\mathcal{J},\mathcal{K})\end{array}$$

(4)

where $\mathcal{J}$$=\{1,\cdots ,k_t\}$ depicts the set of sensor blocks; $\mathcal{K}\in \{\delta \times \delta \}$ depicts the set of edges.

In order to obtain the sef of adjacent blocks, $\mathcal{J}$ in (1) is further sketched as

$$\begin{array}{c}\hfill {\mathcal{J}}_i=\{i\in \mathcal{J}|(i,j)\in \mathcal{K}\}\end{array}$$

(5)

where ${\mathcal{J}}_i$ depicts all adjacent blocks of the i-th block.

2.2. Hellinger Distance

Hellinger distance (HD), also known as the Bhattacharyya distance, is a type of f-divergence [36]. The f-divergence is a function that measures the difference between two probability distributions. HD is a statistical technique that evaluates the resemblance of two frequency distributions to each other. Supposing that $n\left(x\right)$ and $m\left(x\right)$ denote two probability density functions (PDFs) and since the probability distributions of the variables are unknown in the general definition, the HD between $n\left(x\right)$ and $m\left(x\right)$ is next defined as

$$\begin{array}{c}\hfill H(n,m)=\sqrt{\frac{1}{2}{\int }_{-\infty }^{+\infty }{(\sqrt{n\left(x\right)}-\sqrt{m\left(x\right)})}^{2}dx}\end{array}$$

(6)

which can be also be expressed in Euclidean norm

$$\begin{array}{c}\hfill H(n,m)=\frac{1}{\sqrt{2}}\Vert \sqrt{n}-\sqrt{m}{\Vert }_2\end{array}$$

(7)

Based on the Cauchy inequality, HD is a symmetric bounded metric which satisfies $0\le H\left(n,m\right)\le 1$ and $H(n,m)=H(m,n)$. In addition, according to the Lebesgue theorem, the squared form of HD in (6) is characterized as

$$\begin{array}{c}\hfill H^2(n,m)=1-{\int }_{-\infty }^{+\infty }\sqrt{n\left(x\right)m\left(x\right)}dx\end{array}$$

(8)

In order to perform FD task, Lemma 1 gives a concise representation of (8), which serves as the basis of the proposed approach.

Lemma 1.

Consider two PDFs $n\left(x\right)\sim \mathcal{N}({\mu }_n,{{\sigma }^2}_n)$ and $m\left(x\right)\sim \mathcal{N}({\mu }_m,{{\sigma }^2}_m)$, $H^2(n,m)$ is further represented as follows:

$$\begin{array}{c}\hfill H^2(n,m)=1-\left(\frac{2{\sigma }_n{\sigma }_m}{{{\sigma }_n}^2+{{\sigma }_m}^2}\right)exp(-\frac{{({\mu }_n-{\mu }_m)}^2}{4({{\sigma }_n}^2+{{\sigma }_m}^2)})\end{array}$$

(9)

2.3. Average Consensus Algorithm

Given a communication network consisting of $k_t$ nodes, the consensus algorithm is a convergence technique to implement consensus calculations. For a data vector ${\mu }_i$ at the i-th sensor block, the consensus technique can be executed as

$$\begin{array}{c}\hfill {\mu }_i(s+1)=v_{i,i}{\mu }_i\left(s\right)+\sum _{j\in J_i}{v}_{i,j}{\mu }_j\left(s\right),i=1,\cdots ,k_t\end{array}$$

(10)

where ${\mu }_i\left(s\right)$ refers to the calculated value of ${\mu }_i$ during the sth iteration; $v_{i,j}$ is the weighting coefficients,. Many studies [37,38,39] have designed algorithms to solve the weighting problem. Among them, the Metropolis–Hastings technique not only speeds up the convergence of the iterative algorithm but also enables computation in a distributed manner. Therefore, the Metropolis–Hastings technique is used to construct the weighting factors in this paper. The weighting factors are assembled as follows

$$v_{ij}=\left\{\begin{array}{cc}\frac{1}{max\{j_i,j_j\}+1}& j\in J_i\\ 0& \hspace{1em}\hspace{1em}\hspace{1em}j\notin J_{i}andj\ne i\\ 1-{\displaystyle \sum _{k\ne i}}{v}_{ik}& i=j\end{array}\right.$$

(11)

where ${\mathcal{J}}_i$ represents all adjacent blocks of the i-th block. $j_i=\left|J_i\right|$ as the cardinality of $J_i$. Let V be

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}{v}_{1,1}& \cdots & v_{1,k_t}\\ ⋮& \ddots & ⋮\\ v_{k_t,1}& \cdots & v_{k_t,k_t}\end{array}\right]\in R^{k_t\times k_t}.\end{array}$$

(12)

The final consensus results can be presented as

$$\begin{array}{c}\hfill \underset{s\to \infty }{lim}\left[\begin{array}{c}{\mu }_1\left(s\right)\\ ⋮\\ {\mu }_{k_t}\left(s\right)\end{array}\right]=1\otimes (\frac{1}{k_t}\sum _{i=1}^{k_t}{\mu }_i\left(0\right))\end{array}$$

(13)

which indicates that the consensus value of each block will converge to the average value of all sensor blocks.

3. Methodology

In this section, the SKR of the residual generator is first introduced. Considering the probability distribution of fault information, a novel FD strategy is then presented and applied to the FD problem.

3.1. SKR

Considering the process model of (1)–(3) above, the left coprime factorizations (LCF) of $H\left(z\right)$ is given as follows:

$$\begin{array}{c}\hfill H\left(z\right)=D+C{(zI-A)}^{-1}B=P^{-1}\left(z\right)Q\left(z\right)\end{array}$$

(14)

where $\left(P\left(z\right),Q\left(z\right)\right)$ is called the left coprime pair. A key feature of the LCF under noise-free and fault-free conditions is displayed as

$$\begin{array}{c}\hfill \left[\begin{array}{cc}-P\left(z\right)& Q\left(z\right)\end{array}\right]\left[\begin{array}{c}u\left(z\right)\\ y\left(z\right)\end{array}\right]=0\end{array}$$

(15)

where $\left[\begin{array}{cc}-P\left(z\right)& Q\left(z\right)\end{array}\right]$ is denoted as the SKR of (2) and (3) [23]. Therefore, all LTI residual generators can be parameterized as

$$\begin{array}{c}\hfill r\left(z\right)=\left[\begin{array}{cc}-P\left(z\right)& Q\left(z\right)\end{array}\right]\left[\begin{array}{c}u\left(z\right)\\ y\left(z\right)\end{array}\right]=0\end{array}$$

(16)

3.2. Data-Driven Distributed Fault Detection

Given that a sensor network is integrated into the dynamic system, the state-space model of the system and sensor measurements with noise are displayed as

$$\begin{array}{c}\hfill x(k+1)=Ax\left(k\right)+Bu\left(k\right)+\sigma \left(k\right)\end{array}$$

(17)

$$\begin{array}{c}\hfill y_i\left(k\right)=\left\{\begin{array}{c}{L}_{i}x\left(k\right)+{\tau }_i\left(k\right), \mathrm{fault}-\mathrm{free}\hfill \\ L_i\left(x\left(k\right)+f\left(k\right)\right)+{\tau }_i\left(k\right), \mathrm{faulty}\hfill \end{array}\right.\end{array}$$

(18)

where $x\left(k\right)\in R^{k_l}$, $u\left(k\right)\in R^{k_m}$ denote the system state and the process input; $y_i\left(k\right)\in R^{k_n}$ is the output vector at the i-th sensor sub-block; $f\left(k\right)\in R^{k_f}$ denotes the unknown faults; $\sigma \left(k\right)\sim N(0,{\Sigma }_{\sigma })$ and ${\tau }_i\left(k\right)\sim N(0,{\Sigma }_{{\tau }_i})$ represent process and measurement noise, respectively. In addition, $\sigma \left(k\right)$ and ${\tau }_i\left(k\right)$ are assumed to be Gaussian distributions.

Considering that measurement data $y_i$, $i=1,\cdots k_t$, can be sent to a sensor block, a global model is then constructed as follows:

$$\begin{array}{c}\hfill y\left(k\right)=\left\{\begin{array}{c}Lx\left(k\right)+\tau \left(k\right), \mathrm{fault}-\mathrm{free}\hfill \\ L\left(x\left(k\right)+f\left(k\right)\right)+\tau \left(k\right), \mathrm{faulty}\hfill \end{array}\right.\end{array}$$

(19)

where

$$y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ ⋮\\ y_{k_t}\left(k\right)\end{array}\right], \mathrm{L}=\left[\begin{array}{c}{L}_1\\ ⋮\\ L_{k_t}\end{array}\right], \tau \left(k\right)=\left[\begin{array}{c}{\tau }_1\left(k\right)\\ ⋮\\ {\tau }_{k_t}\left(k\right)\end{array}\right].$$

To complete the algorithm implementation of SKR, data models are indispensable in the design processes [40]. Assuming there exists a data variable ${\kappa }_s\left(k\right)\in R^{k_{\kappa }}$, it can be further depicted as

$$\begin{array}{c}\hfill {\kappa }_s\left(k\right)={\left[\begin{array}{ccc}{\kappa }^T\left(k\right)& \cdots & {\kappa }^T(k+s)\end{array}\right]}^T\in R^{(s+1)k_{\kappa }}\end{array}$$

(20)

$$\begin{array}{c}\hfill {{\rm Y}}_k={\left[\begin{array}{ccc}\kappa \left(k\right)& \cdots & \kappa (k+N-1)\end{array}\right]}^T\in R^{k_{\kappa }\times N}\end{array}$$

(21)

$$\begin{array}{c}\hfill {{\rm Y}}_{k,s}={\left[\begin{array}{ccc}{\kappa }_s\left(k\right)& \cdots & {\kappa }_s(k+N-1)\end{array}\right]}^T\in R^{(s+1)k_{\kappa }\times N}\end{array}$$

(22)

where k denotes sampling instants, and s and N are some integers.

According to the extended models of $H\left(z\right)$ in (17) and (18), a data model is derived by iterative computation at each node:

$$\begin{array}{c}\hfill Y_{s,i}=F_{s,i}{X}_k+{\mathrm{G}}_{s,i}{U}_{k,s}+{\mathrm{H}}_{s,i}{\Theta }_{k,s}+{\Xi }_{s,i}\in R^{(s+1)k_n\times N}\end{array}$$

(23)

where

$$F_{s,i}=\left[\begin{array}{c}{L}_i\\ L_{i}A\\ ⋮\\ L_i{A}^s\end{array}\right],$$

$$G_{s,i}=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ L_{i}B& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ L_i{A}^{s-1}B& \cdots & L_{i}B& 0\end{array}\right],$$

$${\Xi }_{s,i}=\left[\begin{array}{ccc}{\tau }_i\left(k\right)& \cdots & {\tau }_i(k+N-1)\\ ⋮& \ddots & ⋮\\ {\tau }_i(k+s)& \cdots & {\tau }_i(k+s+N-1)\end{array}\right].$$

In order to remove the unobservable variable $X_k$, (23) is re-modeled as

$$\begin{array}{c}\hfill \left[\begin{array}{c}{U}_{k,s}\\ Y_{s,i}\end{array}\right]=K_{s,i}\left[\begin{array}{c}{U}_{k,s}\\ X_k\end{array}\right]+\left[\begin{array}{c}0\\ H_{s,i}{\Theta }_{k,s}+{\Xi }_{s,i}\end{array}\right]\end{array}$$

(24)

where $K_{s,i}=\left[\begin{array}{cc}I& 0\\ G_{s,i}& F_{s,i}\end{array}\right]\in R^{(s+1)(k_n+k_m)\times (n+(s+1)k_m)}$. When s is large enough, there must be a left nullspace of $K_{s,i}$ as

$${K_{s,i}}^{\perp }{K}_{s,i}=0,{K_{s,i}}^{\perp }\in R^{((s+1)k_n-k_l)\times (s+1)(k_n+k_m)}\Rightarrow$$

$$\begin{array}{c}\hfill {K_{s,i}}^{\perp }\left[\begin{array}{c}{U}_{k,s}\\ Y_{s,i}\end{array}\right]=\left[\begin{array}{cc}{K_{s,i,u}}^{\perp }& {K_{s,i,y}}^{\perp }\end{array}\right]\left[\begin{array}{c}{U}_{k,s}\\ Y_{s,i}\end{array}\right]=\left[\begin{array}{c}0\\ {K_{s,i}}^{\perp }(H_{s,i}{\Theta }_{k,s}+{\Xi }_{s,i})\end{array}\right].\end{array}$$

(25)

${K_{s,i}}^{\perp }$ is called the SKR of the system. Due to the excellent reliability and robustness of QR algorithm, the data-driven implementation of ${K_{s,i}}^{\perp }$ is able to be executed via QR decomposition and SVD:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{T}_{k-s_p,s_p-1,i}\\ U_{k,s}\\ Y_{s,i}\end{array}\right]=\left[\begin{array}{ccc}{R}_{1,1}& 0& 0\\ R_{2,1}& R_{2,2}& 0\\ R_{3,1}& R_{3,2}& R_{3,3}\end{array}\right]\left[\begin{array}{c}{Q}_{1,1}\\ Q_{2,1}\\ Q_{3,1}\end{array}\right]\end{array}$$

(26)

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{R}_{2,1}& R_{2,2}\\ R_{3,1}& R_{3,2}\end{array}\right]=\left[\begin{array}{cc}{V}_1& V_2\end{array}\right]\left[\begin{array}{cc}{\Psi }_{1,1}& 0\\ 0& {\Psi }_{2,2}\end{array}\right]\left[\begin{array}{c}{D_1}^T\\ {D_2}^T\end{array}\right]\end{array}$$

(27)

where $s_p$ represents the past moment; $T_{k-s_p,s_p-1,i}=\left[\begin{array}{c}{U}_{k-s_p,s_p-1}\\ Y_{k-s_p,s_p-1,i}\end{array}\right]$ is an I/O dataset of the past $s_p$ moments.

In addition, the noise terms can be identified by the proof in [20]:

$$\begin{array}{c}\hfill R_{3,3}{Q}_{3,1}=H_{s,i}{\Theta }_{k,s}+{\Xi }_{s,i}\end{array}$$

(28)

Observe that $\left[\begin{array}{cc}{R}_{2,1}& R_{2,2}\\ R_{3,1}& R_{3,2}\end{array}\right]$ and $K_{s,i}$ have same null space. It thus holds that

$$\begin{array}{c}\hfill {\Psi }_{2,2}\approx 0,{K_{s,i}}^{\perp }={V_2}^T\in R^{((s+1)k_n-k_l)\times (s+1)(k_n+km)}.\end{array}$$

(29)

In order to identify the residual signal, ${K_{s,i,y}}^{\perp }$ in (25) needs to be obtained in a data-driven manner. It has been demonstrated [40] that

$$\begin{array}{c}\hfill {K_{s,i,y}}^{\perp }{F}_{s,i}=0,{K_{s,i,u}}^{\perp }=-{K_{s,i,y}}^{\perp }{G}_{s,i}\end{array}$$

(30)

where ${K_{s,i,y}}^{\perp }={F_{s,i}}^{\perp }$. Therefore, the residual generator can be obtained as follows:

$$\begin{array}{c}\hfill {K_{s,i,y}}^{\perp }(H_{s,i}{\Theta }_{k,s}+{\Xi }_{s,i})={K_{s,i,y}}^{\perp }{R}_{3,3}{Q}_{3,1}\end{array}$$

(31)

$$\begin{array}{c}\hfill r_i\left(k\right)={K_{s,i,y}}^{\perp }{y}_s\left(k\right)-{K_{s,i,u}}^{\perp }{u}_s\left(k\right).\end{array}$$

(32)

Although the residual signal $r_i\left(k\right)$ generated by SKR has the advantages of having a simple design and a low computational effort, the robustness of its FD results often becomes weak under actual operating conditions. In order to improve the robustness of the SKR framework for FD applications, the probability distribution of the residual signal deserves further investigation. Based on the idea of HD, an approach to evaluate the similarity between two PDFs is introduced into the SKR framework.

Considering that the above process noise and measurement noise obey Gaussian distributions, the residual signal is $r_i\sim (0,{{\Lambda }_i}^2)$ in the normal (fault-free) historical dataset. In addition, for the actual fault dataset, the residual signal is ${\widehat{r}}_i\sim ({\widehat{\vartheta }}_i,{{\widehat{\Lambda }}_i}^2)$.

Therefore, a HD metric for the reisdual signal $r_i$ at the i-th sensor block can be represented as follows:

$$\begin{array}{c}\hfill {h_i}^2=1-{\int }_{-\infty }^{+\infty }\sqrt{f\left(r_i\right)f\left({\widehat{r}}_i\right)d{r}_i}\end{array}$$

(33)

where $f\left(r_i\right)$ and $f\left({\widehat{r}}_i\right)$ denote the PDFs of $r_i$ and ${\widehat{r}}_i$, respectively.

According to the property of Lemma 1, (33) is able to be rewritten as

$$\begin{array}{c}\hfill {h_i}^2=1-{\left(\frac{2{\Lambda }_i{\widehat{\Lambda }}_i}{{{\Lambda }_i}^2+{{\widehat{\Lambda }}_i}^2}\right)}^{1/2}exp\left(-\frac{{({\vartheta }_i-{\widehat{\vartheta }}_i)}^2}{4({{\Lambda }_i}^2+{{\widehat{\Lambda }}_i}^2)}\right)\end{array}$$

(34)

Based on the proposed FD algorithm, the $T^2$ statistic at each sensor block can be displayed as

$$\begin{array}{c}\hfill {T_i}^2=({{\widehat{h}}_i}^2-\widetilde{{h_i}^2}){\Psi }^{-1}{({{\widehat{h}}_i}^2-\widetilde{{h_i}^2})}^T\end{array}$$

(35)

where ${{\widehat{h}}_i}^2\in R^{k_n}$ can be obtained by (34) under the actual fault dataset; $\widetilde{{h_i}^2}\in R^{k_n}$ denotes the mean term of ${h_i}^2$ under the normal historical dataset; $\Psi$ is the covariance matrix of ${{\widehat{h}}_i}^2-\widetilde{{h_i}^2}$.

In order to implement distributed FD, each sensor node needs to perform the identical $T^2$ test statistic. Based on the above purpose, the average consensus technique is introduced in this framework. The consensus algorithm for ${\varpi }_i={{\widehat{h}}_i}^2-\widetilde{{h_i}^2}$ is

$$\begin{array}{c}\hfill {\varpi }_i(s+1)=v_{i,i}{\varpi }_i\left(s\right)+\sum _{j\in J_i}{v}_{i,j}{\varpi }_j\left(s\right),\mathrm{i}=1,\cdots ,{\mathrm{k}}_t\end{array}$$

(36)

where s denotes the iteration number. The initial value is ${\varpi }_i\left(0\right)={\varpi }_i$. Furthermore, it holds that [5]

$$\begin{array}{c}\hfill \underset{s\to \infty }{lim}{\varpi }_i(s+1)=\overline{\varpi }=\frac{1}{k_t}\sum _{i=1}^{k_t}{\varpi }_i=\frac{1}{k_t}\sum _{i=1}^{k_t}\left({{\widehat{h}}_i}^2-\widetilde{{h_i}^2}\right)\end{array}$$

(37)

As the algorithm runs until convergence, the consensus result can be obtained at each block:

$$\begin{array}{c}\hfill {\varpi }_i=\frac{1}{k_n}\sum _{i=1}^{k_n}({{\widehat{h}}_i}^2-\widetilde{{h_i}^2})\end{array}$$

(38)

Based on the consensus techniques above, identical ${\varpi }_i$ is obtained at each node. Therefore, (35) can be rewritten as

$$\begin{array}{c}\hfill {T_i}^2={\varpi }_i{\varphi }^{-1}{{\varpi }_i}^T\end{array}$$

(39)

where $\varphi$ denotes the the covariance matrix of ${\varpi }_i$.

As a result, the $T^2$ statistic can be executed in parallel at each block. When the amount of data is sufficient, the used $T^2$ statistic obeys a chi-square distribution ($T^2\sim {{\chi }_{\beta }}^2\left(k_n\right)$). Specifically, ${\chi }_{\beta }$ can be determined by a ${\chi }^2$ distribution with degrees $k_n$ of freedom as follows:

$$\begin{array}{c}\hfill prob\left\{\chi >{\chi }_{\beta }\right\}=\beta \iff prob\left\{\chi \le {\chi }_{\beta }\right\}=1-\beta \end{array}$$

(40)

Based on this, the fault detection threshold can be calculated as

$$\begin{array}{c}\hfill J_{th}={{\chi }_{\beta }}^2\left(k_n\right)\end{array}$$

(41)

where $k_n$ is the dimension of the residual data; and $\beta$ is the confidence level (acceptable false alarm rate).

In addition, the FD logic for each node is represented as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{T_i}^2<J_{th}\Rightarrow \mathrm{fault}-\mathrm{free},\hfill \\ {T_i}^2\ge J_{th}\Rightarrow \mathrm{faulty}.\hfill \end{array}\right.\end{array}$$

(42)

In summary, with the help of the SKR framework, Hellinger distance, and average consensus algorithm, the distributed FD scheme is summarized in Algorithms 1 and 2. In addition, the flow chart of the proposed FD algorithm is shown in Figure 2.

Algorithm 1: Off-Line Phase.

S1. Load the normal (fault-free) data.; S2: Set two indices $i,j$ where $i=1,\cdots ,k_t$ and $j\in J_i\cup {\overline{J}}_i$; S3: while $i\le k_t$ do S4: Constuct I/O data model at each node via (23); S5: Perform QR-decomposition (26) and SVD (27); S6: Identify SKR ${K_{s,i}}^{\perp }$ via (29) and (30); S7: Obtain the residual signals $r_i$ at each sensor node; S8: Calculate Hellinger distance ${h_i}^2$ for each residual signal via (34); S9: end while S10: Constuct weight matrix V via (11).

Algorithm 2: Online Phase.

S1. Load the actual test data.; S2: while $i\le k_t$ do S3: Obtain the residual signals ${\widehat{r}}_i$ under the acutal fault case; S4: Calculate Hellinger distance ${{\widehat{h}}_i}^2$ for each residual signal via (34); S5: end while S6: Calculate $\overline{\varpi }=\frac{1}{k_t}{\displaystyle \sum _{i=1}^{k_t}}\tilde{({{\widehat{h}}_i}^2-{h_i}^2})$ using the average consensus algorithm (36); S7: Obtain the identical $T^2$ test statistic via (39); S8: Make a FD logic decision whether a fault has occurred based on (41).

Figure 2. The flow chart of the proposed fault detection algorithm.

Figure 2. The flow chart of the proposed fault detection algorithm.

4. Case Study

4.1. Facility Description

In this study, a data set obtained from a multi-phase flow plant [35] is used to validate the proposed FD algorithm. The multiphase flow plant can achieve gas–liquid separation at a given flow rate. This device takes into account various working conditions during operation, so it can generate abundant process data from different operating conditions. In addition, the generated process data contain dynamic behavior by changing the set point of the flow rate. It is depicted in Figure 3, and its schematic diagram is presented in Figure 4. Specifically, the device comprises geometrically designed pipes and a $1.2$ m high liquid–gas splitter. It is capable of providing separate air, oil, and water, as well as mixtures of these fluids. During the operation of the plant, the mixtures are split in a horizontal splitter. The air is returned to the environment, while the water–oil mixture is returned to their respective tanks (T100 and T200). The water coalescers ensure complete separation of oil and water before returning to their respective tanks. The flow conditions of air, water, and oil can be regulated by control valves. In addition, the relevant control valves can be operated continuously between closing and opening. In terms of sensor distribution, there are sensors measuring pressure at the air delivery line (PT417) and inside the three-phase splitter (PT501). Other sensors are located at the water delivery line (FIC101, FT102), at the bottom of the two-phase splitter(FT406), at the top of the two-phase splitter (FT404), at the top of the three-phase splitter (PIC501), at the air delivery line (FIC302, FT302), at the top of the water tank (LI101), at the bottom of the three-phase splitter (LI502), and at the top of the water coalescer (LI503).

4.2. Fault Injection and Distributed Fault Detection

In order to gather the necessary historical data for this experiment, a SCADA platform can be utilized at a sample rate of 1 Hz. The data parameters utilized in this validation are outlined in

Table 1. Parameters used in this experiment.

Sensor Location	Description	Units
PIC501	Placement of PIC501 valve	(%)
PT417	Pressure in the blending area	barg
FIC302	Placement of FIC302 valve	(%)
FIC101	Placement of FIC101 valve	(%)
FT404	Air flow rate from 2-phase splitter	m${}^3$/h
FT406	Water flow rate from 2-phase splitter	kg/s
PT501	Pressure in 3-phase splitter	barg
LI101	Level of water tank	m
LI502	Level of 3-phase splitter	(%)
LI503	Level of watercoalescer	(%)
FT302	Air intake velocity	Sm${}^3$/h
FT102	Water intake velocity	kg/s

. A communication topology of the sensor network is represented in Figure 5.

To evaluate the effectiveness of the proposed method, the off-line part is first executed. The input flow rate of the training dataset is shown in Figure 6. Then, two typical faults are used to verify FD performance. The first fault scenario involves an incipient fault that arises due to the obstruction of the top separator input, leading to the shutdown of VC404 between 1136 s and 8352 s. The input situation of the faulty dataset is depicted in Figure 7. The residual signal of fault 1 and the HD-based $h^2$ statistic at node 1 are shown in Figure 8. The residual signal $r_1\left(k\right)$ in Figure 8 indicates the overall trend of the system. However, when the fault amplitude is low, the residual signal is often not effective in capturing abnormal cases. The blue curve in Figure 8 does not changed significantly after the fault occurred.

Based on the above problem, a HD-based metric is implemented into the proposed distributed SKR framework to enhance the sensitivity of the residual signal. When the fault intensity is minor, the information on the probability distribution also changes remarkably. As a result, the sensitivity of FD is significantly increased by measuring the HD of the residual signal. The green curve in Figure 8 indicates that the HD-based $h^2$ statistics have changed significantly after the fault occurred. The performance evaluation focuses on two key metrics: the Missed Alarm Rate (MAR) and the False Alarm Rate (FAR). In this study, the consensus algorithm is embedded in the information interaction among sensor blocks. Each node in the sensor network is used to execute FD algorithm through average consensus techniques. As a consensus result, each sensor node can obtain an identical FD performance. The distributed FD diagrams for fault 1 are displayed in Figure 9 and Figure 10. J(2) and J(8) in Figure 9 and Figure 10 represent the distributed FD results at node 2 and node 8, respectively. In terms of performance metrics, the MAR is 0.0323 and the FAR is 0.0391 at node 2. In addition, the MAR is also 0.0323 and the FAR is also 0.0391 at node 8. The second fault scenario is an intermittent fault, also known as segment plugging in practical engineering terms. This type of fault commonly occurs in the riser of multiphase flow when the flow rate of liquid and gas is low. The fault was introduced by deliberately reducing the air and water flow rates to levels that induce plugging. In this dataset, two plugging conditions were introduced and eliminated, from 686 s to 1172 s and from 1772 s to 2253 s, during the experiment. Specifically, the plugging fault was first formed at 686 s by continuously reducing the flow rate of air and water. The plugging fault was then removed at 1172 s when the air flow rate gradually increases. Additionally, the plugging fault was introduced again from 1772 s to 2253 s by changing the input flow rate. The input situation and detection results for fault 2 are shown in Figure 11, Figure 12 and Figure 13, respectively. The calculated MAR and FAR at node 1 and node 5 are 0.0610 and 0.0454, respectively.

4.3. Comparison Results

In order to show the enhanced FD performance,

Table 2. Performance comparisons for the multiphase flow plant.

Detection Framework	FD Strategies	${\mathit{f}}_1:$ Blockage in the Input Block Separator			${\mathit{f}}_2:$ Slugging Situation		Obtain FD Result’s Way
		MAR	FAR		MAR	FAR
Centralized	Traditional SKR [40]	0.0452	0.0729		0.5036	0.1693	Central node
	Dynamic principal component analysis [36]	0.6417	0.3461		0.5826	0.2754	Central node
Distributed	The developed HSKR	0.0323	0.0391		0.0610	0.0454	Any node
	Distributed CCA [29]	0.3881	0.2668		0.2076	0.3045	Any node

provides four sets of evaluations using MARs and FARs as assessment indicators. In Table 2, traditional SKR and dynamic principal component analysis are centralized designs; The distributed CCA and the proposed scheme are distributed designs.

According to the performance indicators in Table 2, both FAR and MAR of the proposed algorithm are significantly lower than other FD algorithms. This indicates that the accuracy and effectiveness of the proposed algorithm are better than the traditional SKR algorithm and other comparison algorithms. The excellent FD performance is mainly the result of the introduction of the Hellinger distance in the SKR framework. Specifically, the Hellinger distance is first introduced into the traditional SKR framework to further analyze the fault features of the residual signal. Since the Hellinger distance can accurately measure the difference between two probability distributions, the proposed algorithm is more sensitive to the fault information of the residual signal. As a result, the proposed algorithm can capture the fault information in the residual signal more effectively, which further improves the reliability and accuracy of the FD algorithm.

5. Conclusions

This study proposes a novel distributed FD method by introducing Hellinger distance and average consensus algorithm in the SKR framework. The proposed algorithm has the following three main advantages and differences over existing FD methods. This study introduces the first Hellinger distance in the traditional SKR framework to further analyze the fault features of the residual signal. Since the Hellinger distance can accurately measure the difference between two probability distributions, the proposed algorithm is more sensitive to the fault information of the residual signal. In addition, the consensus algorithm is embedded in the information interaction among sensor blocks. Based on this idea, each block can obtain FD results without performing global fusion operations. Finally, the proposed algorithm can identify noise terms and the residual signals directly from the process data. It has superior flexibility in the design of the detection framework, particularly when the system models are not accurately obtained. The accuracy and validity of the proposed FD algorithm have been verified via a multiphase flow facility. In addition, fault-tolerant control based on data-driven SKR is an open problem that can avoid the complex design of control systems. Based on this study, distributed FD with external disturbances and fault-tolerant control based on data-driven SKR will be explored in our future work.