CN113656906B - 一种面向燃气轮机的非平稳多变量因果关系分析方法 - Google Patents
一种面向燃气轮机的非平稳多变量因果关系分析方法 Download PDFInfo
- Publication number
- CN113656906B CN113656906B CN202110863191.4A CN202110863191A CN113656906B CN 113656906 B CN113656906 B CN 113656906B CN 202110863191 A CN202110863191 A CN 202110863191A CN 113656906 B CN113656906 B CN 113656906B
- Authority
- CN
- China
- Prior art keywords
- causal
- variables
- variable
- stationary
- gas turbine
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000001364 causal effect Effects 0.000 title claims abstract description 156
- 238000004458 analytical method Methods 0.000 title claims abstract description 42
- 238000000034 method Methods 0.000 claims abstract description 68
- 238000012546 transfer Methods 0.000 claims abstract description 52
- 239000013598 vector Substances 0.000 claims abstract description 28
- 238000012360 testing method Methods 0.000 claims abstract description 21
- 230000008569 process Effects 0.000 claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims description 37
- 238000009826 distribution Methods 0.000 claims description 20
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000002485 combustion reaction Methods 0.000 claims description 11
- 239000000446 fuel Substances 0.000 claims description 10
- 230000001052 transient effect Effects 0.000 claims description 5
- 230000001174 ascending effect Effects 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 238000012216 screening Methods 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 2
- 238000012545 processing Methods 0.000 claims description 2
- 230000007246 mechanism Effects 0.000 abstract description 9
- 230000005654 stationary process Effects 0.000 abstract description 4
- 238000011156 evaluation Methods 0.000 abstract 1
- 230000002452 interceptive effect Effects 0.000 abstract 1
- 239000007789 gas Substances 0.000 description 58
- 230000000694 effects Effects 0.000 description 11
- 238000010248 power generation Methods 0.000 description 6
- 238000004519 manufacturing process Methods 0.000 description 5
- 230000009286 beneficial effect Effects 0.000 description 4
- 238000012952 Resampling Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 3
- 238000012047 cause and effect analysis Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000010276 construction Methods 0.000 description 3
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000012423 maintenance Methods 0.000 description 3
- 230000001105 regulatory effect Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 239000003245 coal Substances 0.000 description 2
- 238000010924 continuous production Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000007710 freezing Methods 0.000 description 2
- 238000012353 t test Methods 0.000 description 2
- 239000002912 waste gas Substances 0.000 description 2
- LHASLBSEALHFGO-ASZAQJJISA-N 1-[(4s,5r)-4-hydroxy-5-(hydroxymethyl)oxolan-2-yl]-5-[[(2r,3r,4s,5s,6r)-3,4,5-trihydroxy-6-(hydroxymethyl)oxan-2-yl]oxymethyl]pyrimidine-2,4-dione Chemical compound C1[C@H](O)[C@@H](CO)OC1N1C(=O)NC(=O)C(CO[C@H]2[C@@H]([C@@H](O)[C@H](O)[C@@H](CO)O2)O)=C1 LHASLBSEALHFGO-ASZAQJJISA-N 0.000 description 1
- UGFAIRIUMAVXCW-UHFFFAOYSA-N Carbon monoxide Chemical compound [O+]#[C-] UGFAIRIUMAVXCW-UHFFFAOYSA-N 0.000 description 1
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 206010028980 Neoplasm Diseases 0.000 description 1
- 238000000692 Student's t-test Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 201000011510 cancer Diseases 0.000 description 1
- 239000000919 ceramic Substances 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000001311 chemical methods and process Methods 0.000 description 1
- 238000004140 cleaning Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000013329 compounding Methods 0.000 description 1
- 239000008358 core component Substances 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 239000003546 flue gas Substances 0.000 description 1
- 210000004907 gland Anatomy 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000013439 planning Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002195 synergetic effect Effects 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
- 238000012384 transportation and delivery Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 239000002918 waste heat Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种面向燃气轮机的非平稳多变量因果关系分析方法。该方法针对在燃气轮机运行过程中持续出现的大范围非平稳运行工况，通过秩向量抽取数据波动趋势转换为有限的符号序列集合，利用多变量符号转移熵对变量间信息传递的方向与强度进行量化评估。同时设计统计假设检验流程，通过对变量时序信息随机重组作为参照组，判断实际变量因果关系的显著性，实现了对变量间冗杂的虚假因果关系的简化。最后建立多变量因果连接网络。本发明不仅能够揭示大范围非平稳过程的变量因果关系，还具有对恶劣工业环境中噪声干扰的鲁棒性，以数据驱动的方式建立系统底层交互机制模型，为燃气轮机系统的状态评估、故障传播路径分析提供高效可靠的新思路。
Description
技术领域
本发明属于大范围非平稳系统的多变量因果分析领域，特别是涉及了一种将非平稳时间序列波动趋势转化为符号序列用于变量间信息传递量化度量的指标，通过时序重构的重采样方法设计因果显著性检验，以获得显著的多变量条件因果连接关系，构建了非平稳系统的因果网络图，实现对燃气轮机系统底层机理的因果分析。
背景技术
在使用大规模数据进行复杂工业系统分析以及故障传播识别的问题中，因果关系已被公认为是要挖掘的主要特征，因为它可以定性地表示过程单元的互联关系，揭示复杂系统的内部运作机理。不同于线性的相关关系，因果关系是有明确的时间先后顺序(先有因后有果)，同时信息传递方向确定的变量间耦合关系。在实际工业过程的设备工作原理非常复杂，许多变量之间具有复合链接关系，不能单纯通过线性表达式一以概之，很多维护检修的专家也无法明确整个系统各变量间的机理关系，所以需要根据运行数据建立可靠的因果网络直观表示系统的机理。此外，模型之间的耦合性使得一处故障的发生可以迅速在不同的工作环节传递，导致故障的影响不断扩大，最终造成严重的损失。精确的根因定位工作也可以在故障发生早期或甫一监测到异常信号的时候进行，通过分析信号之间的相似性关系，以及协方差等数学性质将最早发现异常的信号进行锁定，缩小关注范围。在故障早期，信号的异常波动往往非常微弱，很容易被噪声或正常振动掩盖，对于根因变量的确准有助于早发现早采取行动，对整体的发展趋势进行预判，能够最大可能降低故障的影响。
现有的因果溯源方法可以按照建模方式与适用对象分为基于系统模型的方法与数据驱动方法。
基于系统模型的方法起源较早，以现代控制理论为基础，建立描述整个系统关系的方程，根据实际连接关系确定变量之间的因果关系。此类方法适用于系统结构清晰，数学模型确定的情况，主要的研究方向是基于邻接矩阵与有向图进行模型的构建，进一步分析系统中各个变量之间的因果关系。对变量多，系统耦合关系复杂的变量无法进行有效的因果分析。
数据驱动的根因溯源方法随着工业大数据平台的发展成熟，逐渐受到更多人的重视，主要有贝叶斯网络，格兰杰因果分析，传递熵等方法。传递熵是一种基于过程历史定量的方法，分析结果是定性过程模型，其形式为有向图，显示变量之间的因果关系，在化工过程的扰动源分析领域应用广泛。格兰杰因果分析通过建立变量之间的预测模型，判断因变量的引入能否减小结果变量的预测误差，进而确定两个变量之间是否存在因果关系。贝叶斯网络则是以贝叶斯公式为基础，根据对于过程的先验知识建立网络结构，结合实际数据进行因果推断。上述方法通常与信号间的信息传递关系有关，可以实现在概率、预测精度、信息流方向等角度对于数据变量之间的因果关系进行合理推断与预测。现有的数据驱动方法通常假设变量符合高斯分布，用以对因果传递过程进行简化约束。
然而，在以燃气轮机为例的复杂工业系统中，由于操作切换、需求复合随着计划灵活调整、工作条件恶劣、环境噪声造成整个系统出现大范围非平稳现象，表现为变量波动范围大，波动时间久。非平稳现象最早由经济学领域的研究者提出，其数学表征是数据的分布随时间变化而改变，通常采用ADF检验判断一个时间序列是否属于平稳序列。基本的传递熵、格兰杰因果分析方法受到了弱平稳性假设的约束。如果将这些因果分析方法直接应用于非平稳变量的因果分析可能带来很多问题。一方面，非平稳变量可能具有相同的变化趋势，此时如果应用通过预测误差判断因果的方法极有可能受到共同趋势的影响，判断出的变量只是偶然具有相同的趋势，或者有未察觉到的中间变量导致了两个测试变量在数值上表现出共同趋势，这种因果通常被定义为虚假因果或者无效因果。另一方面，非平稳的变量方差波动大，会掩盖一些真实的因果信息，如果没有具有较好鲁棒性的方法进行因果分析，往往难以识别出显著的因果关系，或者识别出错误的因果关系。
近年来，关于变量间的非平稳因果分析的研究受到越来越多的关注，较为典型的有时变格兰杰因果分析、协整因果分析和符号转移熵。时变格兰杰因果分析要求已知系统的传递函数，对系统的先验知识要求较高，不适用于复杂多变量系统。基于协整分析的因果方法则要求非平稳时间序列必须是同阶单整的序列，很多实际的工业数据无法满足这一需求。其中符号转移熵因为其非参数建模的灵活特性，对非平稳序列和平稳序列的兼容性，成为非平稳因果分析的重要参考指标。但是一些现有的方法对因果显著性的判断过于简单，导致经常有冗余的因果关系产生，不便于对设备机理特性的分析。所以应该设计与因果指标相配合的显著性检验方法，形成清晰的系统因果网络。
发明内容
本发明的目的在于针对现有面向大范围非平稳多变量因果推断技术的不足，提供一种面向大范围非平稳过程的多变量因果网络构建方法。该方法通过引入一种滑窗式符号集合转化方法，将时间维度上的大范围非平稳变量重构为表示序列趋势的有限符号序列。在此基础上，该方法引入了条件符号转移熵，能够捕获多个变量间信息传递的特性，包括信息传递的量化评估和传递方向的基本确定。该方法针对符号传递熵的判断结果设计了时间洗牌的重采样方法进行因果关系显著性检验，进一步构建燃机系统的因果网络。本发明为大范围非平稳过程因果分析提供了一种新颖的分析角度，不仅跳出了原有方法对变量平稳性假设的约束，而且通过符号化的方式增强了因果分析的可靠性和可信度，有助于工业工程师对燃气轮机大范围非平稳瞬态连续过程运行状态做出准确判断，了解系统底层机理，在发生错误时追根溯源，从而保证燃气轮机发电过程的安全。
本发明的目的是通过以下技术方案来实现：一种基于燃气轮机大范围非平稳运行工况下的多变量因果分析方法，该方法包括以下步骤：
(1)获取在燃机系统非平稳瞬变连续运行过程的J个被测变量的N个样本，表示为二维数据矩阵X(N×J)；
(2)变量间符号转移熵计算，具体包括如下子步骤：
(2.1)设滑动窗口长度为lsw，以步长为s，从每个被测变量的N个样本序列的第一个时间点开始滑动，每个窗口取出的序列为一行，重构形成嵌入向量矩阵其中Nsw＝N-(lsw-1)×s，j∈[1，J]；矩阵的每一行代表一个嵌入向量；
(2.2)根据每个嵌入向量τ表示系统采样间隔，其中的观测值样本数据点按照从小到大的顺序设置0，1，...，(lsw-1)的符号/>将嵌入向量转化为秩向量/>获得符号矩阵/>其中，下标j表示变量序号，t表示时刻，lsw表示样本数据点在嵌入向量中的位置；
(2.3)计算任意被测变量Xm对变量Xn的条件符号转移熵：
其中，h为预测范围，为正整数；为系统中的其余被测变量的秩向量的拼接，n，m，k∈[1，J]且n≠m≠k，/>为变量的联合概率分布，为条件概率分布；
(3)根据步骤2得到与/>分别代表从变量Xn到变量Xm传递的信息量和反方向的信息量，计算因果方向指示量/>
其中，表示Xn是Xm的原因，/>则说明因果关系是从Xm到Xn。
进一步地，所述步骤1中，被测变量包括压气机出口压力、压气机出口温度、压气机轴承温度、压气机推力瓦轴承发电机端温度、压气机推力瓦轴承燃机端温度、压气机轴承振动、压气机侧大轴振动、压气机防冰冻装置进气电动调节阀阀位、燃机透平排气平均温度、燃机透平侧大轴振动、燃烧室压差、燃机嗡鸣、燃机功率和压气机压比。
进一步地，所述步骤2之前，还包括对系统非平稳程度估计步骤，具体如下：
对每个变量进行ADF平稳性假设检验，如果非平稳变量个数Jnon-sta≥J/3，则认为系统具有较为鲜明的大范围非平稳特性，使用本发明方法效果更佳继续执行步骤2，否则停止。
进一步地，还包括对被测变量间相关关系的初步筛选步骤，辅助验证因果关系的准确性。
具体如下：
对二维数据矩阵X进行标准化处理并计算两两被测变量之间的相关系数：
其中E(*)，D(*)分别代表被测变量的均值和方差；其中，corr(Xm，Xn)∈[-1，1]，正值代表被测变量间正相关，负值代表负相关，相关系数的绝对值越接近于1，被测变量间的线性相关性越强。
进一步地，所述步骤2中，联合概率分布和条件概率分布计算方法如下：
其中，为各自符号拼接形成的序列，/>为符号矩阵以及/>拼接而成的多元符号矩阵中符号序列/>出现的频次。
进一步地，所述步骤3中，当接近于0时，还包括如下步骤：
将被测变量Xn和Xm中任意时间点前后的观测值颠倒，打乱时间顺序，去除原有的因果自相关的依赖，得到经过洗牌后的替代序列和/>并利用步骤2方法计算替代序列的符号转移熵/>重复上述过程M次，得到替代时间序列的符号转移熵与原始数据计算得到的符号转移熵/> 组成符号转移熵集合/>该集合符合自由分布。
将因果显著性假设如下：
对于分布自由的单侧检验，检验的p-value计算公式为：
其中r0为在符号转移熵集合中按照值的大小升序排列得到的顺序；
根据pval修正Xn和Xm的因果关系：当pval≤0.10，有90％的概率拒绝H0假设，证明Xn和Xm有显著因果关系。pval值越小，因果关系越显著。
进一步地，还包括建立系统的变量间因果网络步骤，具体如下：
对系统中的两两变量分别进行因果显著性检验，有因果设为1，无因果设为0；遍历J个基本变量得因果关系邻接矩阵C(J×J)，其中行代表原因变量，列代表结果变量；
剔除因果关系邻接矩阵C中行列元素都为0的孤立变量点，得到与观测变量有因果关系的所有变量点；
添加因果连边：将因果关系邻接矩阵C中变量顺序进行调整，转化为下三角矩阵，按照方向链接因果关系为1的两个变量；
去除间接因果关系和冗余连边：若A对B，B对C有因果关系，则发现的A对C的因果连边将被去除，认为A到C是间接的因果关系。
与现有技术相比，本发明的有益效果是：本发明为大范围非平稳多变量系统的因果分析提供了新的研究思路。通过设计一种基于符号集的数据重构策略，首次将时间维度下的大范围非平稳序列巧妙地转化为了包含数据趋势信息的符号矩阵，从而为进一步的分析提供了基础；利用条件符号转移熵量化度量两个变量在系统现有观测条件下的信息传递方向及传递强度；为了避免传递熵值过小所产生的模糊意义和因果方向判断的偶然性，设计了适用于符号转移熵的因果显著性检验，利用时间洗牌的思路进行序列重采样，对原始变量间的传递熵值是否代表显著的因果关系做出判断，避免了冗余、虚假的因果关系。在此基础上，整理所有变量的因果关系形成变量间因果关系表示的二元矩阵，按照点-边-路的三个步骤进行因果图网络的构建。所提出的方法在实际工业过程中做了详细的实验研究，获得了成功应用。该方法通过对大范围非平稳序列的符号化映射，增强了对变量趋势信息的了解和对非平稳噪声扰动的鲁棒性，引入合适的多变量因果判据，配合自适应多变量因果网络构建方法，提高了因果分析对于大范围非平稳变化的适应性和准确性，最终可应用于燃气轮机联合循环发电现场，确保燃机运行过程的安全可靠和检修维护的高效便捷。
附图说明
图1是本发明所述非平稳多变量因果关系分析方法的流程图；
图2是本发明所述燃气轮机典型变量的标号即平稳性说明示意图；
图3是本发明所述典型非平稳变量原始波动趋势图；
图4是本发明所述燃气轮机主要变量的相关关系图；
图5是具体实施例中的本发明方法计算的两两变量间因果显著性结果图；
图6是具体实施例中的本发明方法推断的最显著的10个因果关系对结果图；
图7是本发明具体实施例中的本发明方法(a)的因果判断以及与传统格兰杰因果方法(b)的对比结果图。
图8是本发明建立的燃气轮机以压气机为中心观测的因果传递网络示意图(a)及实际燃气发电热力循环示意图(b)。
具体实施方式
下面结合附图及具体实例，对本发明做进一步说明。
燃气轮机发电机组是一个高度复杂的工业过程，具有时变，动态特性和非平稳性。发电用重型燃气轮机由压气机、燃烧室、燃气透平三大核心部件组成，其中压气机和透平为多级轴流设计，压气机共15级，压比16.9，透平4级，燃烧室为环型结构，逆时针布置，内置陶瓷遮热板，整个机组及燃料控制阀组布置在燃气轮机罩壳内。其工作原理为：压气机从外部吸收空气，空气从燃气轮机进气口进入，通过压气机多级叶片将其压力升高，压缩后送入燃烧室，同时燃料(气体或液体燃料)也喷入燃烧室与高温压缩空气混合，在定压下进行燃烧。生成的高温高压烟气燃烧受热后膨胀，进入透平区经过多级的叶片，推动动力叶片高速旋转，直至从出气口排出，成为废气，废气排入大气中或再加利用(如利用余热锅炉进行联合循环)。对燃气轮机而言，机组负荷的变化使得给煤量会不断根据需求进行调整，从而导致磨煤机的运行状态也不断改变，表现出典型的非平稳瞬态特性。
图1所示为本发明非平稳多变量因果网络构建的流程图，本发明方法包括以下步骤：
(1)获取待分析数据：在燃机系统非平稳瞬变连续过程中，共有J个变量被在线测量，共计采集了K个样本，得到直接采集的二维数据矩阵Xorigin(K×J)；
在本实例中，从浙江某发电厂实地采集了约1500个正常运行时的样本用于因果分析，测量变量为14个：压气机出口压力、压气机出口温度、压气机轴承温度、压气机推力瓦轴承发电机端温度、压气机推力瓦轴承燃机端温度、压气机轴承振动、压气机侧大轴振动、压气机防冰冻装置进气电动调节阀阀位、燃机透平排气平均温度、燃机透平侧大轴振动、燃烧室压差、燃机嗡鸣、燃机功率、压气机压比。
(2)数据预处理：
(2.1)停机数据清洗：对于步骤1所述原始二维数据矩阵Xorigin进行筛选，以燃机功率变量作为指示变量清除停机段的数据，得到实验数据矩阵X(N×J)；
在实例中去除停机数据后实际正常运行的数据有1024个，即N＝1024，J＝14。
(2.2)变量关系粗选：对实验数据X进行标准化处理计算两两变量之间的相关系数，作为对变量间相关关系的初步筛选，其中皮尔逊相关系数公式如下；
其中E(*)，D(*)分别代表变量的均值和方差。其中，corr(Xm，Xn)∈[-1，1]，正值代表变量间正相关，负值代表负相关，相关系数的绝对值越接近于1，变量间的线性相关性越强，初步关注变量间的基础关联，对系统中不同类型的变量的耦合关系进行初步分析。
(3)系统非平稳程度甄别，需要对每个变量进行ADF平稳性假设检验，该步骤由以下子步骤来实现：
(3.1)建立每个变量的向量自回归模型：
xt＝φ1xt-τ+…+φpxt-pτ+εt#(2)
其中，xt为某个变量第t个时间点的观测值，t＝(p+1)τ，(p+2)τ，...，Nτ，其中τ为系统采样间隔；φp为第p个自回归变量的回归系数，εt表示高斯随机噪声；p为正整数，一般取5-10之间，遍历回归获得回归系数。
(3.2)设定平稳性判断量：通过判断自回归方程系数之和是否为1来判断整体序列是否平稳，设计判断量ρ：
ρ＝φ1+φ2+…+φp-1#(3)
由单位根定理可以推证，平稳序列判断标准如下：
根据公式(4)设计假设检验：
(3.3)构造ADF假设t检验统计量为：
其中是随机打乱原始变量的时间顺序求得的判断量ρ与原始变量时间顺序求得的判断量/>的集合的标准差。通过蒙特卡洛法可以得到t检验统计的临界值表，查表得到在95％置信度下的假设结果，检验p-value大于0.05，维持H0假设认为变量属于平稳变量。
(3.4)系统非平稳程度估计：
对J个观测变量重复上述步骤(3.1)-(3.3)，以此确定每个变量的平稳性，如果非平稳变量个数Jnon-sta≥J/3，可认为系统具有较鲜明的大范围非平稳特性，使用发明方法效果更佳。
经过步骤3的非平稳检验，共有7个非平稳变量，数量超过观测变量的总数的三分之一，如附图2所示，可证明燃汽轮机具有大范围非平稳的运行特性。将部分非平稳变量随时间的变化趋势图画出如附图3所示，发现变量的方差和均值都随着时间剧烈波动，数据分布不恒定，数据之间没有稳定的共同趋势。附图4分析了变量的线性相关关系，对变量的基本联系进行分析。
(4)变量间符号转移熵计算：
(4.1)利用滑动窗口重新排列数据；
设滑动窗口长度为lsw，优选地，滑窗长度选择3-7较好。从变量序列的第一个时间点开始，以步长为s，依次滑动完整个序列，每个窗口取出的序列为一行，以第一列的变量序列即变量1为例，X1(N×1)经过滑窗重构形成嵌入向量矩阵其中Nsw＝N-(lsw-1)×s，矩阵的每一行代表一个嵌入向量；滑窗步长不超过滑窗长度的1/3，一般取1，本实施例取值为1。
(4.2)将嵌入向量矩阵符号化：
初始化符号集r＝{0，1，...，(lsw-1)}，对于某个时刻t对应的一个嵌入向量 其中τ表示系统采样间隔，按照嵌入向量中x观测值的大小升序排列并依次从小到大设置符号/>其中，下标1表示变量序号，t表示时刻，lsw表示观测数据点在嵌入向量中的位置；每个数据点对应符号值越大，其代表的原始数值也越大。由此，将嵌入向量转化为可以表示数据波动趋势的秩向量/>秩向量共有lsw！种可能的排列方式，组合成的符号排列集合记作/>遍历Nsw个嵌入向量，将嵌入向量矩阵转化为符号矩阵/>
(4.3)变量符号转移熵计算：对于J个变量依次进行步骤(4.1)-(4.2)的操作得到所有观测变量的符号序列。
变量X2对变量X1的条件符号转移熵公式为：
其中，h为正整数，表示最佳预测范围，通过多组h值择优获得，本实施中取值为1，为t时刻系统中的其余观测变量的秩向量的拼接，/>为变量的联合概率分布，/>为条件概率分布。
首先计算变量X1的边缘概率，统训的每种秩向量排列方式出现的概率，计算公式为：
其中为观测时间内排列/>出现的频次。
同理可得，联合概率分布的计算公式为：
其中为各自符号的拼接，/>为观测时间内即/> 以及/>拼接而成的多元符号矩阵中符号序列出现的频次。
利用公式(8)和(9)的结果，条件概率分布的计算公式依据贝叶斯原理，可以推导为：
由此可以计算出变量间的符号转移熵值。
根据对变量间传递熵的计算，可以去除在系统中其他变量对于所关注的变量的影响，这是传统的因果分析方法无法做到的。
(5)基于符号转移熵的因果判断指标与显著性检验，该步骤由以下子步骤来实现：
(5.1)因果方向判断指标：根据步骤4得到与/>分别代表从变量X1到变量X2传递的信息量和反方向的信息量，真实的因果方向信息传递更多，所以设计因果方向指示量/>
其中，表示X1是X2的原因，/>则说明因果关系是从X2到X1的；
(5.2)因果显著性检验：在两个方向的传递熵都接近于0的时候判断因果容易出现混淆，而且对于因果关系的显著性没有明确的判断。为此，引入因果显著性检验，保持其他条件变量不变，对因果变量随机取一个时间点trand，将此时间点前后的序列颠倒，打乱时间顺序，去除原有的因果自相关的依赖，得到经过洗牌后的替代序列同理可得/>利用步骤4计算替代序列的符号转移熵记为/>重复上述过程M次，得到替代时间序列的传递熵集合
进一步可得因果显著性假设：
将原始数据计算得到的传递熵，记作与打乱时间顺序的序列得到的传递熵组合，为自由分布。对于分布自由的单侧检验，检验的p-value计算公式为：
其中r0为在符号转移熵集合中按照值得大小升序排列得到的顺序。
当pval≤0.10，有90％的概率拒绝H′0假设，证明X1对X2有显著因果关系。pval值越小，因果关系越显著。
遍历系统中的两两变量分别进行因果显著性检验，有因果设为1，无因果设为0；遍历J个基本变量可得因果关系邻接矩阵C(J×J)，其中行代表结果变量，列代表原因变量，结果如附图5所示。
(6)建立系统的变量间因果网络，该步骤由以下子步骤来实现：
(6.1)简化因果关系邻接矩阵，初始化网络节点：首先剔除因果关系邻接矩阵中行列元素都为0的孤立变量点，得到与观测变量有因果关系的变量点；
(6.2)添加因果连边：将因果关系邻接矩阵中变量顺序进行调整，转化为下三角矩阵，便于网络连边添加，按照方向链接因果关系为1的两个变量；
(6.3)去除间接因果关系和冗余连边：当发现A对B，B对C有因果关系，再发现A对C的因果连边将被去除，认为A到C是间接的因果关系；结果如附图6所示。
根据燃汽轮机联合循环的发电设备设计，压气机、透平、发电机依靠中心大轴轴承相互连接，气流从压气机被压缩依次经过燃烧室和透平所以压气机侧大轴振动快慢影响透平侧大轴震动，即X6是X9的原因；振动的快慢代表燃机反应室内燃气燃烧的剧烈程度，振动越快燃机推力瓦轴承温度越高，而推力瓦轴承同时连接着发电机，按照反映的次序与热力循环的流程可知，如图7(a)本方法所判断的因果链路恰好代表已知观测变量在燃气轮机联合循环的主要因果路径。但是如图7(b)普通的格兰杰因果分析(Barnett，Lionel，and AnilK.Seth.″The MVGC multivariate Granger causality toolbox：a new approach toGranger-causal inference.〞Journal of neuroscience methods223(2014)：50-68.)反而更像受到了相关关系的诱导，判断出大量冗余不合乎机理的因果关系，例如变量X6(压气机侧大轴振动)、X7(压气机防冰冻装置进气电动调节阀阀位)并没有直接因果关系，但是格兰杰因果分析的结果误判了二者的因果关系。根据本发明的结果与普通格兰杰因果分析方法的结果进行对比后可以得出，本方法所监测的因果关系更加清晰，也更符合燃机发电机组的运行机理，值得注意的是本实施例以燃气轮机主轴为重点观测对象展开的，因果网络也是表现的燃机热力循环轴向的因果传递关系。如图8所示，根据本发明的因果关系检测，可以得出压气机轴向因果关系传递网络，符号化的方法提高了多变量系统的因果网络的可靠性，有助于工业工程师对过程运行状态做出准确判断，保证实际生产过程的安全可靠运行。
本发明所提面向燃气轮机大范围非平稳特性下的多变量系统因果分析方法，是在考虑了大范围非平稳时间序列的复杂波动特性的基础上，通过引入有限符号化方法进行变量的非平稳转换，最大限度地保留变量的趋势信息，进行变量间信息传递情况分析，进一步建立稀疏因果网络，实现对燃气轮机的多变量因果关系分析。基于条件段划分结果建立的细粒度监测系统可以为实际燃气轮机发电现场的技术管理部门提供清晰的因果分析网络，有助于运维人员快速了解燃气轮机的运行机理，为后续故障发生时迅速锁定故障根因，隔离传播路径提供重要参考，并最终为生产的安全可靠运行和产品的高质量追求奠定了基础。
应该理解，本发明并不局限于上述具体实施的燃气轮机运行过程，凡是熟悉本领域的技术人员在不违背本发明精神的前提下还可做出等同变形或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围。
Claims (5)
1.一种面向燃气轮机的非平稳多变量因果关系分析方法，其特征在于，包括以下步骤：
(1)获取在燃机系统非平稳瞬变连续运行过程的J个被测变量的N个样本，表示为二维数据矩阵X(N×J)；
(2)变量间符号转移熵计算，具体包括如下子步骤：
(2.1)设滑动窗口长度为lsw，以步长为s，从每个被测变量的N个样本序列的第一个时间点开始滑动，每个窗口取出的序列为一行，重构形成嵌入向量矩阵其中Nsw＝N-(lsw-1)×s，j∈[1，J]；矩阵的每一行代表一个嵌入向量；
(2.2)根据每个嵌入向量τ表示系统采样间隔，其中的观测值样本数据点按照从小到大的顺序设置0，1，...，(lsw-1)的符号/>将嵌入向量转化为秩向量/>获得符号矩阵/>其中，下标j表示变量序号，t表示时刻，lsw表示样本数据点在嵌入向量中的位置；
(2.3)计算任意被测变量Xm对变量Xn的条件符号转移熵：
其中，h为预测范围，为正整数；为系统中的其余被测变量的秩向量的拼接，n，m，k∈[1，J]且n≠m≠k，/>为变量的联合概率分布，为条件概率分布；
(3)根据步骤(2)得到与/>分别代表从变量Xn到变量Xm传递的信息量和反方向的信息量，计算因果方向指示量/>
其中，表示Xn是Xm的原因，/>则说明因果关系是从Xm到Xn；
当接近于0时，还包括如下步骤：
将被测变量Xn和Xm中任意时间点前后的观测值颠倒，打乱时间顺序，去除原有的因果自相关的依赖，得到经过洗牌后的替代序列和/>并利用步骤(2)方法计算替代序列的符号转移熵/>重复上述过程M次，得到替代时间序列的符号转移熵与原始数据计算得到的符号转移熵/> 组成符号转移熵集合/>该集合符合自由分布；
将因果显著性假设如下：
对于分布自由的单侧检验，检验的p-value计算公式为：
其中r0为在符号转移熵集合中按照值的大小升序排列得到的顺序；
根据pval修正Xn和Xm的因果关系：当pval≤0.10，有90％的概率拒绝H0假设，证明Xn和Xm有显著因果关系；pval值越小，因果关系越显著；
还包括建立系统的变量间因果网络步骤，具体如下：
对系统中的两两变量分别进行因果显著性检验，有因果设为1，无因果设为0；遍历J个基本变量得因果关系邻接矩阵C(J×J)，其中行代表原因变量，列代表结果变量；
剔除因果关系邻接矩阵C中行列元素都为0的孤立变量点，得到与观测变量有因果关系的所有变量点；
添加因果连边：将因果关系邻接矩阵C中变量顺序进行调整，转化为下三角矩阵，按照方向链接因果关系为1的两个变量；
去除间接因果关系和冗余连边：若A对B，B对C有因果关系，则发现的A对C的因果连边将被去除，认为A到C是间接的因果关系。
2.根据权利要求1所述一种面向燃气轮机的非平稳多变量因果关系分析方法，其特征在于，所述步骤(1)中，被测变量包括压气机出口压力、压气机出口温度、压气机轴承温度、压气机推力瓦轴承发电机端温度、压气机推力瓦轴承燃机端温度、压气机轴承振动、压气机侧大轴振动、压气机防冰冻装置进气电动调节阀阀位、燃机透平排气平均温度、燃机透平侧大轴振动、燃烧室压差、燃机嗡鸣、燃机功率和压气机压比。
3.根据权利要求1所述一种面向燃气轮机的非平稳多变量因果关系分析方法，其特征在于，所述步骤(2)之前，还包括对系统非平稳程度估计步骤，具体如下：
对每个变量进行ADF平稳性假设检验，如果非平稳变量个数Jnon-sta≥J/3，则认为系统具有较为鲜明的大范围非平稳特性，继续执行步骤(2)，否则停止。
4.根据权利要求1所述一种面向燃气轮机的非平稳多变量因果关系分析方法，其特征在于，还包括对被测变量间相关关系的初步筛选步骤，辅助验证因果关系的准确性；具体如下：
对二维数据矩阵X进行标准化处理并计算两两被测变量之间的相关系数：
其中E(*)，D(*)分别代表被测变量的均值和方差；其中，corr(Xm，Xn)∈[-1，1]，正值代表被测变量间正相关，负值代表负相关，相关系数的绝对值越接近于1，被测变量间的线性相关性越强。
5.根据权利要求1所述一种面向燃气轮机的非平稳多变量因果关系分析方法，其特征在于，所述步骤(2)中，联合概率分布和条件概率分布计算方法如下：
其中，为各自符号拼接形成的序列，/>为符号矩阵以及/>拼接而成的多元符号矩阵中符号序列/>出现的频次。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110863191.4A CN113656906B (zh) | 2021-07-29 | 2021-07-29 | 一种面向燃气轮机的非平稳多变量因果关系分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110863191.4A CN113656906B (zh) | 2021-07-29 | 2021-07-29 | 一种面向燃气轮机的非平稳多变量因果关系分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113656906A CN113656906A (zh) | 2021-11-16 |
| CN113656906B true CN113656906B (zh) | 2023-10-03 |
Family
ID=78490844
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110863191.4A Active CN113656906B (zh) | 2021-07-29 | 2021-07-29 | 一种面向燃气轮机的非平稳多变量因果关系分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113656906B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP4231108A1 (en) | 2022-02-18 | 2023-08-23 | Tata Consultancy Services Limited | Method and system for root cause identification of faults in manufacturing and process industries |
| CN114240264A (zh) * | 2022-02-24 | 2022-03-25 | 成都四方伟业软件股份有限公司 | 一种城管事件指标间的因果关系检验方法及装置 |
| CN114881124B (zh) * | 2022-04-21 | 2023-07-25 | 北京百度网讯科技有限公司 | 因果关系图的构建方法、装置、电子设备和介质 |
| CN114925330B (zh) * | 2022-05-19 | 2023-03-21 | 自然资源部第一海洋研究所 | 一种不依赖重构参数的自适应跨尺度因果分析方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102541017A (zh) * | 2012-01-12 | 2012-07-04 | 浙江大学 | 一种复杂化工过程中振荡信号的快速定位方法 |
| CN104359685A (zh) * | 2014-11-24 | 2015-02-18 | 沈阳化工大学 | 一种柴油机故障识别方法 |
| CN108627345A (zh) * | 2018-05-11 | 2018-10-09 | 浙江师范大学 | 一种汽轮机系统级故障的诊断方法及系统 |
| CN109359662A (zh) * | 2018-08-20 | 2019-02-19 | 浙江大学 | 一种面向百万千瓦超超临界机组非平稳特性的基于因果分析的多层贝叶斯网络故障诊断方法 |
| CN110755062A (zh) * | 2019-10-29 | 2020-02-07 | 电子科技大学 | 基于符号转移熵的生理器官网络非平衡性量化分析方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7599840B2 (en) * | 2005-07-15 | 2009-10-06 | Microsoft Corporation | Selectively using multiple entropy models in adaptive coding and decoding |
| CN110390396B (zh) * | 2018-04-16 | 2024-03-19 | 日本电气株式会社 | 用于估计观测变量之间的因果关系的方法、装置和系统 |
-
2021
- 2021-07-29 CN CN202110863191.4A patent/CN113656906B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102541017A (zh) * | 2012-01-12 | 2012-07-04 | 浙江大学 | 一种复杂化工过程中振荡信号的快速定位方法 |
| CN104359685A (zh) * | 2014-11-24 | 2015-02-18 | 沈阳化工大学 | 一种柴油机故障识别方法 |
| CN108627345A (zh) * | 2018-05-11 | 2018-10-09 | 浙江师范大学 | 一种汽轮机系统级故障的诊断方法及系统 |
| CN109359662A (zh) * | 2018-08-20 | 2019-02-19 | 浙江大学 | 一种面向百万千瓦超超临界机组非平稳特性的基于因果分析的多层贝叶斯网络故障诊断方法 |
| CN110755062A (zh) * | 2019-10-29 | 2020-02-07 | 电子科技大学 | 基于符号转移熵的生理器官网络非平衡性量化分析方法 |
Non-Patent Citations (2)
| Title |
|---|
| 郜传厚 ; 渐令 ; 陈积明 ; 孙优贤 ; .复杂高炉炼铁过程的数据驱动建模及预测算法.自动化学报.2009,(第06期),全文. * |
| 黄信林 ; 高建民 ; 高智勇 ; 吕延庆.应用变量因果序分析的符号有向图建模方法.西安交通大学学报. 自然科学版.2010,第44卷(第5期),全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113656906A (zh) | 2021-11-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113656906B (zh) | 一种面向燃气轮机的非平稳多变量因果关系分析方法 | |
| Li et al. | Gas turbine gas path diagnosis under transient operating conditions: A steady state performance model based local optimization approach | |
| Hanachi et al. | A physics-based modeling approach for performance monitoring in gas turbine engines | |
| EP2495631B1 (en) | A system for analysis of turbo machinery | |
| US8751423B2 (en) | Turbine performance diagnostic system and methods | |
| Mathioudakis et al. | Performance analysis of industrial gas turbines for engine condition monitoring | |
| Chen et al. | A time-series turbofan engine successive fault diagnosis under both steady-state and dynamic conditions | |
| CN114757380A (zh) | 一种火电厂故障预警系统、方法、电子设备及存储介质 | |
| CN109359662A (zh) | 一种面向百万千瓦超超临界机组非平稳特性的基于因果分析的多层贝叶斯网络故障诊断方法 | |
| Loboda | Gas turbine condition monitoring and diagnostics | |
| Feng et al. | Gas turbine blade fracturing fault diagnosis based on broadband casing vibration | |
| CN115618592A (zh) | 一种电厂燃气轮机气路故障诊断方法、系统、设备及终端 | |
| Adamowicz et al. | Advanced gas turbines health monitoring systems | |
| Hu et al. | Intelligent condition assessment of industry machinery using multiple type of signal from monitoring system | |
| CN114936520A (zh) | 基于改进mset的汽轮机故障预警分析方法 | |
| Wang et al. | Research on anomaly detection and positioning of marine nuclear power steam turbine unit based on isolated forest | |
| Li et al. | Gas turbine gas‐path fault diagnosis in power plant under transient operating condition with variable geometry compressor | |
| Mathioudakis et al. | Uncertainty reduction in gas turbine performance diagnostics by accounting for humidity effects | |
| Song et al. | Improved CEEMDAN-based aero-engine gas-path parameter forecasting using SCINet | |
| CN115422714A (zh) | 面向燃气轮机的知识条件混合驱动的运行状态监测方法 | |
| CN115526238A (zh) | 一种基于全工况卷积特征记忆的引风机故障预警方法 | |
| Eustace | A real-world application of fuzzy logic and influence coefficients for gas turbine performance diagnostics | |
| Hanachi | Gas turbine engine performance estimation and prediction | |
| Hanachi et al. | A physics-based performance indicator for gas turbine engines under variable operating conditions | |
| Goyal et al. | Prediction Enhancement of Machine Learning Using Time Series Modeling in Gas Turbines |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |