基于无监督特征学习的水电机组健康状态实时评价方法

胡 晓，肖志怀，2，刘 东，吴道平，查海涛，廖志芳

（1.武汉大学 动力与机械学院，湖北 武汉 430072；2.武汉大学 水力机械过渡过程教育部重点实验室，湖北 武汉 430072；3.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；4.国网江西省电力有限公司电力科学研究院，江西 南昌 330096；5.国网江西省电力有限公司柘林水电厂，江西 南昌 330096；6.天津市输水系统水锤阀门控制技术企业重点实验室，天津 300051）

1 研究背景

水电机组健康状态影响着整个水电站的正常运行，随着机组不断向大型化、复杂化方向发展，机组健康状态实时评价的重要性日益突出［1］。为预防水电机组故障，水电行业开发了状态监测系统，在线收集机组运行数据，能实现越限报警及分析诊断等功能。越限报警通常采用状态监测值与相关行业标准规定的单一限值比较，存在不能体现机组个性特色、适用机组运行工况有限的缺点［2］。为弥补传统方法的不足，充分挖掘机组状态监测数据中的有效信息，近年来，通过信号分析方法提取水电机组状态特征成为水电机组健康评价与故障诊断领域研究热点。例如，采用快速傅里叶变换（Fast Fourier Transform，FFT）、小波变换、经验模态分解等方法将信号由时域转换到频域或时频域，再计算峭度，熵值等参数作为特征输入神经网络，支持向量机等分类器进行故障诊断［3-4］。这些方法通常需要借助人工给出标签或者反馈等指导信息来选取最能反映故障的特征。但在水电机组实际运行中，故障概率较小，且故障机理复杂，故障类型不可能全部被预先掌握，用于鉴别一种故障的特征参数可能对其他故障不敏感。因此，故障样本不足，特征选择耗时耗力降低了有监督特征学习方法的实用性。

鉴于水电机组状态监测系统积累了大量机组正常运行数据，从研究机组正常状态特征入手，通过建立健全的正常样本特征库来实现水电机组健康状态实时评价具有很强的工程实用性［5］。无监督特征学习是机器学习方法之一，它能够从无标签训练数据中挖掘出有效特征或表示，发现数据中隐藏的有用信息，在水电机组信号特征提取方面有极大应用潜力。

奇异值分解（Singular Value Decomposition，SVD）和自编码器是无监督特征学习领域的重要方法，可有效提取信号特征。SVD 将包含信号信息的矩阵分解为一系列子矩阵和对应的奇异值，两者反映了矩阵蕴含的信息量，且奇异值具有比例不变性和旋转不变性，是较为稳定的信号特征［6-8］。自编码器（Auto-Encoder，AE）通过无监督方式来学习一组数据的关键特征，采用正常信号训练得到的自编码器，其网络参数和特性都只适合重构正常信号，当有异常信号输入时，模型重构误差将发生较大变化，通过分析自编码器重构误差的变化，可确定当前输入信号是否属于正常状态［9-10］。

为了充分利用水电机组正常数据，量化评价机组健康状态，本文提出了基于奇异值分解的时域劣化指标和基于自编码器重构误差的频域劣化指标。同时，考虑到单纯的时、频域劣化指标描述水电机组状态变化不够全面，提出了融合时、频域劣化信息的综合劣化指标，以反映机组健康状态的劣化趋势和劣化程度，并采用实测数据验证了所提方法的有效性。

2 基于典型工况划分的水电机组健康状态实时评价

水电机组运行工况复杂，且机组状态监测数据与运行工况与紧密相关［2］，因此在构建劣化指标时应考虑工况差异，对不同工况建立不同标准特征库，通过计算机组特征实际值与同一工况标准值的差异，量化评价机组当前健康状态。基于典型工况划分的水电机组健康状态实时评价分为离线和在线两个阶段［11］，具体流程如图1所示。

在离线阶段，从水电机组状态监测系统中选取机组运行状态良好时段内采集的数据建立含工况参数的正常记录，同时划分机组典型工况，建立典型工况正常记录库。当两个工况的水头之差ΔH和有功之差ΔP足够小时，可近似作为相同工况，据此可将水电机组运行工况按照水头、有功划分为M个典型工况。对于某一典型工况（Pm，Hm）将所有符合 |pi -Pm|＜ΔP/2 且|hj -Hm|＜ΔH/2 的（pi，hj，xk）处理成（Pm，Hm，xk）。根据国标GB/T 32584-2016 中定义水电机组稳态运行工况为流量（有功）、水头等主要参数保持在±1.5%额定值以内，可取ΔH= 3%额定水头、ΔP= 3%额定有功［12］。之后将各典型工况正常记录库中的数据分成两部分，分别构成标准数据子库Sm和正常数据子库Nm。利用Sm中的数据计算得到标准时、频域特征，利用Nm中的数据分别求出时、频域差异特征的置信区间和正常阈值。

在线阶段，首先根据pi和hj求出机组当前归属的典型工况，分别从标准数据子库调出标准时、频域特征，从正常数据子库中调出时、频域差异特征正常阈值。然后结合标准时，频域特征计算当前状态监测信号xk的时域差异特征tdfk和频域差异特征fdfk，根据正常阈值计算当前状态监测信号xk的时域劣化指标TI和频域劣化指标FI。最后由TI和FI加权计算得到综合劣化指标CI，CI值小于等于1时，表明机组处于健康状态；CI值越大，机组健康状态劣化程度越高。

3 基于无监督特征学习的时、频域劣化指标及综合劣化指标

3.1 时域劣化指标时域劣化指标是在时间域内对信号进行特征提取，研究对象是信号幅值。对包含信号幅值信息的矩阵进行SVD，可将原矩阵分解到不同子空间中，突出信号本质特征。

3.1.1 基本原理 SVD 相当于将1 个秩为r的实矩阵Am×n分解为r个秩为1 的m×n阶矩阵的加权和，权值即为矩阵A 的奇异值，如下式所示：

式中：U、V 为酉矩阵，且U ∈Cm×m，V ∈Cn×n；λi为ATA 的特征值，λi ＞0 ；vi为A 的右奇异向量； ui为A 的左奇异向量； D 是一个对角矩阵，且［13］。

图1 基于典型工况划分的水电机组健康状态实时评价

由式（1）可知，经过SVD 得到的一系列子矩阵及其对应的奇异值包含了矩阵的信息，在一定程度上代表了矩阵的本质特征，故在模式识别中，矩阵奇异值常被作为信号特征量［8］。

SVD 的步骤为：（1）将矩阵A 和A 的转置做矩阵乘法，得到ATA 和AAT；（2）对ATA 进行特征分解，得到n个特征值和对应的n个特征向量，将ATA 的所有特征向量张成一个n×n的矩阵，即式（1）中的V 矩阵，V 中的每个特征向量叫做A 的右奇异向量，记为vi；（3）对AAT进行特征分解，得到m个特征值和对应的m个特征向量，将AAT的所有特征向量张成一个m×m的矩阵，即式（1）中的U 矩阵，U中的每个特征向量叫做A 的左奇异向量，记为ui；（4）根据σi =Avi/ui，求出每个奇异值σi，进而求出奇异值矩阵D。

3.1.2 基于特征空间重构SVD 的时域特征提取 水电机组状态监测信号通常为一维时间序列，而SVD 的处理分析对象是矩阵，为构造SVD 的输入矩阵，需要对时域信号进行特征空间重构。假设X=[x1，x2，…，xL]表示一个给定信号序列，L为数据长度，定义X的特征空间重构矩阵为：

式中：p、q均为正整数，

为了确定p、q的最佳取值，引入类内距离评估指数作为指标［14］，类内距离评估指数越大，说明同类样本间的差异越大，样本越分散；反之，样本越集中。

以典型工况（Pm，Hm）为例，假设该工况标准数据子库Sm中有M个信号，其中第i个信号的奇异值为si =[si(1)si(2)…si(p)]，第j个信号的奇异值为sj =[sj(1)sj(2)…sj(p)]，且i≠j，类内距离评估指数算法步骤为：

（1）计算Sm中不同信号奇异值的欧氏距离，如下式所示：

然后将该值对所有样本求平均值，即：

（2）计算的变差因子，如下式所示：

（3）计算类内距离评估指数，如下式所示：

确定p的最佳值后，将一维时域信号重构为p×q的矩阵H，对H 进行SVD，得到p个奇异值，作为信号时域特征。

3.1.3 时域差异特征 运行条件波动，测量仪器灵敏度以及随机干扰都会导致水电机组状态监测信号在测量时出现随机误差。因此，同一工况下，即使机组始终处于正常运行状态，信号特征也会发生改变，但这种变化较小，而机组运行状态发生异常时，信号特征会发生较大改变。为了量化特征变化程度，反映机组健康状态，本文提出了时域差异特征。以典型工况（Pm，Hm）为例，记该工况标准数据子库为Sm，正常数据子库为Nm，时域差异特征计算步骤如下。

（1）对Sm中各样本进行中心化处理。中心化公式下式所示：

式中：X=[x1，x2，…，xL]表示一个给定时域信号；L为数据长度；X *为中心化处理后的时域信号；

（2）根据类内距离评估指数确定p的最佳值，并计算Sm中每个样本的时域特征。

（3）采用聚类算法计算标准时域特征。K 均值聚类算法以最小化簇内所有样本到簇均值向量的欧氏距离之和为目标寻找聚类中心，具有快速，实现简单等优点［15］，故本文采用K 均值聚类算法求出Sm中全部时域特征值的聚类中心作为标准时域特征

（4）根据步骤（1）、（2）计算Nm中各样本的时域特征。然后根据式（3）计算Nm中各样本时域特征与标准时域特征的欧氏距离作为时域差异特征。

3.1.4 时域劣化指标 工况不变，机组正常运行时，时域差异特征值会在小范围内波动，如何确定正常状态与异常状态的时域差异特征阈值是对机组进行健康状态评价的关键。由统计学方法可知，对于一个有确定样本个数的未知分布，可用置信区间表示样本参数真实值以一定概率落在测量结果周围的程度［16-17］。根据中心极限定理，对于包含K个样本时域差异特征的集合Td=[tdf1，tdf2，…tdfK]， 如果K足够大的话，那么时域差异特征取值的置信区间可以根据下式计算获得［18］：

以典型工况（Pm，Hm）为例，根据式（8）计算得到正常时域差异特征置信区间上限为tdfθ，假设该工况待测样本时域差异特征值为tdfm，若样本为正常数据，那么tdfm应小于tdfθ，故将tdfθ作为时域差异特征正常阈值。定义时域劣化指标TI，计算公式如下式所示：

在离线阶段，利用水电机组状态监测系统中累积的海量正常数据计算各典型工况下时域差异特征正常阈值，在线阶段，根据式（9）计算时域劣化指标，即可实时评价机组健康状态。

3.2 频域劣化指标为了显示出某些时域上难以观察到的信号特征，可对信号进行FFT 得到其频谱，然后利用自编码器的无监督特征学习能力提取信号频谱中的有效信息。

3.2.1 基本原理 假设有一组B维样本y(n)∈ℜB，1 ≤n≤N，N为样本总数，自编码器将这组数据映射到特征空间，得到每个样本的编码z(n)∈ℜW，1 ≤n≤N，并且希望这组编码可以重构出原来的样本［19］。自编码器的结构可分为两部分：

（1）编码器：

（2）解码器：

自编码器的学习目标是最小化重构误差：

如果特征空间维度W小于原始空间维度B，自编码器相当于是一种降维或特征抽取方法。最简单的自编码器是如图2所示的两层神经网络，输入层到隐藏层用来编码，隐藏层到输出层用来解码，层与层之间互相全连接。

图2 两层网络结构的自编码器

3.2.2 频域差异特征 以典型工况（Pm，Hm）为例，记该工况标准数据子库为Sm，正常数据子库为Nm，频域差异特征的计算步骤为：

（1）对Sm中的样本进行中心化处理。中心化公式如式（7）所示。

（2）对Sm中经过中心化处理后的各样本进行FFT，得到对应的频谱。

（3）将步骤（2）得到的所有频谱划分为训练集和验证集输入自编码器进行训练，得到验证集损失函数收敛的自编码器。

（4）对于Nm中的样本，根据步骤（1）、（2）得到对应频谱，然后输入步骤（3）训练好的自编码器，按照式（13）计算重构误差re，并将re作为各样本频域差异特征。

式中：d为频谱F=[f1，f2，…，fd]的数据长度，F′=[f′1，f′2，…，f′d]为自编码器的输出。

3.2.3 频域劣化指标 按照计算时域劣化指标的思路，定义频域劣化指标FI，计算公式如下式所示：

式中：fdfθ为典型工况（Pm，Hm）的频域差异特征正常阈值；fdfm为典型工况（Pm，Hm）的待测样本频域差异特征值。

3.3 综合劣化指标时域劣化指标和频域劣化指标从不同角度评价机组运行状态，都能有效地反映机组健康状态劣化趋势和劣化程度。然而，只采用时域劣化指标或频域劣化指标不能充分反应机组状态变化情况，同时监测两种劣化指标的变化情况又增加了健康评价的复杂度。因此，为了更全面地描述水电机组健康状态变化情况，同时简化评价方式，增强实用性，将时域劣化指标和频域劣化指标融合得到综合劣化指标CI，如下式所示：

式中：ω1为时域劣化指标TI的权重；ω2为频域劣化指标FI的权重。

为提高CI对异常状态的敏感性，本文取

4 实例分析

4.1 数据说明以某电站3 号机组为研究对象，该机组水轮机型号为ZZA315-LJ-800，额定转速107.1 r/min ，额定功率200 MW，额定水头47 m。2015年8月28日工作人员发现3 号机组运行时有明显异常声音，后停机检查发现该机组转轮室中环钢板出现脱落，中环、下环出现严重裂纹，桨叶裙边损伤严重。从电站状态监测系统中获取了3 号机组2015年8月份历史监测数据355 组，每组数据包含轴向振动A 波形数据和两个工况参数（有功功率和水头）。每个波形包含16 个键相，共4096 个点，采样频率为458 Hz。这355 组数据的工况参数变化范围如表1所示。

表1 3 号机组工况参数变化范围

由于水头变化范围ΔH和有功变化范围ΔP均在3%额定值以内，可将这355 组数据的工况视为同一工况。

4.2 试验设计与结果分析在前235 组数据的采样时间内，现场未发现机组异常，可将前235 组数据视为正常数据。选取前50 组数据的轴向振动波形构成标准集，用来计算标准时、频域特征，第51 组至第150 组数据的轴向振动波形构成正常集，用来计算时、频域差异特征正常阈值。第151 组至第355 组数据的轴向振动波形构成测试集，用来检验所提劣化指标的有效性。

计算时域劣化指标时，首先根据类内距离评估指数得到特征空间重构参数p和q的最佳取值，分别为64、64，然后对各数据集中的样本进行特征空间重构和SVD。计算频域劣化指标时，对轴向振动波形进行快速傅里叶变换得到频谱数据长度为2048，基于Keras 框架下设计自编码器结构如表2所示。标准集中70%的数据用于训练，余下30%的数据用于验证，并采用自适应动量优化算法训练自编码器网络参数来最小化损失函数。时域劣化指标和频域劣化指标计算流程分别如图3和图4所示。

表2 自编码器结构参数

图3 时域劣化指标计算流程

由于本文试验中采用的样本数量远少于实际监测系统中的海量数据，故计算的正常时、频域差异特征值涵盖范围有限，为降低正常状态被误判为故障的概率，在计算时域、频域差异特征正常阈值时，置信区间应取较高的置信度，本文取置信度为99%。

图5至图7为测试集中样本时域劣化指标，频域劣化指标和综合劣化指标。可以看出，3 个指标都能反映机组机组健康状态的劣化过程，起到故障预警作用。从20150822T04∶38∶00 开始，综合劣化指标超过了阈值1，表明机组健康状态出现异常。在20150826T09∶32∶00，综合劣化指标突然增大，超过了10，表明此时机组健康状态明显异常。综合劣化指标超限时间与现场工作人员发现故障的时间相比，提前了6 天，说明所提出的综合劣化指标对故障敏感性较强。此外，在计算机配置为AMD Ryzen 7 4800H+Radeon Graphics 2.90 GHz 的情况下，计算单个样本的综合劣化指标耗时0.014 s，可满足健康评价实时性要求。

4.3 对比试验为了证明综合劣化指标的优越性，针对测试集中的样本，本文计算了5 种常用时域统计指标：平均值，标准差，脉冲指标，歪度指标，峭度指标，计算公式如表3所示，记样本X=[x1，x2，…，xN]，xi表示X第i个点的值，N为X的总点数。计算结果如图8所示。

由图8可得出以下结论：（1）从指标变化趋势来看，在20150822T04∶38∶00 到20150826T09∶32∶00采样时间内，标准差，脉冲指标和峭度指标有增大趋势，与综合劣化指标变化趋势相同，但变化程度不如综合劣化指标明显，而平均值和偏度指标无明显变化，不能体现出机组健康状态的劣化趋势。（2）从指标突变时间来看，20150824T00∶46∶00，脉冲指标和偏度指标幅值发生突变，20150824T05∶45∶00 及20150825T10∶00∶00，平均值发生突变，20150826T09∶32∶00，5 个时域统计指标均发生突变，其中平均值和峭度指标减小，其他3 个指标增大，与综合劣化指标突变时间一致，说明用综合劣化指标进行异常状态预警是可靠的。（3）从指标突变程度来看，综合劣化指标的突变程度最大，最能反映劣化严重程度，从而引起工作人员注意。各指标突变程度按照下式计算：

图4 频域劣化指标计算流程

图5 测试集样本时域劣化指标

图6 测试集样本频域劣化指标

图7 测试集样本综合劣化指标

表3 时域统计指标计算公式

式中，INvn为第n个样本的指标值，INvn-1为第n-1 个样本的指标值。综合劣化指标的ΔINv为268.09%，平均值的ΔINv为3.91%，标准差的ΔINv为55.59%， 脉冲指标的ΔINv为26.52%， 偏度指标的ΔINv为87.42%，峭度指标的ΔINv为29.61%。

由以上结论可知，综合劣化指标与常用时域统计指标相比具有显著优势。为进一步说明综合劣化指标的实用性，本文将其与传统频谱法进行对比。对测试集中部分样本进行FFT 得到其频谱如图9所示。

从图9可以看出，故障特征频率为主要是6f0，12f0，18f0及其谐波。文献［20］将（0.18 ～0.2）f0，（1/6 ～1/2）f0，1f0，2f0和3f0（f0为转频）频段分量信号的能量作为特征进行故障识别，由此可知，对本文案例若采用文献［20］的方法则无法准确识别故障。

综上，综合劣化指标能有效评价机组健康状态，比常用时域统计指标更明显地反映出劣化趋势和劣化程度，且比基于频谱法的水电机组故障诊断实用性更强。

图8 测试集中样本时域统计指标

图9 测试集中样本频谱图

5 结论

本文提出了基于无监督特征学习构建水电机组劣化指标的新方法，针对水电机组状态监测系统中积累的海量正常数据，分别采用基于特征空间重构的奇异值分解算法和自编码器重构误差算法提取状态监测信号的时域特征和频域特征，并建立典型工况下机组正常时域特征库和正常频域特征库，根据正常特征的波动情况确定阈值，计算时域，频域劣化指标及综合劣化指标。应用实例表明，所提出的劣化指标均能有效反映机组健康状态。此外，在表现机组健康状态的劣化趋势和劣化程度方面，综合劣化指标比常用时域统计指标效果更好，实用性更强。同时，本文给出了用综合劣化指标实时评价机组健康状态的方法，通过划分典型工况，使劣化指标适用工况范围更广，评价结果更准确。