一种新颖的深度因果图建模及其故障诊断方法

现代工业过程是一个高度复杂的系统. 由于实际的物理连接、控制回路的作用, 工业过程中的设备、部件、过程变量相互耦合, 构成了复杂的互连网络. 这种互联耦合关系使得系统某一部位一旦发生异常, 将会随着系统网络传播并演变演化, 进而导致更加严重的故障. 采用先进的故障检测和诊断技术是保证工业过程安全有效稳定运行的重要手段. 传统的基于知识或者模型的方法很难构建大规模系统变量间的复杂关系. 随着工业自动化、信息化的快速发展, 工业过程收集了越来越多的传感器数据, 这为数据驱动方法提供强有力的支撑. 数据驱动的故障检测和诊断方法也因此受到了广泛关注和研究, 并大量应用到化工、半导体制造等过程, 尤其是多元统计过程监测 (Multivariate statistical process monitoring, MSPM) 方法[1].

MSPM利用多元投影技术将高维观测数据投影到低维主元子空间和残差子空间, 并设计相应的多元统计量 (如T2$T^2$, SPE$SPE$) 及其控制限来监测数据是否超出正常工作范围, 常用的多元统计方法有主元分析 (Principal component analysis, PCA)[2]、偏最小二乘 (Partial least squares, PLS)[3]、独立主元分析 (Independent component analysis, ICA)[4] 和典型相关分析(Canonical correlation analysis, CCA)[5] 等. 一旦检测出故障, 需要采用贡献图[6]、重构贡献图[7] 等故障隔离方法来辨识故障相关变量. 这些方法由于拖尾效应[8]的影响并不能准确地辨识出所有的故障变量. 随后, 格兰杰因果分析[9]、传递熵[10] 等因果分析方法被用来进行故障根源诊断和传播路径辨识, 然而这些算法很难达到预期效果, 且需要较长的分析时间. 深度学习能够处理大规模数据, 自动地提取深层次特征, 在工业领域也取得了较成功的应用. 近些年来, 自动编码器 (Autoencoder, AE)[11]、深度信念网(Deep belief network, DBN)[12]、变分自编码器(Variational autoencoder, VAE)[13] 等深度学习技术越来越多地应用到过程监测中, 但是由于神经网络的黑箱特性, 变量贡献率难以计算, 限制了其在故障隔离、根源诊断和传播路径辨识方面的应用.

图论技术利用由若干节点和连接节点的线构成的图模型来定性或者定量地表征变量之间关系, 能够较好地描述工业过程的复杂网络结构[14-15]. 其中贝叶斯网络 (Bayesian network, BN) 作为一种概率有向图模型, 已经应用在风险分析、可靠性可维护性分析等领域[16], 在故障检测和诊断领域也取得了较成功的应用. Mehranbod等使用BN进行稳定和过渡阶段的故障检测和隔离[17-18]. Azhdari和Mehranbod验证了BN在田纳西−伊斯曼 (Tennessee Eastman, TE) 过程中故障检测与诊断应用的有效性[19]. Gonzalez等将BN应用在过程监测中的维数化简过程[20]. Chen和Ge提出了一个分层BN (Hierachical Bayesian network, HBN) 建模框架[21], 通过对工业过程进行分解, 构建局部单元以及单元间的贝叶斯网络, 实现大规模工业过程的故障检测与诊断. 随后, Chen和Ge又针对过程监测中低质量数据建模问题, 提出了鲁棒BN (Robust Bayesian networks, RBN)[22]. 针对工业过程中存在的动态特性, Yu和Rashid采用动态贝叶斯网络 (Dynamic Bayesian networks, DBN), 用异常似然指标 (Abnormality likelihood index, ALI) 和动态贝叶斯概率指标 (Dynamic Bayesian probability index, DBPI) 分别进行故障检测和根源诊断[23]. Zhang和Dong将高斯混合模型(Gaussian mixture model, GMM)和三时间切片DBN整合来解决过程监测中存在的数据缺失和非高斯问题[24].

BN为过程变量的因果关系提供了条件概率表示, 但是, 该条件概率关系一般是线性的, 无法描述过程中存在的非线性特性. 同时, 大多数基于BN的方法需要通过过程知识得到BN的网络结构, 这在很多复杂工业过程中是很难确定的. BN作为一种有向无环图模型, 也较难表示系统中存在的闭环结构. 为解决上述问题, 并考虑系统中存在的动态特性, 本文提出了一种深度因果图 (Deep causality graph, DCG) 建模方法, 利用多层感知器 (Multilayer perceptron, MLP) 和门控循环单元 (Gate recurrent unit, GRU) 对每一个过程变量建立概率预测模型. 在模型训练过程中, 引入组稀疏惩罚项, 自动地检测变量间因果关系, 从而得到过程变量的因果有向图结构以及定量的条件概率表征; 然后基于DCG模型的条件后验概率分布建立单变量监测统计指标, 并通过贝叶斯推理融合, 构建综合的监测统计量, 实现工业过程的整体监测. 进一步通过计算变量贡献度指标, 隔离出故障相关变量; 最后根据深度因果图模型获得的有向图结构, 诊断出故障根源, 并辨识故障的传播路径.

论文的结构如下: 第1$1$节详细介绍了深度因果图推导和建模过程; 第2$2$节提出了基于深度因果图模型的故障检测和诊断框架; 随后, 在第3$3$节利用TE过程数据对所提算法进行了验证, 并在最后一节中进行了总结.

1.   深度因果图建模方法

1.1   图结构已知的因果关系建模方法

考虑一个有向图模型, 其网络结构为G=⟨V,E⟩$G=⟨V,E⟩$, V$V$表示节点的集合, E$E$表示连接节点的有向边集合. 在工业过程中, 节点可以指过程变量特征, 有向边对应着连接节点之间的因果关系, 有向边的首和尾连接的节点分别表示为父节点和子节点.

对一个有n$n$维观测变量的工业过程, t$t$时刻观测变量表示为xt=[x1,t,x2,t,⋯,xn,t]${\mathit{x}}_t=[x_{1,t},x_{2,t},\cdots ,x_{n,t}]$, 则节点i$i$在t$t$时刻观测变量为xi,t$x_{i,t}$, 其父节点集在t$t$时刻的状态用xpa(i),t${\mathit{x}}_{pa(i),t}$表示. 考虑系统的动态特性, 观测变量xi,t$x_{i,t}$不仅与当前时刻父节点状态有关, 而且与历史时刻自节点和父节点的状态有关. t$t$时刻观测变量xi,t$x_{i,t}$用节点i$i$的历史观测数据xi,t−T:t−1${\mathit{x}}_{i,t-T:t-1}$和其父节点集的观测数据xpa(i),t−T:t${\mathit{x}}_{pa(i),t-T:t}$非线性表示:

式中, εi,t${\epsilon }_{i,t}$为随机噪声项. 利用概率形式来表示这种变量间的依赖关系:

式中, xc,t−T:t${\mathit{x}}_{c,t-T:t}$为xpa(i),t−T:t${\mathit{x}}_{pa(i),t-T:t}$和xi,t−T:t−1${\mathit{x}}_{i,t-T:t-1}$组合的观测值向量; μi(xc,t−T:t)${\mu }_i({\mathit{x}}_{c,t-T:t})$和σ2i(xc,t−T:t)${\sigma }_{{}_i}^2({\mathit{x}}_{c,t-T:t})$分别为xi,t$x_{i,t}$的后验概率分布的均值和方差, 均为非线性函数, 其非线性关系可以用神经网络表示.

通过神经网络构建节点i$i$预测模型, 其条件概率对数似然表示为lnp(xi,t|xc,t−T:t)$\ln p(x_{i,t}|{\mathit{x}}_{c,t-T:t})$. 设预测模型参数集为Θi${\mathrm{\Theta }}_i$, 模型参数Θi${\mathrm{\Theta }}_i$学习的目标是最大化对数似然期望值, 目标函数表示为:

由上述讨论可知观测变量xi,t$x_{i,t}$由节点i$i$与其父节点的历史状态以及父节点的t$t$时刻状态共同决定. 将节点i$i$及其父节点的历史状态中与变量xi,t$x_{i,t}$相关的动态特征信息表示为zi,t−1${\mathit{z}}_{i,t-1}$, 父节点当前时刻的相关特征信息表示为hi,t${\mathit{h}}_{i,t}$, 式 (1) 中的非线性预测模型于是可以转变为以下形式:

当前时刻相关特征hi,t${\mathit{h}}_{i,t}$的信息提取可以利用多层感知器实现, 历史时刻的动态特征信息zi,t−1${\mathit{z}}_{i,t-1}$采用GRU进行特征提取. GRU作为循环神经网络 (Recurrent neural network, RNN) 的一种变体能够提取时间序列的动态特征, 同时有效地解决标准RNN算法存在的长期记忆和反向传播中的梯度问题. GRU的详细介绍可以参考文献[25].

1.2   图结构未知的因果图建模方法

上述节点预测模型是建立在因果有向图结构已知的情况. 而在很多复杂工业过程中, 变量之间的因果关系是未知的. 为了寻找每个变量的因变量, 本文首先将除自变量外的其他变量都默认为因变量, 然后在模型训练过程中引入Group Lasso惩罚项, 使得节点预测模型中与输入变量相关的连接稀疏化, 使用尽可能少的输入获取尽可能准确的预测, 进而实现变量之间的因果关系自动检测, 其单一节点i$i$的后验概率预测模型如图1所示. 节点的概率分布预测模型由两个动态特征提取单元和一个预测输出单元构成. 在第一个动态特征提取单元中, 输入为t$t$时刻观测变量值xt={xi,t,xi−,t}${\mathit{x}}_t=\{x_{i,t},{\mathit{x}}_{i-,t}\}$, xi−,t${\mathit{x}}_{i-,t}$表示除节点i$i$之外其他节点的观测值. h1i−,t∈Rm1${\mathit{h}}_{{}_{i-,t}}^1\in {\mathbf{R}}^{m^1}$为在t$t$时刻从输入xi−,t${\mathit{x}}_{i-,t}$提取的m1$m^1$维隐藏特征, 其过程表示为:

图 1  深度因果图的单节点预测模型网络结构

Fig. 1  The network structure of single node prediction model for deep causality graph

式中, σ$\sigma$为激活函数; W1${\mathit{W}}^1$表示第一层全连接层的连接权重, 是一个m1×(n−1)$m^1\times (n-1)$大小的权重矩阵. 将W1${\mathit{W}}^1$展开为{w1i,1,⋯,w1i,i−1,w1i,i+1,⋯,w1i,n}$\{{\mathit{w}}_{i,1}^1,\cdots ,{\mathit{w}}_{i,i-1}^1,{\mathit{w}}_{i,i+1}^1,\cdots ,{\mathit{w}}_{i,n}^1\}$, 则式(5) 重写为:

将提取的隐藏特征h1i−,t${\mathit{h}}_{i-,t}^1$和观测变量xi,t$x_{i,t}$合并, 作为GRU单元的输入, 可得到t$t$时刻的动态特征z1i,t${\mathit{z}}_{i,t}^1$.

为了获取深层次的特征, 增大节点特征提取的视野域, 以便用较少的变量连接获取准确的预测结果, 将除自节点外的其他节点提取的动态特征合并作为第二个动态特征提取单元的输入, 再进一步地提取节点的动态特征. 其网络结构与第一个动态特征提取单元一致, 第二层全连接层连接权重为W2={W2i,1,⋯,W2i,i−1,W2i,i+1,⋯,W2i,n}${\mathit{W}}^2=\{{\mathit{W}}_{i,1}^2,\cdots ,{\mathit{W}}_{i,i-1}^2,{\mathit{W}}_{i,i+1}^2,\cdots ,{\mathit{W}}_{i,n}^2\}$, 其中W2∈Rm2×n−1${\mathit{W}}^2\in {\mathbf{R}}^{m^2\times n-1}$.

预测输出单元将t−1$t-1$时刻的动态特征z1i,t−1${\mathit{z}}_{i,t-1}^1$和z2i,t−1${\mathit{z}}_{i,t-1}^2$以及t$t$时刻隐藏特征h1i−,t${\mathit{h}}_{i-,t}^1$合并作为输入, 通过多层感知器输出节点i$i$在t$t$时刻观测变量预测概率分布的均值μi,t${\mu }_{i,t}$和方差对数ln(σ2i,t)$\ln ({\sigma }_{i,t}^2)$.

组最小绝对收缩和选择算子(Group least absolute shrinkage and selection operator, Group Lasso) 惩罚项被选择为稀疏惩罚项加入到单节点的概率预测模型目标函数中. Group Lasso是Yuan在2006$2006$年在Lasso稀疏约束方法上到组上的推广[26]. 它根据输入变量对权重矩阵进行分组, 然后在目标函数中惩罚每一组的L2范数, 这样就可以将部分组的全部系数同时消成零, 即抹掉整组的变量, 这种方法叫做Group Lasso分组最小角回归算法. 对节点i$i$的预测模型网络权重矩阵W1${\mathit{W}}^1$和W2${\mathit{W}}^2$, 根据连接的过程变量节点划分成n−1$n-1$组. 与过程变量节点j$j$对应的权重范数为∥∥w1i,j∥∥22+∥∥W2i,j∥∥2F−−−−−−−−−−−−−−√$\sqrt{{\Vert {\mathit{w}}_{i,j}^1\Vert }_2^2+{\Vert {\mathit{W}}_{i,j}^2\Vert }_{\mathrm{F}}^2}$. 加入Group Lasso惩罚项后, 节点i$i$的预测模型通过最小化目标函数(7) 来迭代更新网络参数:

超参数λ$\lambda$控制着连接权的稀疏性, λ$\lambda$越大, 越多组的权重会被消零, 使其节点的因果关系连接呈现为稀疏集.

由于在单一节点预测过程中需要利用其他节点提出的动态特征信息, 因此在模型训练过程中需将n$n$个过程变量节点的预测模型作为整体进行训练, 故深度因果图整体模型的目标函数是n$n$个目标函数的和形式:

在深度因果图模型训练过程中, 首先要用大小为T+1$T+1$的滑动窗口将标准化后的时间序列训练数据集规范为Xtrain∈RNtrain×(T+1)×n${\mathit{X}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\in {\mathbf{R}}^{N_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\times (T+1)\times n}$. 由于目标函数中存在Group Lasso惩罚项, 其在权重范数为0$0$处不可微. 为解决目标函数存在不可微的问题, 采用了近端梯度下降 (Proximal gradient descent, PGD) 算法进行模型优化训练. 对包含可微和不可微函数的目标函数f(θ)$f(\theta )$, 优化问题分解为:

式中, g(θ)$g(\theta )$是可微函数, h(θ)$h(\theta )$是不可微函数.

定义一个近端映射函数:

参数θ$\theta$的迭代过程可以表示为:

对给定的一个批次的训练数据, PGD的一次迭代算法可以理解为: 给定优化参数起点θ(k−1)${\theta }^{(k-1)}$, 可微函数g$g$沿着起点的负梯度方向, 做步长为tk$t_k$的梯度下降得到一个预更新值, 然后使用近端映射寻找一个z$z$, 使得不可微函数h(θ)$h(\theta )$足够小, 并接近预更新值, 用z$z$作为本次迭代的更新值θ(k)${\theta }^{(k)}$.

在模型训练中, 选择软阈值算子作为PGD的近端算子. 该优化算法可以使非因果变量的连接权重精确地收敛到0$0$, 以便能够用权重值来解释变量之间的因果性.

2.   基于深度因果图模型的故障检测与诊断

2.1   故障检测

完成深度因果图模型构建后, 可以利用该图模型得到的条件后验概率分布设计过程监测统计量. 首先, 对每一个变量节点设计单变量统计指标. 考虑观测变量为n$n$维的动态连续过程, t$t$时刻观测变量表示为xx1:n,t={x1,t,x2,t,⋯,xn,t}${xx}_{1:n,t}=\{x_{1,t},x_{2,t},\cdots ,x_{n,t}\}$. 节点i$i$有实际观测值xi,t$x_{i,t}$. 根据深度因果图模型, 计算其预测值概率分布为x~i,t∼N(μi,t,σ2i,t)${\tilde{x}}_{i,t}\sim \mathrm{N}({\mu }_{i,t},{\sigma }_{i,t}^2)$, 其平方马氏距离D2i,t$D_{i,t}^2$服从自由度为1$1$的χ2${\chi }^2$分布:

D2i,t$D_{i,t}^2$的控制限可以根据用户指定的置信度α$\alpha$来确定:

当D2i,t>D2lim$D_{i,t}^2>D_{lim}^2$时, 认为观测变量xi,t$x_{i,t}$为异常状态.

为了获取一个综合评价指标来监测整个过程, 本文利用贝叶斯融合对各节点检测信息进行融合. 根据贝叶斯推理, 变量节点i$i$的统计量D2i,t$D_{i,t}^2$为故障的概率表示为:

式中, PD2i(False)$P_{D_i^2}(\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e})$表示为统计量D2i,t$D_{i,t}^2$异常状态的先验概率, 在本文中PD2(False)=1−α$P_{D^2}(\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e})=1-\alpha$, PD2(True)=α$P_{D^2}(\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e})=\alpha$. 根据贝叶斯规则, PD2i(xi,t)$P_{D_i^2}(x_{i,t})$用以下公式计算:

式中, PD2i(xi,t|False)$P_{D_i^2}(x_{i,t}|\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e})$和PD2i(xi,t|True)$P_{D_i^2}(x_{i,t}|\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e})$设定成如下形式:

调节因子ηp${\eta }_p$可以调整概率值对故障状态的灵敏度, 在本文中ηp${\eta }_p$设置为1$1$.

根据上述推导过程, 可以将单节点的监测统计指标转换为概率形式PD2i(False|xi,t)$P_{D_i^2}(\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}|x_{i,t})$表示. 当PD2i(False|xi,t)>1−α$P_{D_i^2}(\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}|x_{i,t})>1-\alpha$时, 节点i$i$发生异常. 然后对单变量概率监测指标加权平均, 来计算综合监测指标PD2${PD}^2$: