岳健国 郭瑞 傅行军 田新启 宗曜犇 王旭

（1.大型发电装备安全运行与智能测控国家工程研究中心；2.东南大学能源与环境学院）

0 引言

目前，我国正朝着实现“双碳”目标大步迈进，构建“以新能源为主体的新型电力系统”是实现“双碳”目标的重要手段。新型电力系统中高比例可再生能源和高比例电力电子设备接入电力系统的特征会导致电网运行波动性升高，从而对火电机组深度调峰能力提出了更高的要求。同时，火电机组负荷的多变使得各种辅机的运行状况发生较大变化，深入研究辅机参数间内在规律对维持其稳定健康运行及保证发电机组运行效率有重要意义。

当前发电设备运行状态分析方法主要分为两类，即基于模型的方法和基于数据驱动的方法。其中，基于数据驱动的方法在信息技术和电力设备运行数据采集技术的大力推动下得到快速发展，此类方法能揭示设备运行参数间的复杂耦合关系，并能对设备运行状态和故障进行有效监测及预警。

目前数据分析技术主要依靠相关性概念分析设备运行数据间的关联关系，在设备强相关关系的基础上，若能结合参数间因果不对称性，利用因果推断技术研究设备运行变量间的因果关系，将能更加深入挖掘设备内在规律性，为机组辅机安全稳定运行提供帮助，并能为现场运行调整提供有力指导。

因果推断理论已在医学、生物制药、社会实验、气候变化等问题的研究上已获得广泛应用[1-3]。同时，已有部分将因果推断技术应用到工业领域中去的尝试，穆钢[4]在信息几何因果推理方法的基础上提出一种逆信息熵因果推理的方法用于电力系统运行数据中因果关系分析，Zhang[5]将个体平均因果效应ITE用到了室内空调系统的故障诊断中，通过正常异常情况下数据分类对比结合因果关系判断故障类型。Sun[6]提出了一种将因果图和滑动窗口KPC结合的方法用于复杂工业过程的故障监测中。

但现有研究中将因果发现模型用于实际设备中的实践较少，其中因果推断技术中不同因果发现模型所要求的不同数据生成方式和噪声假设是其用于实际设备时进行模型选用及准确性和有效性判断的一大障碍。

本文对因果关系发现及评估算法进行了改进，并将专家知识引入数据驱动的因果推断技术中。同时，以某一次风机为研究对象验证了文中所提方法的有效性，利用本方法构建了其合理可信的因果关系网络、并通过计算因果效应强度反应了其状态恶化过程和状态突变节点。

1 基于因果函数的改进因果推断方法

1.1 改进因果关系发现及评估方法

因果关系发现是因果推断技术的一部分，目前常见的因果发现算法主要分为两类，分别是基于约束的方法和基于因果函数模型的方法，基于约束的方法由于无法解决Markov等难题所以存在部分无法判断的因果关系方向[7]。而基于因果函数模型的方法通常对数据的生成方式和噪声有较强的假设性。计算机领域往往有标准集对模型因果关系识别结果有效性和准确度进行验证，但实际生产设备中不同参数间因果关系通常只能依靠物理模型或随机试验进行验证，而物理模型的复杂性和随机试验的高成本或不可实施性使得模型对于实际设备因果关系识别结果的验证往往无法顺利进行。

因此，为了提高从观测数据识别因果关系的准确度，本文提出一种改进因果关系发现及评估方法，首先对原始数据进行分组，在分组大小保证所选模型可正常识别的情况下，分别使用所选模型对不同组的数据进行因果关系发现。据此，使用因果保持率CRR（Causality retention rate）描述同一因果方向在多组数据中被识别出的概率:

式中，CRRi代表第i条因果关系方向的保持率；pi代表第i条因果关系出现的次数；N代表组的个数。

为了校验所选模型准确度，本文提出设定强因果关系和非因果关系。结合专家知识将已知具有明确物理意义的参数间因果方向设置为强制因果方向，已知绝对不是因果方向的路径设置为非因果方向。据此判断每组因果图识别结果的置信度Si,并使用置信度修改每组数据得到的因果图中除强因果方向和非因果方向外其他方向在总体因果图的邻接矩阵中的权重。因果图置信度Si计算公式如下:

式中，Si为每组识别得到的因果图的置信度；m为分组的总数量；hi为错误或未识别出的强因果方向的数量；li为错误识别的非因果方向的数量；M为强因果方向和非因果方向的总数量。

之后根据Si对因果关系保持率进行修改得:

除此之外，本文提出利用CRRi计算因果信息不对称度DDAi（Degree of directional asymmetry），用以描述同一因果关系两方向之间差异，判断两参数间因果方向识别结果的可靠性。

式中，CRRi+和CRRi-代表同一关系两个方向的因果关系保持率。

基于因果函数的模型识别准确率都较为依赖不同模型所提出的数据生成方式假设和两参数间的因果信息不对称度。其中，数据生成方式假设可以依靠专家经验进行辨别，而因果信息不对称度DDAi可以通过式（4）进行计算。对CRRi和DDAi设置一定的阈值，即可最终确定基本正确的因果关系网络。完整改进因果关系发现及评估方法流程图如图1所示。

图1 完整改进因果关系发现及评估方法流程图Fig.1 The flow diagram of improved causal relationship discovery and evaluation method

1.2 线性非高斯无环模型

因果关系发现作为因果推断技术的一部分，针对基于约束的方法存在的马尔可夫等价类问题，许多学者从因果作用机制出发，以结构因果模型为基础结合数据分布特性提出了不同的因果函数模型。

线性非高斯无环模型LiNGAM[8]（Linear Non-Gaussian Acyclic Model）属于基于因果函数的因果发现模型，此模型对数据生成方式提出了三点假设:

a）数据生成方式是线性的，变量xi是对应的原因变量xj,噪声变量ei和常数ci的线性求和，即:

b）模型中不存在未观测到的混杂因子；

c）噪声变量ei独立非高斯，即ei服从非零方差的非高斯分布。

通过数据预处理减去每个变量的均值，消去公式（5）中的常数项ci，从而得到以下矩阵方程:

式中，X为i维随机变量；B为i×i的连接矩阵；N为i维的非高斯随机噪声变量。

对式（6）求解X可得:

式中，A=(I-B)-1，A矩阵可以通过适当的置换算法置换为下三角矩阵或近似为下三角矩阵。

式（7）和噪声变量ei的独立非高斯性构成了标准的线性非高斯无环模型。通过独立成分分析ICA算法[9]（Independent Component Analysis）等算法可以对LiNGAM 模型进行求解,获得包含具有所寻求的连接强度bij的矩阵B，从而即可得到具体因果方向，最后通过修剪边枝算法[10]得到最终因果图网络。

1.3 因果效应强度计算

在得到各参数间的因果关系网络后，因果推断技术可以量化各因果方向“因”对“果”的影响程度即因果效应强度。对原因X按是否进行了干预分为x=1和x=0组，Y1(x)和Y0(x)分别代表单个个体接受与未接受x干预后的结果。虽然无法同时观察到单个个体的Y1(x)和Y0(x)从而计算x的作用结果，但是可以识别总体的因果效应，总体的因果效应强度为进行干预与未进行干预时结果差值的期望ATE(average treatment effect)，即:

当数据维度升高时，在两组间寻找除X外其他参数均相同的个体几乎是不可能的，因此本文选用倾向得分[11]作为多维数据匹配依据。

倾向得分匹配的原理是，通过以所计算因果关系中的“因”为因变量，以协变量作为自变量，利用Logistic回归建立回归模型估计每个样本接受实验的可能性作为倾向性得分，然后，根据倾向得分使用近邻匹配法在实验组和对照组间进行匹配。由此，就可以计算具有多维可观测特征样本中实验组收到参数X的干预后的因果效应ATT(x)(average treatment effect for thextreated group)，计算公式为式（9）。

最后，对式（9）所有ATT(x)进行加权平均得到ATT式（10）:

2 一次风机样本实验验证

2.1 一次风机因果发现模型选择

本文选取某火电厂一次风机整年运行数据为研究对象，测点包括调节挡板阀位、风机电流、轴承温度（前、中、后）、轴承X和Y方向振动信号，根据一次风机物理性质结合专家知识，本文设置了一条强因果关系和六条非因果关系，如图2所示。

图2 强因果关系与非因果关系示意图Fig.2 The diagram of strong causality and non-causality

从原始数据抽取1000 组样本，每组样本量10 万个数据点，分别使用不同因果发现模型进行识别，由公式（1.3）计算置信度，最终计算结果如表1所示：

表1 不同模型对此一次风机因果图识别结果置信度Tab.1 The Confidence coefficient of cause-and-effect diagram recognition result of primary air fan with different model

由表1 可见，对此一次风机，线性非高斯无环模型LiNGAM[8]的识别结果具有极高置信度，因此本文选用LiNGAM模型对一次风机进行因果关系发现。

2.2 因果函数模型可识别性分析

基于因果函数模型使用前需要验证参数不对称性和模型假设。一次风机的参数间固然具有较强的因果信息不对称性，电流增加会导致挡板阀门开度增大，但是反过来不成立。挡板阀门开度增大会对前轴承温度产生影响，虽然不能确定正向因果关系是否成立，但是反方向一定不成立，不对称性显著。

对于LiNGAM模型的三条假设:1）在一次风机稳定运行阶段，各参数间相互作用近似处于线性范围内，运行数据基本服从线性变化。2）本文选取数据中包含一次风机主要参数测点，不存在重要的未观测的混杂因子。3）假设风机运行过程中的噪声项为非高斯分布。

由此可见，一次风机中参数间不对称性清晰，参数情况满足LiNGAM模型假设。因此，对一次风机运行变量使用LiNGAM模型进行因果关系分析是可行的。

2.3 构建一次风机参数间因果关系网络

对此一次风机运行数据重新抽取1000 组样本，每组样本量10万个数据点，经过LiNGAM模型识别，计算每两个变量间的因果关系保持率CRR。最终结果因果方向保持率CRRi如表2所示。

表2 因果保持率CRRi 矩阵Tab.2 The causal retention rate CRRi matrix

由表2 可以看出，电流指向阀门开度的因果方向CRR为100%，对于六条非因果方向，除了后轴承温度对挡板阀门开度出现两次误判，其他五条在所有1000 次计算中均未被识别，说明此模型成功识别出了事前设定的一组强因果关系和六组非因果关系，表明LiNGAM模型可以识别出一次风机中基于物理规律的总体框架，并能够对其他复杂耦合方向进行识别。

因果关系网络中某一关系具体因果方向的识别结果可靠性主要依赖于各方向间参数的因果不对称度，在得到各方向间的因果方向保持率之后，通过计算其不对称度来确定识别结果的可信度。不对称度越高，因果方向识别结果可信度越高。设置CRRi的阈值为20%，对超过此阈值的因果关系方向计算因果信息不对称度DDAi。因果信息不对称度的计算结果如表3 所示，不对称度DDAi的置信值选取为0.2。之后，对于因果信息不对称度DDAi≥0.2 的边，选取CRRi较大的方向作为因果图中的因果方向，最终得到的因果关系网络如图3所示。

表3 因果关系不对称度Tab.3 The causality asymmetry

因果关系网络可以为更深层次地认识系统物理过程提供有利帮助。如图3所示，对于温度变化和振动的耦合问题，本系统中前中后轴承的温度变化都是轴承X和Y方向振动改变的“因”。这些因果方向的CRR基本均在50%以上且后轴承温度对轴承Y向振动的CRR达到了77.15%，说明此一次风机运行过程中温度变化是振动情况改变的一个原因，为风机的温度振动耦合作用的解耦指出具体的方向。

图3 通过LiNGAM模型识别得到的一次风机因果关系网络Fig.3 The causality network of primary air fan identified through LiNGAM model

2.4 因果效应强度计算

得到因果关系网络后，“因”对“果”的影响程度即因果效应强度也是因果推断的关注点之一。本文使用滑动窗口结合倾向得分匹配的方法量化不同因果关系方向之间的因果效应强度随时间的变化情况。本文选取窗口长度为24小时，滑动距离一小时。

通过对具有因果关系两参数的值按升高和降低以及不变进行分组处理后计算倾向得分，再根据倾向得分进行匹配。计算参数变化组即实验组的平均干预效应即为此因果方向的因果关系强度。

为进一步验证因果关系网络的准确性，本文通过三组反驳实验测试2.3 节得到的因果关系网络的准确性，分别是随机添加混杂因子，数据子集验证和安慰剂干预。虽然无法通过观测数据对因果关系网络进行证真，但可以通过反驳测试中模型出现的异常行为对其证伪。若模型结构正确，则随即添加混杂因子和数据子集测试时的因果效应强度与原始结果应近似相同，安慰剂干预测试是用独立的随机变量替换真正的原因变量计算因果效应，若模型正确，其结果应该接近0,三组测试结果如图4所示。

由三组反驳测试结果图4可以看出，随机添加混杂因子和数据子集测试时的因果效应强度均与原始效应差别很小，两者残差分别在0.1和0.12以内，说明在此因果图中随机添加混杂因子对其无较大影响，且用随机选择的子集替换给定的数据集时也具有较高的准确性。安慰剂干预测试的结果表明作为原因的独立随机变量对结果的因果效应强度都为0。在三种测试情况下此因果关系网络均没有异常产生，进一步证明了此因果关系网络拥有。

图4 反驳测试结果Fig.4 The refuted test results

2.5 因果状态监测

根据2.3 节得到的因果图和1.3 节所述的因果效应强度计算方法，计算得各因果方向间的主要因果效应变化曲线如图5所示。

图5 各因果方向因果效应强度随时间窗口变化曲线Fig.5 The causal effect intensity of each causal direction varies with the time window

由图5可以看出：（1）如图5-1所示，预设的强因果关系的效应变化情况符合实际物理模型。此一次风机全年电流-挡板阀位开度的因果效应强度基本一直维持在0.8附近，说明电流大小的变化将在很大程度上影响挡板开度的变化，符合实际物理模型。1月100点处的突降对应为停机，但在7月上旬750点处，此值由0.8突降至0附近，观察实际数据发现此时负荷增大，挡板阀门开度开至最大，电流对其作用效果达到饱和，进一步增大电流也不会引起阀门开度的变化，因此因果效应降为0。

（2）如图5-2 所示，电流大小对轴承X向振动的因果效应变化情况大概经历了200，750，1600三个窗口节点，其中，在窗口200 点，系统在长时间停机检修后投产，此效应缓慢爬升后稳定在0.5。在750点窗口即7月2 日附近，此值发生突变由0.5 升至0.7 高位附近保持，说明此时电流大小变化对轴承X向振动的影响增大。而在此处，轴承X向振动对轴承Y向振动的因果效应强度如图5-4 也发生了突变，由原来在0 波动跃升到0.4附近，说明此时两个方向的振动不再独立，上述讨论说明7 月2 日的高负荷导致系统振动情况发生突变性恶化。此后，1600点窗口对应11月初的停机检修，从原始数据及图5可以看出，此次停机检修后系统振动情况发生了较大改善。翻阅检修报告知，此风机失效轴承润滑油脂过少，在7 月中旬后逐渐出现轴承振动异常，并于11月1日进行停机检修替换失效轴承，系统状态变化过程与本文结果相符。

此例证明本方法成功跟踪了系统状态的变化，并发现原始数据无法得到的7 月2 日高负荷运行导致的系统状态突变情况。

3 结论

1）提出了一种改进因果关系发现和评估方法，通过引入因果关系保持率CRR、因果信息不对称度DDA，并结合专家知识设置强因果关系和非因果关系，可以对不同的因果函数模型因果关系识别准确度给出具体的评估结果，优化因果发现过程，提高因果关系网络的可靠性。

2）以某一次风机为对象进行因果关系发现，证明了改进的因果发现评估方法能成功判断不同模型的有效性和准确度并对因果关系发现过程进行修正。CRR评分表明LiNGAM模型对此一次风机强、非因果关系的识别准确率达到了99.95%。三组反驳测试结果均表明了此方法得到的因果关系网络的正确性。此方法为对发电设备进行因果关系发现提供了新的模型评估方法，解决了因果模型难以选择和准确度无法判断的问题。

3）因果状态监测提供了一种发电机组辅机运行状态分析的新角度，本文以某一次风机为例证明其不仅能够成功跟踪反映系统状态，还能够发现原始数据无法发现的系统状态恶化突变节点从而进行精确高效的故障告警。结合因果推断技术的数据分析方法可以成为研究系统更深层次本质规律和进行运行状态分析的有效手段。