基于改进LSTM的多工况机电设备剩余寿命预测

葛 阳，马家欣，任 勇，秦健聪

（1.常熟理工学院 机械工程学院，江苏 常熟 215500；2.东南电梯股份有限公司，江苏 苏州 215200）

复杂设备的剩余寿命（remaining useful life, RUL）的高精度预测是当前机电设备可靠性领域的难题[1]．目前常见的寿命预测方法有基于模型的方法和数据驱动方法[2]．基于模型的方法是通过建立数学模型来描述设备剩余寿命与衰退特征的关系，利用观测数据实时更新模型参数，从而预测剩余寿命，例如粒子滤波[3]、马尔可夫方程[4]、指数退化模型[5]、维纳模型[6]等．基于数据驱动的方法则是通过机器学习建立运行数据与退化状态之间的联系，长短期记忆（long short-term memory，LSTM）网络[7]、贝叶斯网络[8]、生成对抗网络[9]等深度学习模型图像识别领域的方法被大量应用于机械设备的寿命预测中．

多工况设备由于不同工况条件下数据分布不尽相同，剩余寿命预测难度大．Cai[10]提出了一种新的基于相似度的剩余使用寿命预测算法和机器预测方法．Li[11]提出了一种利用生成对抗网络预测RUL的方法．Zhu[12]采用隐态马尔可夫模型自动检测状态变化，从而定位故障发生时间．Kundu[13]提出了一种在模型参数估计过程中同时考虑工况参数和状态监测信号的Weibull加速失效时间回归模型．对于工况频繁变化的设备而言，不同工况下传感数据可能存在较大差异，且很多情况下的样本是不均衡的，失效阈值不固定，现有研究对此关注不足．

本文提出了一种健康指标融合方法，用于消除多工况对特征的影响；提出一种自适应权重赋值方法，解决非均衡训练数据带来的过拟合问题；最后利用改进的双向长短期网络进行多工况航空发动机的剩余寿命预测．

1 数据预处理

1.1 工况识别

复杂多变工况往往会掩盖设备的性能退化．不同工况下设备表现一般不同，不利于对其进行剩余寿命预测．本文提出先识别出工况，然后再进行数据的标准化处理．采用K均值聚类算法进行工况识别，具体如下：

设X={X1,X2, …,Xn}为已知数据集，其中X1,X2, …,Xn为n个数据，且每个数据都是N维的，即Xi={xi1,xi2, …,xiN}．K均值聚类算法就是找到含有K个聚类中心的集合C={C1,C2, …,CK}={(c11,c12, …,c1N), (c21,c22, …,c2N), … , (ck1,ck2, …,ckN)使得目标函数

其中：ni是被归为类Ci的数据对象点数；d(Ci,Xj)表示聚类中心与数据对象的欧几里得距离，计算公式为

K均值聚类算法的核心思想是把数据集划分成使目标函数达到最小值的K类，具体步骤如图1所示．首先利用训练数据集通过聚类得到聚类中心和半径，然后用测试数据计算实时工况参数与各个聚类中心的距离，距离较近者即为该工况，识别流程如图2所示．

图1 K均值聚类算法

图2 工况识别流程

1.2 数据标准化

不同工况下设备性能状态参数一般存在很大差异．为了克服工况不同带来的差异，可采用正态标准化方法将不同工况的特征数据映射到一个标定的范围，从而有助于去除差异后体现设备的退化趋势情况，具体如下：

假设第p种工况条件下设备的性能状态参数的均值和标准差s p分别为：

工况p下设备状态参数向量可标准化为

其中i=1,2,…,np，j=1,2,…,N，np表示设备检测点数，N表示性能参数维度．

1.3 特征融合

为了提高预测算法的效率，一般都会设法降低样本数据的维度．其中主成分分析法（PCA）是使用较多的一种降维方法．但是PCA方法难免存在信息丢失问题[12]，为此本文提出了一种特征融合方法，具体如下：

假设第i个设备在k时刻的特征向量为，利用线性加权模型获得第i个设备在k时刻的健康指标

其中 表示权值向量．在具备历史数据的情况下，可以利用最小二乘法求解其值．为了使健康指标HI有明显退化趋势，特将HI的值定义为

其中Ti为第i个设备的寿命．可以看出，在全寿命过程中，设备的健康指标HI值的变化范围是[0,1]．当设备是新品时，HI=0．随着设备服役时间的增加，HI值不断增加，直至单元失效，HI=1．

1.4 时间序列

时间序列是指按照时间先后顺序排列而成的数列，包含设备性能随时间变化的退化信息．对于采用统计方法或神经网络方法进行复杂设备寿命预测而言，时间序列样本是重要前提．为了构造预测时间序列样本，本文采用滑动窗口技术，时间窗口长度取值为40，每次向前滑动1个窗格，从而形成若干训练样本，如图3所示．

图3 滑动窗口示意图

2 深度预测网络

2.1 双向 LSTM（BLSTM）

LSTM是一种特殊的循环神经网络（RNN）类型，对于给定序列x=(x1,x2,…,xn)，利用RNN模型通过迭代计算（8）（9）可以得到一个预测序列y=(y1,y2,…,yn)．

其中：h=(h1,h2,…,hn)为隐藏层序列；Wxh表示输入层到隐藏层的权重系数矩阵；Why表示隐藏层到输出层的权重系数矩阵；bh表示隐藏层的偏置向量；by表示输出层的偏置向量；f表示激活函数，一般是非线性的，如tanh或ReLU函数．

经典的RNN（simple-RNN）很容易出现梯度消失或梯度爆炸的问题．LSTM模型可以学习长期的依赖信息，同时避免梯度消失问题．LSTM在RNN隐藏层的神经节点中增加了一种记忆单元用来记录历史信息，并增加了3种门（Input、Forget、Output）来控制历史信息的使用．

设i、f、c、o分别代表输入门、遗忘门、单元状态、输出门，W和b为对应的权重系数矩阵和偏置向量，σ和tanh分别代表sigmoid和双曲正切激活函数．其迭代计算公式如下：

设备性能的退化对时间存在较强的依赖关系，数据变换不仅与之前时间有关，还反映在后续时间上，BLSTM是在LSTM的基础上优化而来的，具备双向传播能力，可以同时提取前后时间的信息．本文所建立的BLSTM结构如图4所示．

图4 BLSTM结构图

2.2 回归损失函数

（1）均方误差（MSE）

均方误差指的就是模型预测值f(x)与样本真实值y之间距离平方的平均值．其公式为

其中：yi和f(xi)分别表示第i个样本的真实值和预测值；m为样本个数．

（2）平均绝对误差（MAE）

（3）本文提出的损失函数

Huber损失相比于平方损失来说对于异常值不敏感，但它同样保持了可微的特性．它基于绝对误差但在误差很小的时候变成了平方误差．我们可以使用超参数δ来调节这一误差的阈值．当δ趋向于0时它就退化成了MAE，而当δ趋向于无穷时则退化为MSE．其表达式见式（17），是一个连续可微的分段函数．

MSE计算简便，但MAE对异常点有更好的鲁棒性．使用MAE训练神经网络的一个大问题是它的梯度总是很大．这可能导致在梯度下降训练模型结束时丢失最小值．对于MSE，随着损失值逐渐接近最小值，梯度会逐渐减小，从而使其更加准确．在这些情况下，Huber损失函数真的会非常有帮助，因为它围绕的最小值会减小梯度．而且相比MSE，它对异常值更具鲁棒性．因此，它同时具备MSE和MAE这两种损失函数的优点．不过，Huber损失函数也存在一个问题，我们可能需要训练超参数δ，而且这个过程需要不断迭代．鉴于Huber损失函数的优点，本文在其基础上建立回归损失函数

其中：wk为第k个样本的权重；δ=1.35．本文采用KNN（K-NearestNeighbor）方法来确定训练样本的权重，具体步骤如下：

（1）准备数据，对数据进行预处理；

（2）计算测试样本点（也就是待分类点）到其他每个训练样本点的距离；

（3）对每个距离进行排序，然后选择出距离最小的K个样本；

（4）对选出的K个样本进行权重赋值，距离最小的赋值K，依次递减，距离最大的赋值1．

3 剩余寿命预测方法

图5为本文提出的多工况设备剩余寿命预测流程，具体步骤包括：

图5 剩余寿命预测流程

（1）输入训练样本集和测试样本集，利用K均值聚类算法进行训练数据集工况聚类，然后利用欧几里得距离公式计算测试数据到训练集的聚类中心的距离，进而识别出对应的工况．

（2）对于训练样本，针对每一种工况，利用1.2节中的数据标准化方法进行训练数据集的标准化处理，降低工况对特征的影响．测试样本可能包含截断数据，可以将其融合到训练数据中进行标准化处理．

（3）利用1.3节中的特征融合方法，将多个特征融合成为一个健康指标，训练样本和测试样本分别处理．

（4）设定步长和窗口长度，利用滑动移窗方法生成多个健康指标时间序列，并将每个序列的剩余寿命作为标签，训练样本和测试样本分别处理．

（5）利用线性拟合方法对指标时间序列进行拟合，本文采用3次曲线拟合函数，形成新的健康指标，训练样本和测试样本分别处理．

（6）由于训练样本集可能存在非均衡问题，利用KNN方法，为每个测试样本序列匹配K个训练样本，组成新的训练样本集．

（7）将带标签的新训练样本集和测试样本集输入本文构建的BLSTM网络进行剩余寿命预测．

4 实验验证

本文数据来源于美国国家航空航天局（NASA）提供的C-MAPSS涡轮风扇发动机仿真数据集．该数据集每次采样记录航空发动机的24种特征数据，其中包含21个传感器数据和由不同飞行高度、马赫数以及油门杆解算角度组合的运行状态参数．数据集包含4个子数据集（FD001～FD004），每个子数据集均有一个测试集和训练集．训练集数据记录了发动机从正常运行到失效的完整过程数据，测试集数据为发动机失效前的若干循环数据．本文以FD002数据集来验证所提方法．FD002数据集有6种工况，1种故障模式．训练样本数量为260个，测试样本数量为259个．

4.1 工况识别

从FD002训练数据集中挑出飞行高度、马赫数以及油门杆解算角度3个运行状态参数，利用K均值聚类算法识别出6种工况，其中10号和72号发动机的全寿命周期内的工况识别结果如图6所示．

图6 10号和72号发动机寿命周期内工况变化情况

由图6可以看出10号发动机的工况在整个寿命周期内一直处于变动状态．两台发动机的工况变化各不相同，且没有规律性．当然，其他发动机的工况也都是随机变化的．

4.2 故障特征分析

为了比较正常状态和故障出现前发动机运行状态的区别，采用t-SNE二维可视化传感数据，结果如图7所示．其中正常数据取train_FD002数据260个训练样本的最开始的10个周期样本，故障数据取故障发生前的10个周期数据样本．由图7可以看出，正常数据和故障数据的特征聚集性较好、区分明显．这说明发动机故障状态和正常状态的特征区别明显，故可以利用故障时的特征进行发动机的寿命预测．同时也能看出故障特征聚集区不止一处，这说明利用单一阈值预测剩余寿命不合理．

图7 发动机正常状态和故障状态特征t-SNE降维结果

4.3 健康指标

利用1.2节的数据标准化方法，分别进行21个传感数据的标准化．按照1.3节的数据融合方法，将260组训练样本和259组测试样本的21个标准化后的特征融合成1个健康指标，结果如图8所示．可以看出，健康指标从正常状态开始呈逐渐增大趋势，故障时达到最大值，整体上可以体现发动机的退化趋势．

图8 训练样本和测试样本全寿命特征融合结果

图8所示健康指标整体上可以呈现发动机的退化趋势，但从相邻观测时间可以看出健康指标的波动较大，会给寿命预测带来影响．为此，本文采用3次曲线拟合方法对融合后的健康指标进行平滑处理，结果如图9所示．

图9 健康指标的3次曲线拟合结果

4.4 训练样本和测试样本构建

针对训练样本，利用1.4节中的滑动窗口方法，窗口长度设定为40，构建健康指标时间序列样本，共有43 619个样本．以时间序列最后一点对应的剩余寿命作为标签，统计各标签对应的样本数量，如图10所示．可以看出RUL小的样本数量较多，而RUL大的样本数量较少．训练样本不均衡会导致训练过拟合问题．为了克服这个问题，引入了2.2节中提出的利用KNN确定训练样本权重的方法．

图10 样本数量统计

对于测试样本，仅取每个发动机最后40个健康指标作为一个序列样本用来预测RUL．对于时间数量不足40的测试样本可直接忽略．最后共组成测试样本247个．

4.5 评价指标

为了方便比较本文提出的RUL预测方法与现有的一些方法的性能，本文采用平均绝对误差、均方误差以及均方根误差进行度量．

其中：PT表示RUL预测值；RT表示RUL真实值；n表示测试样本个数．以上3个评价指标值越小说明预测精度越高．

4.6 预测结果

实验采用两种方式进行，一是直接利用训练样本集的21个原始传感数据作为特征向量，以对应的RUL作为标签构成训练数据集，以测试样本最后一个观测周期的21个传感数据作为测试数据集．分别利用多层LSTM和多层BLSTM网络，在Kears环境中搭建如表1所示网络结构，配合MSE、MAE以及Huber 3种损失函数进行RUL的预测，预测结果见表2．可以看出，在预测网络结构相同的情况下，使用Huber损失函数的精度最好．相同损失函数条件下，BLSTM网络预测精度不如LSTM．这是因为21个传感数据之间不存在前后因果关系，双向传播网络反而会造成过拟合问题，效果不佳．

表1 预测网络结构

表2 直接利用21个特征预测剩余寿命

二是利用本文提出的健康指标时间序列进行RUL预测，预测网络结构与表1中方法的网络结构相同，RUL预测结果如表3所示．需要说明的是表2中第1～6种方法没有进行训练样本的权重赋值，而第7种方法采用了KNN方法为训练样本赋予权重．在表3的7种方法中，本文提出的预测方法精度最好．同时也能看出，相同损失函数情况下，BLSTM的预测精度好于LSTM．这是因为输入的健康指标时间序列，前后存在一定的退化关系，所以双向移动网络的优势就能体现出来．最后，将测试数据集的RUL预测结果按照真实值大小排序后的对照图见图11．

表3 基于健康指标融合的预测剩余寿命

图11 测试集发动机RUL预测结果

5 结论

本文提出了一种多工况设备剩余寿命预测方法：首先给出了多工况条件下的传感数据的标准化方法，降低工况变化对传感数据的波动影响；其次提出了一种健康指标融合方法，将多个传感特征融合成为一个健康指标，该指标可以体现设备全寿命周期内的劣化趋势；再提出一种基于KNN的训练样本权重赋值方法，以降低因为训练样本不均衡导致的过拟合影响；最后搭建多层BLSTM网络预测测试样本的RUL，并在训练网络中添加了改进的Huber损失函数，以C-MAPSS仿真数据验证本文提出的方法．实验结果表明本文提出的RUL预测方法针对多工况设备具有很好的预测效果．

当然，本文提出的RUL预测方法中，还存在一些有待进一步优化的参数，比如滑动窗口的长度、Huber损失函数中的参数．这些参数对预测精度有直接的影响．如何确定这些参数的最优值以及参数最优值自适应方法等还有待于进一步研究．另外，多故障模式也是设备在实际工作中常见的问题，这类设备的RUL预测也是下一步的研究内容．