CN109359662A - 一种面向百万千瓦超超临界机组非平稳特性的基于因果分析的多层贝叶斯网络故障诊断方法 - Google Patents

Abstract

本发明公开了一种面向百万千瓦超超临界机组非平稳特性的基于因果分析的多层贝叶斯网络故障诊断方法。本发明首先运用了稀疏协整分析和格兰杰因果分析分别对非平稳变量提取局部的精细化的非平稳变量间的因果关系子网，同时对所有变量提取全局性的因果关系。最后，结合贝叶斯网络对故障路径进行追溯并定位根源故障变量。该方法克服了大型燃煤机组因果分析中非平稳性的影响，充分挖掘了非平稳变量间所包含的潜在因果关系，建立了多层诊断模型，有效地解决了复杂非平稳故障过程定位根源故障变量困难的问题，大大提高了非平稳过程的故障传播的诊断性能，有助于火电厂对工厂进行有效及时的排查检修，从而保证了百万千瓦超超临界发电机组的安全可靠运行。

Description

一种面向百万千瓦超超临界机组非平稳特性的基于因果分析 的多层贝叶斯网络故障诊断方法

技术领域

本发明属于大规模非平稳过程故障诊断领域，特别是针对一种面向百万千瓦超超临界机组的非平稳性的多层因果结构的贝叶斯网络故障诊断方法。

背景技术

在高速发展的21世纪，火力发电行业的大规模化、复杂化的方向发展成为主流。随着科学技术的不断发展以及全球能源的急速消耗，火力发电行业对能效的要求越来越高。前期的百万千瓦超临界机组已经不足以满足要求，于是百万千瓦超超临界机组开始占据主导地位，该机组是世界上最先进的高效、大容量燃煤发电设备，具有明显的能效优势，是我国电力工业发展的代表性机组和主流方向。近年来，我国逐渐形成了以大容量、高参数、低能耗的超临界和超超临界机组为主体的电力能源结构，取代了以往高能耗小型火电机组为主导的陈旧模式。其中，我国百万千瓦超超临界机组保有量世界第一，并且发展空间巨大。然而，这种大型燃煤发电过程除了具有具有规模庞大、设备众多、参数多样化且相互影响等方面的特点，由于频繁调峰，过程扰动，设备老化等原因还会具有明显的非平稳特性。

综上所述，火力发电过程是一个生产工艺流程长、设备多样化、变量高度耦合、非平稳性强、安全要求高的工业过程。由此可知，该过程一旦发生故障将会带来人生安全和经济效益的双重损失。常见的典型故障有由于磨损、腐蚀、设备老化、应力拉伸等原因造成的四管泄露，由于堵塞、漏粉、断裂等原因造成的磨煤机故障，由于振动大、温度高等原因造成的送风机、引风机故障等。其中四管四管爆裂故障一旦发生其后果将十分严重。

确保工业过程正常运行、提高产品质量以及人员安全的重要技术被称为故障诊断。在检测到故障后，如何进一步判断故障类型，隔离故障变量，并及时定位故障根源变量，从源头根治故障是故障诊断领域所研究的主要内容。对于火力发电机组，由于其设备分布空间广、数量庞大，其发生的故障往往会随着时间的推移从故障根源变量传播至其他相关流程变量。而又由于火力发电机组存在非平稳性，对故障传播的诊断带来了一定的干扰因素。在实际生产过程中，常常由于非平稳特性的存在，无法正确地提取变量间的因果关系，导致故障发生后不能快速地确定故障传播的路径，浪费了检修时的有效时间，徒劳无功。

针对火力发电机组故障诊断中的上述问题，前人已经尝试从多种角度进行了探讨与研究，并相继提出了一些针对非平稳性过程的故障诊断方法。总体来说，主要包括基于解析模型的方法和基于数据驱动的方法。随着工业自动化程度的加深和信息化水平的不断提高，大量的过程数据得以存储并利用，因此，基于数据驱动的方法得到越来越多研究者的关注。然而，在数据驱动的方法中传统的用差分消除非平稳特性的方法，可能会损失有效的因果信息，导致建立的故障追溯网络中的因果关系不完善甚至错误，最终造成故障根源变量的误判。因此，百万千瓦超超临界机组具有明显的大规模、复杂化和非平稳特性，这给故障诊断带来了很大的挑战。

发明内容

本发明针对百万千瓦超超临界机组这一典型的大规模、非平稳过程提出了一种面向百万千瓦超超临界机组的非平稳性的多层因果结构的贝叶斯网络故障诊断方法。该方法充分提取百万千瓦超超临界机组的非平稳变量间的因果关系，在这种局部精细化的因果关系的上层，又构建了相对完善的全局性的整个系统的因果网络。随后，结合具有因果推理能力的贝叶斯网络模型，对得到的精细化因果关系网进行合理利用，能有效找到故障的传播途径，并准确定位故障根源变量。

本发明的目的通过以下技术方案实现：一种面向百万千瓦超超临界机组的非平稳性的多层因果结构的贝叶斯网络故障诊断方法，该方法包括以下步骤：

(1)获取待分析数据：设一个热力系统生产过程具有J个测量变量和操作变量。则针对正常流程，每一次采样可以得到一个1×J的向量，采样N次后得到的数据表述为一个二维矩阵X(N×J)，所述测量变量为运行过程中可被测量的状态参数，包括流量、温度、速率等；所述操作变量包括进风量、给料量、阀门开度等。建模过程选取一组正常数据，表示为X_i(N_i×J)，i表示第i组正常数据，测试过程选取某一种或几种故障案例，故障案例表示为F_case(N_f×J),N_f表示故障案例casef的采样总次数，其中casef∈{case1,...,casen}，n等于故障案例的总个数。

(2)识别非平稳变量：通过单位根检验方法(Augmented Dickey-Fuller， ADF)识别矩阵X_i(N_i×J)中的非平稳变量，得到非平稳变量数据矩阵V_ns(N_i×J_ns)， J_ns表示非平稳变量个数。

(3)针对非平稳变量建立局部的因果关系模型，该步骤通过以下子步骤实现：

(3.1)提取平稳特征：利用(2)中得到的非平稳变量数据矩阵V_ns(N_i×J_ns) 建立稀疏协整模型。非平稳变量数据矩阵可以表示为 其中v_T表示第T次采样的采样值。建立稀疏协整模型具体包括以下子步骤：

(3.1.1)对v_t建立向量自回归模型，t＝1,2,…J_ns

v_t＝Π₁v_t-1+…+Π_pv_t-p+c+μ_t (1)

其中，Π_i(J_ns×N_i)为系数矩阵，i＝1,2,…p，μ_t(N_i×1)为高斯白噪声，μ_t～N_i(0,Ξ)，Ξ表示方差；c(N_i×1)为常数向量，p为模型阶次；

(3.1.2)在公式(1)中两端减去v_t-1得到误差修正模型

其中， 为(J_ns×J_ns)的单位矩阵；

(3.1.3)将步骤(3.1.2)中的Γ分解为两个列满秩的矩阵Γ＝ΑΒ^T公式(2) 变为

其中Α(N_i×R)，Β(N_i×R)；

(3.1.4)通过极大似然估计方法对公式(3)中的协整向量矩阵Β进行估计

其中，L(*)表示极大似然函数，tr(*)表示矩阵的迹。 t∈{p+1,...,N_i}；Ω＝(Ω₁,...,Ω_p-1)^T，Θ＝Ξ^-1；

(3.1.5)对公式(4)的极大似然估计转化为特征方程求解过程

其中，参数矩阵Θ_i及Φ_i可以通过最小二乘算法求得；

(3.1.6)对公式(4)的目标函数加入惩罚函数得到稀疏协整向量

其中，P₁，P₂，P₃为参数B，Ω，Θ的惩罚函数，这里的惩罚函数采用1范数形式。调整参数λ₁和λ₂采用交叉检验来确定，调整参数λ₃采用贝叶斯信息准则确定。通过对公式(6)的求解可以得到稀疏的协整向量。

(3.2)协整向量分组：利用公式(6)中得到的稀疏协整向量可以根据平稳残差序列的平稳程度将变量进行划分。包括以下子步骤：

(3.2.1)利用ADF检验衡量每个稀疏协整向量的残差序列的平稳程度，并记录残差序列的ADF检验统计量t_i。

(3.2.2)将t_i进行升序排序，最小的检验统计量值对应的稀疏协整向量被保留。稀疏协整向量中的非零元素对应的变量被分到子组中，并记为V_{sf_i}。

(3.2.3)将V_{sf_i}中的变量从原始数据集V_ns中移除。

(3.2.4)重复迭代步骤(3.2.1)-(3.2.3)直到所有变量都被分到不同的子组中。

(3.3)建立局部的非平稳变量间的因果关系，具体子步骤包括如下：

(3.3.1)利用EViews统计软件中的格兰杰因果分析功能，分别处理每一个得到的子组中的非平稳变量，得到每组V_{sf_i}中的两两非平稳变量的F 统计量。m_i是子组中非平稳变量的个数；

(3.3.2)假设变量v_i不是v_j的格兰杰原因(i,j∈{1,...,m_i})，若该假设的F统计量小于0.05，则拒绝该假设，即v_i是v_j的格兰杰原因；如果F统计量大于等于 0.05，则保留原假设不变。

(3.3.3)根据步骤(3.3.2)计算得到的格兰杰原因关系，将每个子组V_{sf_i}中的非平稳变量以带箭头的直线连接，起点为格兰杰原因，箭头方向指向结果变量。

(3.4)建立全局的因果关系网络，具体子步骤如下：

(3.4.1)把每个子组V_{sf_i}对应的稀疏协整向量记为平稳特征。

(3.4.2)将平稳特征和平稳变量一起进行格兰杰因果分析。具体子步骤与步骤(3.3)中相同。

(3.5)构建基于因果关系结构的贝叶斯网络诊断模型，包括以下子步骤：

(3.5.1)贝叶斯网络的结构为上文所得到的因果关系结构，包括全局的因果关系和局部因果关系。

(3.5.2)计算条件概率表：对于数据X(N×J)，有N个观测样本和J个变量。x_j ⁽ⁱ⁾表示第j个观测样本的第i个变量的值。条件概率和全概率的计算公式如下：

其中，x∈{1,...,k}，k表示状态个数。

当x_j ⁽ⁱ⁾＝x时，(x_j ⁽ⁱ⁾＝x)＝1，否则(x_j ⁽ⁱ⁾＝x)＝0。

当x_j ⁽ⁱ⁾＝x，且y_j ⁽ⁱ⁾＝y时，(x_j ⁽ⁱ⁾＝x&y_j ⁽ⁱ⁾＝y)＝1，否则(x_j ⁽ⁱ⁾＝x&y_j ⁽ⁱ⁾＝y)＝0；

当y_j ⁽ⁱ⁾＝y时，(y_j ⁽ⁱ⁾＝y)＝1，否则(y_j ⁽ⁱ⁾＝y)＝0

(3.6)模型建好后，开始进入诊断策略，包括以下几个步骤：

(3.6.1)提取故障变量：通过稀疏FDFDA从不同方向上选出受到故障影响的故障变量，记作V_F{V_{F_1},...,V_{F_l}}，l表示故障变量的数量。

(3.6.2)从(3.6.1)中提取的故障变量V_F{V_{F_1},...,V_{F_l}}中选择故障根源变量追溯的输入证据：如果在V_F中，存在V_{F_i}，i∈{1,...,l}，是一个没有子节点的根节点，那么该变量作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。如果在V_F中，不存在根节点的变量，则选择具有子节点数最少的V_{F_i}，i∈{1,...,l}作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。

(3.6.3)当追溯到平稳特征时，如果非平稳变量不在V_F中，则停止追溯；如果在V_F中存在非平稳变量，则选择具有子节点数最少的非平稳变量作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。

(3.7)当经过步骤(3.6)后，诊断网络中每个变量节点的故障概率得到更新，得到一条完整的故障传播路径。从证据变量开始往回追溯，直到某个故障变量的故障概率小于50％时停止。故障传播路径的起点，即为故障根源变量。

本发明的有益效果在于：本发明针对百万千瓦超超临界机组这一典型的大规模、非平稳过程提出一种面向百万千瓦超超临界机组的非平稳性的多层因果结构的贝叶斯网络故障诊断方法。该方法考虑到火电机组的非平稳特性，首先将非平稳变量分离，充分提取百万千瓦超超临界机组中非平稳变量间的因果关系，避免了伪回归等因非平稳性带来的对因果关系的不利影响。在这种局部精细化的因果关系的上层，又通过构建平稳特征，建立了非平稳变量和平稳变量间的合理化因果联系，挖掘了相对完善的全局性的整个系统的因果网络。随后，结合具有因果推理能力的贝叶斯网络模型，对得到的精细化因果关系网进行合理利用，有效地找到故障的传播途径，并准确定位故障根源变量。该方法克服了大型燃煤机组因非平稳特性而导致的因果关系扭曲的问题，充分挖掘变量间所包含的潜在因果信息，建立了故障根源追溯的诊断准则，大大提高了复杂非平稳过程故障诊断的性能。有助于火电厂对故障设备进行高效准确的定位和检修，从而保证了百万千瓦超超临界发电机组的安全可靠运行。

附图说明：

图1是本发明面向百万千瓦超超临界机组非平稳性的多层因果结构的贝叶斯网络故障诊断方法流程图；

图2是本发明方法中非平稳变量间的局部因果关系图；

图3是本发明方法的多层因果结构的贝叶斯网络故障诊断方法的应用结果图。

具体实施方式

下面结合附图及具体实例，对本发明作进一步详细说明。

本发明以浙能集团下属嘉华电厂7号机组为例，该机组的功率为10000MW，为百万千瓦超超临界机组，包括众多过程设备。本次实例选取火电机组的1个典型磨煤机故障作为本次诊断过程详细说明的故障类型。

如图1所示，本发明是一种面向百万千瓦超超临界机组的非平稳性的多层因果结构的贝叶斯网络故障诊断方法，包括以下步骤：

(1)获取待分析数据：设一个百万千瓦超超临界机组的磨煤机设备生产过程具有J个测量变量和操作变量。则针对正常流程，每一次采样可以得到一个 1×J的向量，采样N次后得到的数据表述为一个二维矩阵X(N×J)。建模过程选取一组正常数据，表示为X_i(N_i×J)，i表示第i组正常数据，测试过程选取磨煤机F的一个典型故障案例，表示为case 1。本实例中，采样周期为1分钟，该故障类型采集1000个样本，12个过程变量，所测变量为运行过程中的润滑油温度、分离器压力、密封风压差等，详见表1。

表1：磨煤机的参数表

(2)识别非平稳变量：通过单位根检验方法(Augmented Dickey-Fuller， ADF)识别矩阵X_i(N_i×J)中的非平稳变量，得到由JHMEAS2，JHMEAS9， JHMEAS10这3个非平稳变量组成的1个非平稳变量数据矩阵V_ns(N×J_ns)，J_ns表示非平稳变量个数，此处J_ns＝3。

(3)针对非平稳变量建立局部的因果关系模型，该步骤通过以下子步骤实现：

(3.1)提取平稳特征：利用(2)中得到的非平稳变量数据矩阵V_ns(N_i×J_ns) 建立稀疏协整模型。非平稳变量数据矩阵可以表示为 V_ns(N×J_ns)＝[JHMEAS_2,JHMEAS_9,JHMEAS_10]。建立稀疏协整模型具体包括以下子步骤：

(3.1.1)对v_t建立向量自回归模型，t＝1,2,…J_ns

v_t＝Π₁v_t-1+…+Π_pv_t-p+c+μ_t (1)

其中，Π_i(J_ns×N_i)为系数矩阵，i＝1,2,…p，μ_t(N_i×1)为高斯白噪声，μ_t～N_i(0,Ξ)，Ξ表示方差；c(N_i×1)为常数向量，p为模型阶次；

(3.1.2)在公式(1)中两端减去v_t-1得到误差修正模型

其中， 为(J_ns×J_ns)的单位矩阵；

(3.1.3)将步骤(3.1.2)中的Γ分解为两个列满秩的矩阵Γ＝ΑΒ^T公式(2) 变为

其中Α(N_i×R)，Β(N_i×R)；

(3.1.4)通过极大似然估计方法对公式(3)中的协整向量矩阵Β进行估计

其中，L(*)表示极大似然函数，tr(*)表示矩阵的迹。 t∈{p+1,...,N_i}；Ω＝(Ω₁,...,Ω_p-1)^T，Θ＝Ξ^-1；

(3.1.5)对公式(4)的极大似然估计转化为特征方程求解过程

其中，参数矩阵Θ_i及Φ_i可以通过最小二乘算法求得；

(3.1.6)对公式(4)的目标函数加入惩罚函数得到稀疏协整向量

其中，P₁，P₂，P₃为参数B，Ω，Θ的惩罚函数，这里的惩罚函数采用1范数形式。调整参数λ₁和λ₂采用交叉检验来确定，调整参数λ₃采用贝叶斯信息准则确定。通过对公式(6)的求解可以得到稀疏的协整向量。

(3.2)协整向量分组：利用公式(6)中得到的稀疏协整向量可以根据平稳残差序列的平稳程度将变量进行划分。包括以下子步骤：

(3.2.1)利用ADF检验衡量每个稀疏协整向量的残差序列的平稳程度，并记录残差序列的ADF检验统计量t_i。

(3.2.2)将t_i进行升序排序，最小的检验统计量值对应的稀疏协整向量被保留。稀疏协整向量中的非零元素对应的变量被分到子组中，并记为V_{sf_i}。

(3.2.3)将V_{sf_i}中的变量从原始数据集V_ns中移除。

(3.2.4)重复迭代步骤(3.2.1)-(3.2.3)直到所有变量都被分到不同的子组中。

(3.3)建立局部的非平稳变量间的因果关系，具体子步骤包括如下：

(3.3.1)利用EViews统计软件中的格兰杰因果分析功能，分别处理每一个得到的子组中的非平稳变量，得到每组V_{sf_i}中的两两非平稳变量的F 统计量。m_i是子组中非平稳变量的个数；

(3.3.2)假设变量v_i不是v_j的格兰杰原因(i,j∈{1,...,m_i})，若该假设的F统计量小于0.05，则拒绝该假设，即v_i是v_j的格兰杰原因；如果F统计量大于等于 0.05，则保留原假设不变。

(3.3.3)根据步骤(3.3.2)计算得到的格兰杰原因关系，将每个子组V_{sf_i}中的非平稳变量以带箭头的直线连接，起点为格兰杰原因，箭头方向指向结果变量。以非平稳变量间的因果关系求算为例，见表2。

表2：变量间因果关系的构建过程，以非平稳变量为例

(3.4)建立全局的因果关系网络，具体子步骤如下：

(3.4.1)把每个子组V_{sf_i}对应的稀疏协整向量记为平稳特征。

(3.4.2)将平稳特征和平稳变量一起进行格兰杰因果分析。具体子步骤与步骤(3.3)中相同。

(3.5)构建基于因果关系结构的贝叶斯网络诊断模型，包括以下子步骤：

(3.5.1)贝叶斯网络的结构为上文所得到的因果关系结构，包括全局的因果关系和局部因果关系。

(3.5.2)计算条件概率表：对于数据X(N×J)，有N个观测样本和J个变量。x_j ⁽ⁱ⁾表示第j个观测样本的第i个变量的值。条件概率和全概率的计算公式如下：

其中，x∈{1,...,k}，k表示状态个数。

当x_j ⁽ⁱ⁾＝x时，(x_j ⁽ⁱ⁾＝x)＝1，否则(x_j ⁽ⁱ⁾＝x)＝0。

当x_j ⁽ⁱ⁾＝x，且y_j ⁽ⁱ⁾＝y时，(x_j ⁽ⁱ⁾＝x&y_j ⁽ⁱ⁾＝y)＝1，否则(x_j ⁽ⁱ⁾＝x&y_j ⁽ⁱ⁾＝y)＝0；

当y_j ⁽ⁱ⁾＝y时，(y_j ⁽ⁱ⁾＝y)＝1，否则(y_j ⁽ⁱ⁾＝y)＝0

(3.6)模型建好后，开始进入诊断策略，包括以下几个步骤：

(3.6.1)提取故障变量：通过稀疏FDFDA从不同方向上选出磨煤机F受到故障影响的故障变量，记作记作V_F{JHMEAS_12,JHMEAS_9,JHMEAS_11}。

(3.6.2)从(3.6.1)中提取的故障变量V_F{V_{F_1},...,V_{F_l}}中选择故障根源变量追溯的输入证据：如果在V_F中，存在V_{F_i}，i∈{1,...,l}，是一个没有子节点的根节点，那么该变量作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。如果在V_F中，不存在根节点的变量，则选择具有子节点数最少的V_{F_i}，i∈{1,...,l}作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。在V_F{JHMEAS_12,JHMEAS_9,JHMEAS_11} 中，选择JHMEAS_12作为证据输入。

(3.6.3)当追溯到平稳特征时，如果非平稳变量不在V_F中，则停止追溯；如果在V_F中存在非平稳变量，则选择具有子节点数最少的非平稳变量作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。在V_F{JHMEAS_12,JHMEAS_9,JHMEAS_11}中，选择JHMEAS_9作为非平稳变量的证据输入。

(3.7)当经过步骤(3.6)后，诊断网络中每个变量节点的故障概率得到更新，得到一条完整的故障传播路径。从证据变量开始往回追溯，直到某个故障变量的故障概率小于50％时停止。故障传播路径的起点，即为故障根源变量。本次案例的诊断结果如图3所示。

本发明中选择的实际案例是嘉华机组磨煤机F由于密封风压差较低导致的故障报警。从图3中可以看出，从证据变量磨煤机F液压油压力(JHMEAS_12) 出发，根据故障概率的大小，可以准确追溯到故障根源变量磨煤机F密封风/炉膛差压(JHMEAS_6)。同时，还能通过局部精细化的因果关系，追溯到非平稳的故障变量磨煤机F润滑油温度(JHMEAS_9)和磨煤机F润滑油压力 (JHMEAS_10)。基于本发明方法的故障诊断方法提高了实际在线故障诊断的准确性和可靠性。本发明方法的故障诊断的优越性体现在通过多层的因果关系构建利用贝叶斯网络的因果推断能力，不仅能分析平稳变量间的因果追溯路径，还能定位非平稳变量的故障根源变量。以此发明方法作为参考，可以帮助现场工程师做准确有效的诊断。此外，在准确诊断后可以集中人力赶赴现场进行有效及时的排查，在早期发现不良情况，避免严重事故的发生，由此保证了实际生产过程的安全性和可靠性。

Claims (1)

1.一种面向百万千瓦超超临界机组非平稳特性的基于因果分析的多层贝叶斯网络故障诊断方法，其特征在于，该方法包括以下步骤:

(1)获取待分析数据：设一个热力系统生产过程具有J个测量变量和操作变量。则针对正常流程，每一次采样可以得到一个1×J的向量，采样N次后得到的数据表述为一个二维矩阵X(N×J)，所述测量变量为运行过程中可被测量的状态参数，包括流量、温度、速率等；所述操作变量包括进风量、给料量、阀门开度等。建模过程选取一组正常数据，表示为X_i(N_i×J)，i表示第i组正常数据，测试过程选取某一种或几种故障案例，故障案例表示为F_case(N_f×J),N_f表示故障案例casef的采样总次数，其中casef∈{case1,...,casen}，n等于故障案例的总个数。

(2)识别非平稳变量：通过单位根检验方法(Augmented Dickey-Fuller，ADF)识别矩阵X_i(N_i×J)中的非平稳变量，得到非平稳变量数据矩阵V_ns(N_i×J_ns)，J_ns表示非平稳变量个数。

(3)针对非平稳变量建立局部的因果关系模型，该步骤通过以下子步骤实现：

(3.1)提取平稳特征：利用(2)中得到的非平稳变量数据矩阵V_ns(N_i×J_ns)建立稀疏协整模型。非平稳变量数据矩阵可以表示为 其中v_T表示第T次采样的采样值。建立稀疏协整模型具体包括以下子步骤：

(3.1.1)对v_t建立向量自回归模型，t＝1,2,…J_ns

v_t＝Π₁v_t-1+…+Π_pv_t-p+c+μ_t (1)

其中，Π_i(J_ns×N_i)为系数矩阵，i＝1,2,…p，μ_t(N_i×1)为高斯白噪声，μ_t～N_i(0,Ξ)，Ξ表示方差；c(N_i×1)为常数向量，p为模型阶次；

(3.1.2)在公式(1)中两端减去v_t-1得到误差修正模型

其中， 为(J_ns×J_ns)的单位矩阵；i＝1,2,...p-1。

(3.1.3)将步骤(3.1.2)中的Γ分解为两个列满秩的矩阵Γ＝ΑΒ^T公式(2)变为

其中Α(N_i×R)，Β(N_i×R)；

(3.1.4)通过极大似然估计方法对公式(3)中的协整向量矩阵Β进行估计

其中，L(*)表示极大似然函数，tr(*)表示矩阵的迹。 t∈{p+1,...,N_i}；Ω＝(Ω₁,...,Ω_p-1)^T，Θ＝Ξ^-1；

(3.1.5)对公式(4)的极大似然估计转化为特征方程求解过程

其中，参数矩阵Θ_i及Φ_i可以通过最小二乘算法求得；

(3.1.6)对公式(4)的目标函数加入惩罚函数得到稀疏协整向量

其中，P₁，P₂，P₃为参数B，Ω，Θ的惩罚函数，这里的惩罚函数采用1范数形式。调整参数λ₁和λ₂采用交叉检验来确定，调整参数λ₃采用贝叶斯信息准则确定。通过对公式(6)的求解可以得到稀疏的协整向量。

(3.2)协整向量分组：利用公式(6)中得到的稀疏协整向量可以根据平稳残差序列的平稳程度将变量进行划分，包括以下子步骤：

(3.2.1)利用ADF检验衡量每个稀疏协整向量的残差序列的平稳程度，并记录残差序列的ADF检验统计量t_i。

(3.2.2)将t_i进行升序排序，最小的检验统计量值对应的稀疏协整向量被保留。稀疏协整向量中的非零元素对应的变量被分到子组中，并记为V_{sf_i}。

(3.2.3)将V_{sf_i}中的变量从原始数据集V_ns中移除。

(3.2.4)重复迭代步骤(3.2.1)-(3.2.3)直到所有变量都被分到不同的子组中。

(3.3)建立局部的非平稳变量间的因果关系，具体子步骤包括如下：

(3.3.1)利用EViews统计软件中的格兰杰因果分析功能，分别处理每一个得到的子组中的非平稳变量，得到每组V_{sf_i}中的两两非平稳变量的F统计量。m_i是子组中非平稳变量的个数；

(3.3.2)假设变量v_i不是v_j的格兰杰原因(i,j∈{1,...,m_i})，若该假设的F统计量小于0.05，则拒绝该假设，即v_i是v_j的格兰杰原因；如果F统计量大于等于0.05，则保留原假设不变。

(3.3.3)根据步骤(3.3.2)计算得到的格兰杰原因关系，将每个子组V_{sf_i}中的非平稳变量以带箭头的直线连接，起点为格兰杰原因，箭头方向指向结果变量。

(3.4)建立全局的因果关系网络，具体子步骤如下：

(3.4.1)把每个子组V_{sf_i}对应的稀疏协整向量记为平稳特征。

(3.4.2)将平稳特征和平稳变量一起进行格兰杰因果分析。具体子步骤与步骤(3.3)中相同。

(3.5)构建基于因果关系结构的贝叶斯网络诊断模型，包括以下子步骤：

(3.5.1)贝叶斯网络的结构为上文所得到的因果关系结构，包括全局的因果关系和局部因果关系。

(3.5.2)计算条件概率表：对于数据X(N×J)，有N个观测样本和J个变量。x_j ⁽ⁱ⁾表示第j个观测样本的第i个变量的值。条件概率和全概率的计算公式如下：

其中，x∈{1,...,k}，k表示状态个数。

当x_j ⁽ⁱ⁾＝x时，(x_j ⁽ⁱ⁾＝x)＝1，否则(x_j ⁽ⁱ⁾＝x)＝0。

当x_j ⁽ⁱ⁾＝x，且y_j ⁽ⁱ⁾＝y时，(x_j ⁽ⁱ⁾＝x&y_j ⁽ⁱ⁾＝y)＝1，否则(x_j ⁽ⁱ⁾＝x&y_j ⁽ⁱ⁾＝y)＝0；

当y_j ⁽ⁱ⁾＝y时，(y_j ⁽ⁱ⁾＝y)＝1，否则(y_j ⁽ⁱ⁾＝y)＝0

(3.6)模型建好后，开始进入诊断策略，包括以下几个步骤：

(3.6.1)提取故障变量：通过稀疏FDFDA从不同方向上选出受到故障影响的故障变量，记作V_F{V_{F_1},...,V_{F_l}}，l表示故障变量的数量。

(3.6.2)从(3.6.1)中提取的故障变量V_F{V_{F_1},...,V_{F_l}}中选择故障根源变量追溯的输入证据：如果在V_F中，存在V_{F_i}，i∈{1,...,l}，是一个没有子节点的根节点，那么该变量作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。如果在V_F中，不存在根节点的变量，则选择具有子节点数最少的V_{F_i}，i∈{1,...,l}作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。

(3.6.3)当追溯到平稳特征时，如果非平稳变量不在V_F中，则停止追溯；如果在V_F中存在非平稳变量，则选择具有子节点数最少的非平稳变量作为输入证据输入到模型中，证据变量的故障概率被置为100％，更新网络的条件概率表。

(3.7)当经过步骤(3.6)后，诊断网络中每个变量节点的故障概率得到更新，得到一条完整的故障传播路径。从证据变量开始往回追溯，直到某个故障变量的故障概率小于50％时停止。故障传播路径的起点，即为故障根源变量。