一种基于BPNN和SVM-PDE的旋转机械变工况预警方法

2021-07-05崔锦淼胡明辉石保虎

测控技术 2021年6期

崔锦淼，胡明辉,2，冯坤,2*，贺雅,2，石保虎

(1.北京化工大学发动机健康监控及网络化教育部重点实验室，北京 100029； 2.北京化工大学高端机械装备健康监控与自愈化北京市重点实验室，北京 100029； 3.中国石化销售股份有限公司华南分公司，广东广州 510180)

旋转机械是工业应用最为广泛的机械设备之一，如电机、离心泵、离心压缩机、燃气轮机和航空发动机等，涉及到国民经济各个发展领域[1]，由于在高温、高压、变负载环境下长时间连续运行，故障多发。为及时发现设备故障以避免安全事故、减小设备维修成本，需实时监测和评估设备运行状态，及时对设备故障预警。

振动监测是目前常用的旋转机械状态监测手段。传统的监测方法是专业技术人员结合相关振动标准和现场经验设定报警阈值实现设备振动故障的预警，一般取所有工况中的最大测量值作为设备振动烈度，以设备可长期运行的基线值上限25%作为报警限值，即固定报警阈值[2]。该方法虽能在一定程度上避免某些恶性故障的发生，但由于判别特征单一且未考虑设备在时变工况下运行导致的振动水平变化，对于燃气轮机、航空发动机等工况变化频繁复杂的旋转机械，容易造成高工况虚警和低工况漏警的问题。

国内外研究学者在设备报警阈值设定方面开展了大量的研究。魏中青等[3]根据旋转机械发生突变故障时故障部件特征频率幅值突变的故障机理，基于高斯分布3σ法则对故障特征频率幅值设定报警阈值，实现了旋转机械突变故障的快速预警。程晶等[4]基于高斯分布置信度为0.95的单侧置信区间分别设定设备振动总能量和1/3倍频程带级能量特征的报警阈值，该阈值可对舰船设备故障进行有效预警。李会鹏等[5]利用Bayes随机加权法估计小样本数据振动有效值的均值和方差，然后基于高斯分布的3σ法则设定有效值阈值，解决了样本数据不足时阈值难设定的问题。然而，上述预警方法中考虑的特征较为单一，难以全面涵盖设备各状态振动信息和实现早期故障预警。因此，Jiang等[6]基于设备正常状态下振动信号峰峰值、有效值等多个参数的波动特性，设定了多个特征报警阈值，结果表明，多特征报警阈值能够有效提高风电齿轮箱故障预警准确率。明廷锋等[7]通过相关分析选取对故障敏感的特征参数并基于高斯分布3σ法则设定报警阈值，应用时依据报警效果自适应调整阈值，相较于固定阈值应用价值更高。上述报警阈值的设定均基于高斯分布3σ法则，但对于大型复杂设备，其故障形式多样，振动特征往往不属于高斯分布且特征分布复杂未知；并且分别对每个特征设定阈值，忽略了各个特征间的交互作用。

核概率密度估计法是常用的非参数估计法，可解决特征分布未知问题，且其中的核函数可对多元参数进行融合，因此该方法可广泛应用于设备预警阈值估计中。刘晗等[8]基于Parzen窗估计多元监测参数的联合概率密度函数以设定报警限，并随数据积累更新调整阈值，实现了风机机组故障预警。Desforges等[9]基于核概率密度估计求解数据概率密度，依据测量数据概率密度与正常数据概率密度偏离度对电机故障有效预警。Zang等[10]提出一种基于多核概率密度估计的多变量报警阈值设定方法，相较于依据3σ法则设定的阈值，有效降低了虚警率和漏警率。但上述概率密度方法中窗函数宽度难以确定，选取过大或过小都易造成较大误差，很难应用于大样本数据。为了解决该问题，张炤等[11]提出了基于支持向量机(Support Vector Machine,SVM)的概率密度估计(Probability Density Estimation,PDE)方法，该方法利用SVM求解样本概率密度，在相同求解精度下解决了窗函数宽度选取的问题，相较于核概率密度估计方法，该方法工程应用价值更高，更适合于估计振动特征分布特性设定振动阈值。

考虑到变工况对设备振动水平的影响，胡雷等[12]利用相关向量机拟合振动特征与工况参数间的函数，然后设定随工况参数变化的自适应阈值，实现了设备不同工况下的故障预警。顾煜炯等[13]手动划分设备运行工况区间区分设备运行工况，然后分别设定各个运行工况的报警阈值，实现了风电齿轮箱在变工况下的故障预警。刘长良等[14]则利用模糊C均值聚类方法识别设备运行工况，然后设定各工况的预警阈值，相较于固定报警阈值降低了虚警率。但上述变工况阈值设定过程中，当设备工况参数个数过多、运行工况区间过多时，上述工况识别方法难以有效拟合工况参数与工况标记间的函数关系。BP神经网络(BP Neural Network,BPNN)具有强大的拟合能力，可拟合任何非线性函数，且网络节点数可依据工况参数个数自适应调整，不受限于工况参数维度与运行工况区间个数，更具应用价值。

因此，针对目前相关预警方法和设备运行工况识别方法存在的问题，提出基于BPNN和SVM-PDE的旋转机械变工况预警方法。首先利用BPNN准确识别设备运行工况；然后提取各工况振动信号多维特征，并基于SVM-PDE设定各工况振动阈值进行预警；最后利用双转子试验台实验数据验证所提出的方法的有效性。

1 方法理论

1.1 BPNN原理

图1 典型三层BPNN结构

BPNN基于BP算法反向更新权重训练网络，原理如下[15]。

① 网络前向计算。

y=f(w·x+b)

(1)

式中,“·”为矩阵乘积;x为网络输入;w为连接权重;b为偏置，可视为输入为1的神经元所对应的权重；f(·)为激励函数；y为网络输出。

② 误差反向计算。

δi=wi+1·δi+1

(2)

式中，δi为第i层误差项；wi+1为第i+1层的权重；δi+1为第i+1层误差项。网络最后一层的误差计算与选取的目标函数相关，如常用的目标函数为最小化均方误差：

(3)

式中，m为训练样本总数；ti为目标值；yi为网络输出；L为网络误差。则网络最后一层误差项δlast为

(4)

③ 权重梯度计算。

(5)

④ 权重更新。

(6)

利用BPNN识别设备运行工况时，输入层为表征设备运行工况的参数，输出层为人为给定的工况标签，通过BPNN拟合输入与输出即工况参数和工况标签间的函数关系，实现工况参数到工况标签的非线性映射，继而实现多工况参数、多工况区间下的设备运行工况自动识别。

1.2 SVM-PDE原理

张炤、温珂伟、徐玉兵等人利用SVM求解线性算子方程以求解样本概率密度，样本的概率密度函数仅与样本的边缘即支持向量有关，相较于核概率密度估计方法计算量更小，且解决了窗函数宽度选取的问题，该方法原理如下[11,16-17]。

利用SVM求解线性算子方程可转换为对样本概率密度进行回归估计，即

(7)

式中，k(·)为核函数;xi为第i个样本;xj为第j个样本;k(xi,xj)为样本xi、xj经核函数投影值;βj为系数向量β对应的第j个系数；m为样本总数；p(xi,β)为样本xi的概率密度。

其中常用的多项式核函数为

(8)

结合式(7)和式(8)可得样本的分布函数为

(9)

(10)

为防止样本中噪声干扰引起过拟合，引入松弛因子ξ(ξ≥0)和不敏感损失常数ε，保证SVM拟合的概率分布F(xi,β)与经验概率分布F(xi)保持一致性，则有：

(11)

因此，利用SVM求解线性算子方程时，SVM相应的二次规划问题可转换为

(12)

式中不敏感损失常数ε可依据经验公式计算，即

(13)

通过求解式(12)中的二次规划问题，可解得系数向量β，结合核函数矩阵K(xi,xj)，代入式(7)即解得各样本的概率密度值。

故障数据特征分布与正常数据特征分布存在一定的差异，即故障数据各维特征与正常数据各维特征间的距离大于各正常数据各维特征间的距离，因此结合式(7)和式(8)，基于正常数据计算的故障数据概率密度值会小于正常数据概率密度值。以正常数据概率密度下边界作为报警限，当待预测样本概率密度值小于概率密度阈值时，即可视为异常，反之则视为正常，基于此可对故障预警。

1.3 特征提取

从设备振动信号时域波形中提取出有效表征设备运行状态的特征参数，是实现设备状态监测、故障预警和诊断的关键[18]。本文共选取了21个时域、频域统计指标，计算公式如表1所示。

表1 振动信号特征计算公式

上述振动信号特征之间难免会存在信息重叠问题，并且特征维度过高，不利于后续处理分析。主成分分析(Principal Component Analysis,PCA)是常用的无监督线性降维方法[19]，其将高维数据向低维空间投影以实现降维，并保留数据的主要信息。所以可利用PCA对提取到的特征参数进行数据降维，提取二次特征，减少数据冗余。

2 提出的方法

提出的基于BPNN和SVM-PDE的旋转机械变工况预警方法原理如图2所示，主要包含以下步骤。

图2 方法原理图

(1) 设备运行工况识别。

确定表征设备运行工况的参数X，并根据设备振动水平随运行工况表征参数的变化，划分设备运行工况区间，标定设备运行各工况Y。将参数X和标定工况Y划分为训练集和测试集。构建BPNN模型，利用训练集迭代训练网络，测试集测试模型工况识别效果，进而实现多工况参数、多工况区间下设备运行工况的快速准确识别。

(2) 振动数据特征提取。

按照表1所示的特征计算公式，分别从振动信号中提取各个特征指标，并构建出高维特征集。利用PCA对特征集降维，得到低维特征矩阵。基于步骤①中标定的工况，划分低维特征矩阵，得到各工况下的低维特征矩阵，实现各工况振动数据的特征提取。

(3) 变工况报警阈值设定。

基于各工况正常振动状态的低维特征矩阵，求解式(12)的二次规划问题，得到正常数据的系数向量β，进而按式(7)和式(8)计算各工况下的正常数据概率密度值，并以概率密度值边界即下限作为报警限。考虑到数据的波动特性，将报警限下调3倍标准差。

(4) 报警状态判别。

基于步骤①的工况识别模型识别待判定数据所处工况，依据式(8)计算待判定数据与该工况正常数据的核函数矩阵，依据步骤③中解得的该工况正常数据系数向量β和式(7)，计算出待判定数据的概率密度值。比较其与该工况振动概率阈值的大小，当该概率密度值大于阈值时，判定设备处于正常状态；反之，判定设备处于报警状态，基于此实现故障预警。

3 实验验证

3.1 双转子试验台

本研究使用的正常数据和故障数据取自北京化工大学诊断与自愈工程研究中心的双转子试验台。试验采集系统如图3所示，试验台主要由电主轴(外)、联轴器、轴承座、1#轴承、2#轴承和轴承座、联轴器、电主轴(内)组成。本实验通过改变与中介轴承内外圈相连主轴的转速模拟旋转机械变工况运行状态，分别采集内外圈转速在0/300 r/min、300/600 r/min、600/1200 r/min、900/1500 r/min、1500/1800 r/min下中介轴承在正常状态、内圈故障、外圈故障、滚动体故障状态下的内、外主轴转速信号和2#轴承座处的加速度信号，用于验证所提出方法的有效性，其中加速度信号的采样频率为20480 Hz。由于本实验未模拟双转子设备随负载变化情况，故仅以内外转轴转速作为工况标定参数划分5个运行工况(如表2所示)，实际应用中还应考虑负载变化情况。

图3 双转子实验采集系统

表2 工况标定表

3.2 固定阈值预警分析

由于本实验主要模拟滚动轴承故障，加速度信号对滚动轴承故障敏感，故采用加速度信号设定报警阈值，实际应用中应结合设备情况选取位移、速度、加速度确定报警阈值[20]。若采用固定阈值预警法，即不考虑设备时变工况影响，各个工况下使用同一阈值进行报警，选取试验台内外主轴转速300/600 r/min无故障运行状态下振动加速度数据有效值平均值的1.25倍作为报警阈值[21]，即报警阈值为2.4654 m/s2。报警效果如图4所示，在低工况(工况1、工况2)下，内圈故障、外圈故障和滚动体故障的振动有效值均低于阈值，会出现漏警；类似的，在高工况(工况3～工况5)下则有严重的虚警现象。显然，固定报警阈值难以满足设备时变工况下故障准确预警的需求。

此外，进一步对各工况分别设定振动加速度有效值报警阈值，即以各工况正常数据振动加速度有效值平均值的1.25倍作为该工况振动阈值。报警效果如图4所示，仍存在故障值小于阈值的漏警现象，说明振动加速度有效值不能全面反映出设备运行状态，仅以振动加速度有效值单一特征设定报警阈值亦无法满足设备预警需求。

图4 振动有效值固定阈值预警效果图

3.3 所提方法预警分析

3.3.1 工况识别

从已人为标定的工况1～工况5的转速信号中各选取2000组数据，并以1∶1的比例进行划分，分别作为BPNN的训练数据和测试数据。

BPNN参数设置如表3所示，为2-10-1的3层简单网络。利用同工况的训练集和测试集训练并测试BPNN模型，工况识别效果如图5所示，BPNN能准确自动识别设备运行工况。由于实验限制，本实验只考虑内外转轴转速，未考虑负载变化，只模拟5个运行工况，但基于BPNN的工况识别方法可通过调节网络节点参数在多工况参数、多工况区间下自动准确识别设备运行工况，故推广性强。

表3 BPNN参数

图5 工况识别效果图

3.3.2 特征提取

选取5个工况下各120组正常状态的振动数据设定振动阈值，根据表1中的计算公式提取特征。其中，本实验信号最大分析频率为10240 Hz，将频谱以2000 Hz为间隔划分为0～2000 Hz、2001～4000 Hz、4001～6000 Hz、6001～8000 Hz、80001～10000 Hz共5个频带，可求得5个频带能量占比特征，共为25维特征。为了降低数据冗余度，利用PCA对提取的25维特征进行降维。按照图5中识别出的工况区间，对各工况特征降维，得到各工况下低维特征矩阵，工况2的低维特征如图6所示(篇幅有限，其他工况低维特征不予展示)，在轴承各状态下，振动特征间差异明显，即提取的特征对轴承状态敏感，可有效反映出双转子试验台运行状态。

图6 工况2降维特征效果图

3.3.3 报警阈值设定

利用上述提取的各工况低维特征，基于SVM-PDE分别计算各工况下正常数据的概率密度值，并设定各工况下的振动概率阈值。首先基于正常数据求解式(12)的二次规划问题，得到正常数据的系数向量β；然后结合式(7)和式(8)计算正常数据概率密度值，最后以正常数据概率密度值边界即下限作为振动报警阈值。考虑到数据的波动特性，将报警限下调3倍标准差，求得的各工况振动阈值如表4所示。

表4 各工况振动阈值表

3.3.4 报警状态判别

基于测试数据，结合正常数据系数向量β及式(7)和式(8)计算测试数据的概率密度值，基于表4中的振动概率阈值和测试数据的概率密度值判别测试数据状态。其中测试数据组成如表5所示，包括各工况下的60组正常数据、180组内圈故障数据、180组外圈故障数据和180组滚动体故障数据。

表5 测试数据列表

各工况下的振动概率阈值与测试数据的概率密度值间的大小关系如图7所示，正常数据概率值均大于阈值，故障数据概率值均小于阈值，按照报警状态的判别规则，待判别数据的概率密度值小于振动概率阈值会被判定为报警状态，反之则为正常状态。从图7可以看出，该方法在各工况下虚警率、漏警率均为0，有效解决了固定阈值造成的虚警、漏警问题，体现出所提出的方法的有效性。