基于PCA/MSET联用模型的烟气轮机故障预警研究*

张伟亚，陈文武，韩 磊，潘 隆

(中石化安全工程研究院有限公司化学品安全控制国家重点实验室，山东青岛 266104)

0 前言

烟气轮机是催化裂化装置的关键机组设备，一旦发生故障，将影响整个催化裂化反应的进行[1]。烟气轮机的驱动介质主要来自催化裂化装置沉降器，经沉降器沉降后的高温烟气仍含有部分催化剂，导致烟气轮机在工作状态下处于高温、高转速、催化剂冲蚀、催化剂沉积结垢等恶劣环境，造成烟气轮机故障频发。因此，开展烟气轮机故障预警，对于预防性维护和保障设备长周期安全运行具有重大意义。

针对装置设备进行故障诊断，根据故障特征表现形式和提取来源的不同，大致可分为基于数学模型的故障诊断、基于经验知识的故障诊断和基于数据驱动的方法[2]。其中，基于数据驱动的故障诊断方法不需要精确数学模型，经验知识应用较少，仅通过过程数据即可对故障进行诊断，具有广泛的应用[3]。主元分析方法(Principal component analysis, PCA)是基于数据驱动故障检测的常用方法之一，适合复杂线性系统的故障检测。为反映模型主元空间和残差空间的模态变化，引入T2统计量和Q统计量，用于故障检测和诊断，然而根据T2和Q统计量超限情况判别是否出现故障存在结果不确定性，因此李传坤，等[2]采用新的综合故障监测指标实现催化装置智能预警。此外，韩万里，等[4]利用PCA主元分析，进行数据降维，结合多元状态估计(Multivariate State Estimation Technique，MSET)，实现引风机的故障早期预警。金秀章[5]发现PCA结合MSET模型对照相关性分析数据降维方法，预警效果更好。MSET是一种非线性的多变量预测估计算法，与其他方法相比具有满足现场实时工况、建模简单等优点。为实现压气机故障预警，陆永卿[6]建立了基于MSET和偏离度的故障预警非参数模型。同样，李大中[7]针对风力发电机中监测参数建立MSET预警模型，并验证了模型有效性。

综上，PCA结合MSET模型在风机设备监测方面可以有效实现监测预警，然而设备特性不同，信号采集频率也不尽相同，本文提出将PCA和MSET相结合的方式应用到催化裂化装置烟气轮机故障预警方面，验证模型的可适性。首先，采用PCA方法提取烟气轮机正常工况下运行的历史振动数据，开展时频域分析和特征计算，并基于筛选后特征构建MSET模型中历史记忆矩阵，计算预测数据与观测数据之间的残差，实现设备故障预警。

1 PCA数据降维

PCA方法本质是通过坐标系旋转和平移转换，使得样本的重心与新坐标系原点重合，同时新坐标系第一轴与数据变异最大的方向平行，第二轴与数据变异第二大方向平行，以此类推。摒弃携带信息量较小的主轴后，剩余k个主轴可以有效反馈原n个变量的情况，从而实现数据降维的目的。

假设烟气轮机正常运行状态下共采集m个数据和n个参数变量，则样本集合矩阵Xm×n如下：

(1)

首先，对Xm×n进行z-score标准化处理：

(2)

(3)

(4)

其中，i的范围为1～m，j的范围为1～n。

标准化后样本集合矩阵Zm×n为：

(5)

然后，计算样本集合矩阵Zm×n的相关系数矩阵L：

(6)

(7)

(8)

其中，k和l分别对应样本集合矩阵Zm×n第k列和第l列。

相关系数矩阵L定义如下：

(9)

对相关系数矩阵L进行特征分解，得到特征值λ和特征向量P。

对特征值λ上下反转使其从大到小排列，并计算各主元累积贡献率Ck：

(10)

若累积贡献率大于90%，则不再增加主元个数，从而实现数据降维。

新的主元计算公式如下：

t=Zm×nPn×k

(11)

为衡量包含在主元模型中的信息大小，引入T2(Hotelling)统计量，对于样本个数为m，主元个数为k的变量，T2统计量服从自由度为k和n-k的F分布，置信度为0.99的T2统计量控制上限为：

(12)

2 MSET数据建模

非参数建模方法MSET本质是将健康状态采集的历史数据与全部监测数据进行对比，当设备性能发生退化，由健康状态训练样本构成的记忆矩阵得到的各主元估计值精度发生变化，与监测对象之间的残差增大，因此采用残差的大小可以反映设备的运行状态。

假定设备k个统计参数在j时刻的观测向量表示为：

(13)

正常历史运行状态采集到m′个样本的观测向量构成记忆矩阵D：

(14)

将MSET模型的输入定义为观测向量Xobs，输出定义为估计向量Xest。估计向量Xest可由记忆矩阵D计算获得：

Xest=DW

(15)

其中W为观测向量对应的m维权值向量，表示当前观测向量与历史状态样本的相似度，表示为：

W=[w1w2wj…wm]

(16)

权值向量W可通过模型输入观测值Xobs与输出估计值Xest之间的残差确定：

ε=Xobs-Xest

(17)

令约束条件‖ε‖2最小，权值向量可表示为：

(18)

由于构造记忆矩阵D的观测向量数远大于实时监测参数量，导致运算不可逆，因此采用非线性运算符代替普通矩阵乘积运算：

(19)

非线性运算距离度量采用欧式距离运算，得到观测向量估计值为：

(20)

3 烟气轮机故障预警案例应用

PCA结合MSET模型在风机设备监测方面可以有效的实现监测预警，然而设备特性不同，信号采集频率也不尽相同，因此将PCA和MSET相结合应用到催化裂化装置烟气轮机故障预警方面，需验证模型的可适性。

3.1 计算流程

基于PCA和MSET相结合的计算方法流程如下：

伴随信息化和现代化的不断发展，我国保险公司的财务风险管控也越来越受到社会关注，这不仅关系着国民经济的发展，也牵扯着社会关系的和谐与否。因此，建立完整的财务风险预警体系是当务之急。利用财务风险预警体系，可以在实际运行中为保险公司提供实时的风险动态，同时也能对潜在的风险进行预警通报以做好防控措施。可以有效降低财务风险发生的可能性，同时也能在风险真的发生之时以最有效的措施应对，减少因财务风险引起的企业财务损失。

a) 采集烟气轮机振动信号，正常运行状态数据作为训练样本。

b) 提取训练样本时频域特征参数，包含均值、均方根值、方根幅值、绝对平均值、偏斜度、峭度、方差、最大值、最小值、峰峰值、波形指标、脉冲指标、裕度指标、偏斜度指标、峭度指标、均值频率等共计29个特征。

c) 对提取到的时频域特征进行主成分分析，计算每一变量的贡献量，选取累积贡献量达到90%时的主要成分，实现数据降维。

d) 计算置信度为0.99的T2统计量控制上限。

e) 将筛选后的特征参数作为列向量构建记忆状态矩阵D，将正常到故障状态的振动信号特征参数作为观测向量输入记忆矩阵，输出观测向量的估计值，计算两者之间的残差值。

f) 通过比较观测向量和估计向量之间的差值大小来判别烟机是否发生故障，实现烟机故障预警。

3.2 数据获取

烟机自2020年12月存在振动值上涨现象，相位同步发生变化，分析转子存在渐变不平衡，2021年3月机组停机检修进行清垢后恢复正常，转子结垢现象如图1所示。

图1 烟气轮机转子结垢

3.3 结果与分析

图2(a)表示烟机联端的时域波形数据，采集时间跨度为28 056×5 min，每个时间采集点包含1 024个数据，共计2.87×107个数据，数据量过于庞大，给分析造成困难，因此有必要进行数据降维分析。

目前企业在线监测系统进行设备故障预警主要采取固定阈值报警的策略，当振动峰峰值达到一定值时即发生报警，通常高报值设定为80 μm，高高报值设定为100 μm。从图2(b)中可以看出，在时间点3 283～3 287出现5次连续超过高报阈值线，定义为报警点，此时推断发生催化剂结垢脱落，转子质量平衡发生变化，反馈为振动值发生突变。在后续运行相当一段时间(66 d)，振动值仍然超过80 μm，随着振动磨损加剧，在22 435时间点出现激增趋势，停机检修处理。如果催化剂结垢大面积脱落，严重者会损伤叶片，因此实现超前预警具有重要意义。

图2 烟机联端振动时域波形和峰峰值随时间变化

根据3.1计算流程，选取正常数据m=2 000个时间点，特征参量n=29进行降维分析。通过标准化处理后建立样本集合矩阵Zm×n并计算相关系数矩阵L，特征分解后得到从大到小排列的特征值λ和特征向量P。图3(a)表示29个特征值对应的贡献率。特征值从大到小进行排列，前6个特征值所占的比重已超过90%，可认为前6个特征值对应主元包含的信息可以反馈29个特征参量的变化。计算载荷矩阵，选取相关系数较大的前6个主元参量进行分析。主元1与特征参量6、特征参量11和特征参量27的相关系数为0.51，0.52，0.34，是峭度、波形指标和频谱集中程度的综合表征。主元2与特征参量22、特征参量26和特征参量29的相关系数为0.32，0.32，0.25，均反映了频谱分散或集中程度。主元3与特征参量15、特征参量19和特征参量20的相关系数为0.30，0.40，0.40，是偏斜度指标和频谱分散或集中程度的综合表征。主元4与特征参量15、特征参量19和特征参量20的相关系数为0.35，0.37，0.40，反映了偏斜度指标和频谱分散或集中程度的综合表征。主元5与特征参量6、特征参量17和特征参量25的相关系数为0.26，0.24，0.28，峭度、均值频率和主频带位置的综合表征。主元6与特征参量2、特征参量7和特征参量18的相关系数为0.29，0.29，0.29，均方根值、方差和频率标准差的综合表征。通过公式(11)计算，可得到新的6个主元随时间变化趋势，如图3(b)所示相对其他主元，主元6在发生故障阶段波动较大，可作为重点监测的对象。

图3 PCA主元特征值贡献率和主元特征参数值随时间变化

正常数据量的选择往往只能定性分析，为探讨数据量的选择是否满足分析要求，图4(a)计算了不同时间点选择情况的T2值变化。从图中可以看出，在时间点数据量从1 000增加至1 750区间，主元99%可信区间阈值从24.73缓慢降至24.55。在数据量1 750～2 000期间，发生了陡降。从峰峰值(图2(b))时间波形图中可以看出，在时间点1 871处发生了振动下降现象，从而影响了阈值幅值。从2 000增加至3 000区间，主元99%可信区间阈值从22.62缓慢降至22.57。时间点数据量2 000满足分析要求。

图4 PCA T2值和MSET残差值随时间点变化

图4(b)表示基于MSET模型计算的各主元观测值与估计值之间的残差。从图中可以看出，由于记忆矩阵选取前2 000个时间点，观测值和估计值之间的残差为0。在时间点2 000～3 000之间，残差开始出现，但设备处于正常运转范围内，残差较小。以主元6监测对象为例，残差阈值设定为5，在时间点3 181～3 185连续5次超过残差阈值，相对于固定高报报警时间点提前102个时间点，超前时间8.5 h，可以实现烟气轮机超前预警。

4 结论

烟气轮机作为催化裂化装置重要的节能设备，其可靠性和稳定性直接影响催化裂化装置的安全性和经济性。采用基于PCA和MSET模型相结合的方式对烟气轮机进行故障预警，该方法不仅可以降低数据维度，同时构建的多维状态估计模型具有建模速度快，精度较高的优点，可以实现烟气轮机超前预警。以某催化裂化装置烟气轮机真实振动监测数据进行模型验证，实现超前预警8.5 h。