风机是化工厂及化学化工实验室的重要设备，保障着安全与环保，例如排风机能将通风橱及室内的有害气体及时排出。

根据工程经验，风机皮带的寿命在半年左右，皮带断裂会导致风机停机，产生安全隐患，也可能加剧风机轴承的磨损^[1]。皮带断裂不是突然发生的，而是在自然老化过程中产生的小裂缝不断生长导致的^[2]。为此，开发风机皮带寿命预测方法是必要的，用于识别皮带的早期缺陷并预测其寿命终止(end of life, EOL)的时间，以降低实验室的安全风险和人力成本。

本文所研究的风机皮带属于旋转机械的一部分，通行的旋转机械运行监测方法是采集并分析振动信号。这种方法具有成本低、精度高、易实施等优点，也被认为是最有效的旋转机械故障诊断手段^[3]。旋转机械常具有与故障相关的特征频率^[4]，因此这种方法的可用性有物理原理保证。

传统的机械设备故障诊断聚焦于机器是否发生故障的二元划分，但这在工业实践中尚不能满足设备维护的需求。由此便衍生出了能对设备故障提前预报的监测方法^[5]。剩余可用寿命(remaining useful life, RUL)预测在近年来吸引了大量关注，是一种新颖的故障诊断思路。RUL的预测相比传统故障诊断，增加了如下难点：

1) 数据收集过程存在根本困难^[6]。工业设备大多寿命长，且不允许运行至寿命终止，因此数据采集的周期长、也无法获得完整的生命周期数据。为克服这一困难，研究者采用了加速寿命实验^[7]、破坏性实验^[8]等方法。加速寿命实验将设备置于极端温度或极端负荷等非正常工况下运行，使设备加速老化、寿命周期缩短，相当于将健康状态在时间轴上压缩，如图 1中蓝线所示。破坏性实验则使用预先破坏过的设备，实验起点的健康状态与普通实验使用的全新设备不同，相当于将健康状态在时间轴上平移后截断，如图 1中绿线所示。

2) 时间维度的引入增加了模型建立的难度。大部分机械的损坏程度都是随时间非线性变化的^[9]，且损坏程度与体现出的特征的关系通常也是非线性的^[10]。

剩余寿命预测的一般步骤是：数据采集、特征提取、构建HI(health indicator)、划分状态、寿命预测^[6]。这5个步骤中，特征提取与划分阶段有时被省略。

1) 数据采集。

对旋转机械设备剩余寿命的研究多使用加速度计采集振动信号^{[7, 11-15]}。皮带断裂标志着皮带的寿命终点，也是数据采集的终点。

2) 特征提取。

特征是由原始数据经过处理后的量，一个特征表征了原始数据的某一方面。时域特征被广为使用，有信号均方根(root mean square, RMS)^[16]、峰值系数等^[6]。振动信号处理中，频域分析同样重要。Yadav等^[17]使用Sumpeak算法，比较了正常与故障工况下信号的频谱，将差异最大的若干频率分量作为特征。也可将时域和频域分析结合，例如考察信号在特定频段内的能量^[18]。

最近也有使用非线性特征的研究，例如借用混沌理论中的Lyapunov指数来表征机械系统振动状态^[19]。

3) 构建HI。

HI是用于预测寿命的变量。按所使用的特征不同，构建HI的思路主要有3种。

一是省略前述的特征提取步骤，将原始数据作为特征使用，例如Chen等^[20]使用深度卷积自动编码器从长度为2 560的振动信号原始数据中直接构建HI。

二是可以将单个特征或可测量的物理量直接用作HI^{[19, 21]}。

三是使用多特征融合构建HI。例如Khazaee等^[7]将统计学特征作用于单层神经网络，为滚动轴承构建了HI。自动编码器(autoencoder, AE)通过对称的编码层与解码层能实现数据降维，因而适用于从特征构建HI。例如She等^[22]使用正则化的稀疏多层自动编码器，Shen等^[23]使用多层自动编码器(stacked autoencoder, SAE)，将原始特征转化为了HI。

4) 划分阶段。

划分阶段的目的是为了在各阶段采用不同的模型，划分依据是HI的变化率^[6]。在破坏性实验中，由于设备已经进入快速损毁状态，可以不划分阶段。也有其他不划分阶段的研究，例如Chen^[20]认为早期损坏速度慢而晚期速度快，因而将时间的二次函数模型应用至全寿命周期，不划分阶段而用模型自身的特性来体现变化率先慢后快的特点。

5) 寿命预测。

寿命预测阶段将HI序列转化成具有推算未来HI能力的模型，给出HI达到EOL阈值的时刻与当前时刻的差，即得到RUL的预测值。

由上可见，采用破坏性实验的研究相当少，且风机皮带体系在前人的研究中也少有涉及。而破坏性实验能在不改变实验环境与工况的情况下缩短生命周期数据的采集，适合部分工业实际条件。风机皮带寿命较长，且皮带损坏模式单一，人为破坏割裂能近似自然损耗过程中的开裂，因此该体系尤为适合以破坏性实验的手段进行研究。

本文开发了一种适合破坏性实验的HI提取方法，并对清华大学化工系EF-3F-02风机皮带这一实际案例进行了模型验证。

1 理论背景 1.1 信号特征 1.1.1 频域特征

频域特征考察特定频率范围内的信号能量占总能量的比，可以表征风机皮带在频率域中的振动模式。频谱的形状与机器本身的物理结构、材料特性、操作条件等有关。

典型的振动信号频谱如图 2所示，将平滑后的频谱进行多正态分布分解，得到的均值和方差用作滤波器参数，并在同组实验数据中均使用相同参数的滤波器。图 2中，平滑后的频谱经多正态分布分解后有一均值为μ=3 386 Hz、标准差σ=266 Hz、幅值A=0.956的强峰，则基于峰参数换算成带通滤波器上下截止频率分别为μ+σ=3 652 Hz和μ-σ=3 120 Hz。

1.1.2 时域特征

假定有序列x=(x[0], x[1], …, x[n]…, x[N-1])，则一些在故障诊断领域常用的时域特征可如下表示。

1.1.2.1 矩特征

常用的矩特征由1、2、∞阶原点矩M₁、M₂、M_∞及2、3、4阶中心矩M_{2, cen}、M_{3, cen}、M_{4, cen}构造。

1) 方均根(RMS)：

$ y_{\mathrm{RMS}}=M_{2}^{\frac{1}{2}}=\left(\frac{1}{N} \sum\limits_{n=0}^{N-1} x^{2}[n]\right)^{\frac{1}{2}}. $

(1)

y_RMS表征信号的能量，在故障诊断领域广为使用，发生故障时机器振动加剧将导致y_RMS增大。

2) 峰值(peak)与峰值系数(crest factor)：

$ y_{\text {Peak }}=M_{\infty}=\max\limits _{0 \leqslant n \leqslant N-1}|x[n]|, $

(2)

$ y_{\text {Crest }}=\frac{M_{\infty}}{M_{2}^{\frac{1}{2}}}=\frac{y_{\text {Peak }}}{y_{\mathrm{RMS}}}. $

(3)

y_Crest表征信号在时域内的尖锐程度。如果发生故障时产生脉冲冲击，则y_Peak和y_Crest将增大。

3) 偏度(skewness)：

$ y_{\text {Skewness }}=\frac{M_{3, \mathrm{cen}}}{M_{2, \mathrm{cen}}^{\frac{3}{2}}}=\frac{\frac{1}{N} \sum\limits_{n=0}^{N-1}(x[n]-\bar{x})^{3}}{\left(\frac{1}{N} \sum\limits_{n=0}^{N-1}(x[n]-\bar{x})^{2}\right)^{\frac{3}{2}}} . $

(4)

y_Skewness刻画了信号的3阶矩，描述了信号分布的不对称程度。

4) 峰度(kurtosis)：

$ y_{\text {Kurrosis }}=\frac{M_{4, \text { cen }}}{M_{2, \text { cen }}^{2}}=\frac{\frac{1}{N} \sum\limits_{n=0}^{N-1}(x[n]-\bar{x})^{4}}{\left(\frac{1}{N} \sum\limits_{n=0}^{N-1}(x[n]-\bar{x})^{2}\right)^{2}} . $

(5)

y_Kurtosis刻画了信号的4阶矩，描述了信号分布的陡峭程度。

1.1.2.2 非线性特征

由于在发生故障或皮带老化时信号可能会出现更多的随机波动，因此近年来一些学者提出使用动力系统中常用的非线性指标作为描述机器健康状态的特征^[19]。

1) Lyapunov指数(Lyapunov exponent)，记作y_LE。含义为相空间中两条无限接近轨迹的分离率。对离散序列而言，该指数为2条相似样本序列间差异性随时间的指数变化率。以嵌入维度(embedding dimension)等于2为例，若在时刻i和j的2条长度为2的序列间的距离充分近，即

$ D((x[i], x[i+1]), (x[j], x[j+1]))<\varepsilon, $

(6)

其中：ε为事先确定的小量，D(·, ·)为2个序列之间的距离函数。则y_LE可估计如下：

$ \begin{array}{c} y_{\mathrm{LE}}=k^{-1} \log (D((x[i+k], x[i+1+k]), \\ (x[i+k], x[i+1+k]))). \end{array} $

(7)

2) 样本熵(sample entropy)，记作y_SE。含义为相空间中两条接近的轨迹片段在下一时刻的值仍然接近的概率。以嵌入维度(embedding dimension)等于2为例，给定

$ D((x[i], x[i+1]), (x[j], x[j+1])) <\varepsilon, $

(8)

记所有满足式(8)的(i, j)集合为Ω，则y_SE可估计如下：

$ {y_{{\rm{SE}}}} = - \log \frac{{\sum\limits_{(i, j) \in \mathit{\Omega }} I \left( {{D^\prime }(x[i + 2], x[j + 2]) < \varepsilon } \right)}}{{\sum\limits_{(i, j) \in \mathit{\Omega }} I \left( {{D^\prime }(x[i + 2], x[j + 2]) > \varepsilon } \right)}}. $

(9)

其中：I(·)为示性函数，当括号内的条件满足时值为1，否则为0；D′(·, ·) 为2点之间的距离函数。

1.2 自动编码器

自动编码器是一类人工神经网络，由输入层、隐藏层、输出层3层组成，其中输入层与隐藏层构成自动编码器的编码(encoder)部分，隐藏层与输出层构成解码(decoder)部分。特殊的结构使得它适用于数据降维与恢复的功能。

多层自动编码器由多个自动编码器嵌套组成，是一个深度神经网络，具有灵活性和适应性。此外，SAE的每一层可以分别训练。

2 实现方法 2.1 流程图

寿命监测面临数据收集周期长的挑战。因此本文进行破坏性实验，即人为损坏皮带，再将破损皮带上机运行直至寿命终止，此方法也与将皮带生命周期分为“正常运行至损坏开始”与“损坏开始至寿命终止”的两阶段法相似，而本文将关注后一阶段。

本文提出的破坏性实验预测剩余寿命的流程如图 3所示。

2.2 特征提取

根据图 2多正态分布分解的结果，对每组实验数据选取能量最大的5个峰值区段，计算该频段内的信号能量，得到5个频域特征。

同时，将每组实验所得数据序列按1.1.2节中的方法处理，得到6个时域特征y_RMS、y_Skewness、y_Kurtosis、y_Peak、y_SE、y_LE。

2.3 构建HI

HI通常被限定在[0, 1]，且用0表示健康状态，1表示寿命终止。在机器正常的生命周期内，典型的HI大致可分为健康阶段和损坏阶段。健康状态下HI始终保持较低的值；损坏阶段HI出现可追踪的增加趋势，直至机器寿命终止，HI达到1。

图 4中，使用SAE构建指示变量。文[22-23]中将SAE中编码层(图 4中单神经元层，即编码层蓝色方框与解码层红色方框相交的层)的输出直接作为HI使用。但本文提出的方法中，自动编码器将最后一层编码层将线性内积层与激活层的运算分离，图 4中紫色圈中神经元仅进行线性内积运算，输出为前一层神经元输出的线性内积，而绿色圈中神经元仅进行激活函数运算，输出为线性内积被激活函数作用后的结果。进行分离的原因是：未被激活函数作用的紫色神经元输出有更宽的取值范围，保留了更多的信息，更有利于解码层重建特征；同时HI需要神经元输出的取值范围为[0, 1]，而作用Sigmoid激活函数后恰好能满足该要求。其他各层间的连接均使用线性内积与CELU(continuously differentiable exponential linear units)激活函数。

损失函数为

$ L(\boldsymbol{h}, \tilde{\boldsymbol{Y}} \mid \boldsymbol{n}, \boldsymbol{Y})=\|\boldsymbol{h}-f(\boldsymbol{n})\|_{2}+k\|\tilde{\boldsymbol{Y}}-\boldsymbol{Y}\|_{2} . $

(10)

其中h=(h[0], h[1], …, h[n], …, h[N-1])为自动编码器构建出的HI序列，即图 4绿色圈中神经元的输出，f(n)=(f[0], f[1], …, f[n], …, f[N-1])为从参考时间序列n=(0, 1, …, n, …, N-1)构建出的HI预估值的序列，Y为模型的输入值，即各时刻各特征y组成的矩阵(y由2.2节中的11个特征组成)，$ \tilde{\boldsymbol{Y}}$为模型的输出值，即重建出的特征矩阵。可见损失函数由2部分组成，第一部分是HI与其期望形式的差异，第二部分是原始输入特征与经过多层编码与解码后的特征的差异。k为权重因子，调整上述二者的比重^[20]。在式(10)中，最为关键的是f(n)的选取，它是对HI函数形式的预估。在接近寿命终止的阶段，损坏程度逐渐加快，因此f(n)也应随时间非线性增大。本文实验中皮带经过破坏处理，初始时刻的HI取0不合适，考虑到皮带中起到骨架作用的是包裹在粘合橡胶层内的绳芯，因此可以按初始时人为破坏绳芯的比例c确定HI初始值，因此f(0)≠0，同时f(0)是c的函数。例如c=3/5表明实验开始时破坏了5根绳芯中的3根。综上f(n)在n∈[0, N-1]应满足以下条件：

$ f(N-1)=1, $

(11)

$ f(0)=g(c) \in(0, 1), $

(12)

$ f^{\prime}(n)>0, $

(13)

$ f^{\prime \prime}(n)>0. $

(14)

构建如下函数：

$ f(n)=\left(\frac{1-c}{N-1} \cdot n+c\right)^{a}. $

(15)

式(15)满足f(0)=c^a，f(n)=1^a=1，可看做幂函数f₀(n)=n^a平移、缩放、截断后的结果，这也与图 1中破坏性实验和正常实验中健康状态的关系类似。

对于HI序列h，常用的量化评价标准有3个^[6]：单调性(monotonality)、趋势性(trendability)、稳健性(robustness)。值的可能范围均为[0, 1]。

1) 单调性。

$ z_{\mathrm{Mon}}=\frac{1}{N-1}\left|\sum\limits_{n=1}^{N-1}\left(I_{\{h[n]>h[n-1]\}}-I_{\{h[n] <h[n-1]\}}\right)\right|. $

(16)

其中I_{·}为示性函数。

2) 趋势性。它是h与参考时间序列n的Pearson相关系数。

$ \begin{array}{l} {\mathit{z}_{{\mathop{\rm Trd}\nolimits} }} = |{\mathop{\rm Corr}\nolimits} (\mathit{\boldsymbol{h}}, \mathit{\boldsymbol{n}})| = \\ \left| {\frac{{N\sum\limits_{n = 0}^{N - 1} n h[n] - \sum\limits_{n = 0}^{N - 1} n \sum\limits_{n = 0}^{N - 1} h [n]}}{{\sqrt {N\sum\limits_{n = 0}^{N - 1} {{n^2}} - {{\left( {\sum\limits_{n = 0}^{N - 1} n } \right)}^2}} \sqrt {N\sum\limits_{n = 0}^{N - 1} {{h^2}} [n] - {{\left( {\sum\limits_{n = 0}^{N - 1} h [n]} \right)}^2}} }}} \right|. \end{array} $

(17)

由于实际应用中，剩余寿命是随时间线性减小的，理想的指示变量应该随时间线性单调变化，此时z_Mon和z_Trd取到最大值1。

3) 稳健性。

$ z_{\text {Rob }}=\frac{1}{N} \sum\limits_{n=0}^{N-1} \exp \left(-\left|\frac{h[n]-\widetilde{h}[n]}{h[n]}\right|\right). $

(18)

其中$\tilde{\mathit{\boldsymbol{h}}} $为平滑后的h，本文使用LOESS平滑。z_Rob考察HI在总体趋势以外的波动程度，理想的指示变量在平滑前后的差别很小即h[n]-$\tilde{h} $[n]≈0，因此z_Rob将接近1。

实际应用时应同时考虑上述三者，本文构造加权平均的综合指标为

$ z=0.2 z_{\mathrm{Mon}}+0.5 z_{\mathrm{Trd}}+0.3 z_{\mathrm{Rob}} \in[0, 1] . $

(19)

2.4 剩余寿命预测

使用二阶多项式对HI在时间上进行回归。考虑到HI的单调且逐渐加快的变化趋势，多项式的系数均限制为正数。系数均随数据的离散编号n实时变化，在n时刻用历史数据得到关于连续时间t的回归结果h(t)=a₂[n]t²+a₁[n]t+a₀[n]后，n时刻的RUL可以预测为

$ \mathrm{RUL}[n]=\inf\limits _{t}(h(t)>1)-n. $

(20)

3 案例分析 3.1 数据介绍

实验使用的风机为清华大学化工系排风系统机群中编号为EF-3F-02的风机。风机配有两根皮带，皮带型号为Mitsuboshi公司生产的SPZ1600/3V360。振动传感器采样率为f_s=25 600 Hz，由数据采集器和输入模块将传感器数据传入计算机。实验硬件配置如图 5所示。

皮带具有4层结构(见图 6)，即覆盖皮带四周的织物层、粘合橡胶层、绳芯层、压缩橡胶层。引起皮带断裂的关键在绳芯层。

使用刻刀破坏皮带，使压缩橡胶层全部断裂，并破坏粘合橡胶层使绳芯断裂3~4根，待风机启机运行平稳后，按每60 s采集10 s数据的方式采集实验数据，因此每段数据序列长度为256 000。由于数据存在噪声，将原始特征按窗宽120进行平滑。

进行3组实验，采集到的数据如表 1所示。

3.2 预测结果

以式(19)为标准优化超参数可得编码器层数取4、学习率取0.01、损失函数取f[n]=n⁴且k=10^-2、滑动平均窗宽取120时误差最小。

训练模型，得到的HI曲线如图 7所示。模型给出的实验1、2、3的HI的z分别为0.799 2、0.787 2、0.798 8。

图 8是使用每一时刻之前的所有历史HI数据，采用2阶多项式拟合HI得到的预测RUL及真实RUL。将图 8中每张子图的预测和真实RUL作差，得到图 9的预测RUL的绝对误差，可见在3组数据后1/3时间段内的预测中，绝对误差均不超过1 d。图 9中y轴混合使用对数和线性坐标，y轴坐标绝对值为[0, 0.1]时为线性坐标，其他范围y轴坐标绝对值采用对数表示。

3.3 训练集外数据的预测

在实验1、2、3组成的训练集外有第4组实验数据，它包含6.26 d共9 015个数据点。图 10为应用训练得到的模型所得的HI。图 11为用2阶多项式拟合HI得到的预测RUL。

从中可见HI较好地落在了[0, 1]内，且呈现出了良好单调增加趋势，这说明使用训练所得模型提取HI是可靠的。图 11中，实验刚开始时由于数据量较少，预测的RUL曲线仍然波动较大，但后期数据渐渐稳定，且在数据采集开始6.2 d后，皮带剩余寿命预测约为5.4 d。

3.4 指示变量与解码器嵌入程度的比较

使用SAE提取HI的文献[22-23]多将最中间层的输出作为HI，是上一层输出的线性内积后通过激活函数的结果，称为传统方法。而本文将线性内积与激活函数分离，将未通过激活函数的内积结果直接输入解码层，而使用独立的Sigmoid层将内积结果转化为HI，且该HI不参与解码。

为比较两种方法在模型训练时的差异，比较了两种方法的损失函数和3组数据的z：使用表 1中的3组数据进行10次模型训练实验，每一次训练实验得到一个模型的损失函数Loss，并将该模型用于同样的3组数据，对每一组数据可得到与之对应的z，如图 12所示。本文方法的损失函数平均值为0.006 446(见图 12b)，小于传统方法的0.007 505(见图 12a)。同时z的差别很小，这是因为2种方法在由原始特征到HI的编码过程相同。模型训练中，更小的损失函数意味着模型拟合效果更好，并且解码层输出的特征与原特征更接近。

对每次模型训练，通过得到的HI可以计算预测RUL，而后得到实验1、2、3的各时刻的预测误差e(即预测RUL减去真实RUL)与反向累积平均误差$ \tilde{e}$，$ \tilde{e}$的定义为

$ \tilde{e}[n]=\frac{1}{N-n} \sum\limits_{i=n}^{N-1}|e[i]| . $

(21)

它考察了预测值在接近寿命终点时的预测精度。图 13是10组模型训练实验给出的全过程的$ \tilde{e}$的均值，本文方法所得的橙色线在过程的后半段均在蓝色线所代表的传统方法下方。表 2是在皮带寿命后70%、50%、30%时间段内计算得到的$ \tilde{e}$，3组实验后50%时间段的$ \tilde{e}$比传统方法分别减小了0.12、0.18、0.58 d。这都表明本文方法在接近寿命终点时有更低的预测误差。

4 结论

本文基于振动信号测量、破坏性实验、自动编码器实现了风机皮带的寿命预测。在使用自动编码器构建HI中，将线性内积运算和激活函数运算分离，输入解码层的值为线性内积运算的结果，而提取的HI为激活函数层运算的结果。这有助于解码器解码，且Sigmoid函数自有的特性将HI限制在了[0, 1]。实验案例证实了数据融合、将激活函数与线性内积分离计算的优势，3组实验中，在后50%时刻RUL的预测误差比传统方法分别减小了0.12、0.18、0.58 d。

本文仍需要更多数据来获得更好的模型。此外，破坏性实验虽然能缩短皮带的生命周期，但它与皮带在正常运行中的自然损坏存在些许差别。因此，应当尽可能多地收集皮带生命周期的数据，并探索剩余寿命预测方法。