引言

往复式压缩机、柴油机等复杂机械作为流程工业中的关键设备，一旦出现故障而未及时发现，将会给企业带来巨大的经济损失。目前，企业主要采用单特征门限报警法(简称SF方法)对这些机械进行状态监测及故障预警，以确保生产的安全稳定。然而复杂机械由于自身激励源多、传递路径复杂，其振动信号呈现较强的非平稳特性，致使SF方法常出现误报、漏报警现象^[1-2]。针对该问题，一些学者对多特征融合的机械故障预警方法进行了研究。文献[3]基于双树复小波变换从机械振动信号中提取故障特征构建故障特征向量，利用其训练出的反向传播神经网络模型实现对柴油机的故障预警。文献[4]利用最小二乘支持向量机对风力发电机的多种状态特征建模，以实现故障预警。与SF方法相比，上述两种方法不仅考虑了多项状态特征，而且考虑了各特征间内在联系的变化，从而大幅提升了对复杂机械的故障预警效果，但以上方法均需要大量故障样本训练预警模型，而实际故障样本的稀缺限制了其应用^[5-6]。

文献[7]提出一种基于多元状态估计的机械故障预警方法，其以概率密度采样的方式稀疏历史振动数据并构建过程记忆矩阵，通过比较观测时段与历史正常时段下对应过程记忆矩阵的相似度，实现对电厂辅助设备的故障预警。该预警方法基于设备历史运行工况相似性原理，摆脱了对故障样本的依赖；然而对于存在较多激励源的复杂机械，该方法由于无法识别机械内部各激励源所激发激励信号的变化，对机械局部故障不敏感，导致预警时效性较差。

无限学生t混合模型(infinite student’s t-mixture model，iSMM)是近年来提出的一种无监督学习聚类模型，相比于传统聚类模型如高斯混合模型，它不仅能以较好的鲁棒性来拟合数据的统计分布，而且能够自学习数据的类数以避免出现过、欠拟合数据问题^[8-9]。凭借上述优点，iSMM在图像分割、盲信号处理等领域已得到应用^[10-11]。iSMM在拟合机械振动信号的统计分布时，可根据历史振动数据自学习出机械内激励源的数量以及各激励信号服从的统计分布，从而为识别机械内各激励信号的变化提供有效方式。基于上述分析，本文提出一种基于iSMM聚类的机械故障预警方法，该方法利用iSMM拟合机械振动信号的统计分布，计算在观测时段和历史正常时段下对应的统计分布模型间距离，并将该距离与自适应确定的报警阈值作比较，实现故障预警。

1 iSMM定义及其参数估计方法 1.1 无限学生t混合模型

iSMM由无限个学生t分布加权线性组合构成，其概率密度函数如式(1)所示。

$ P({x_n}) = \sum\limits_{j = 1}^\infty {{\pi _j}} t({x_n}|{\mathit{\boldsymbol{\mu }}_j},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_j},{\nu _j}) $

(1)

其中，

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t({x_n}|{\mathit{\boldsymbol{\mu }}_j},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_j},{\nu _j}) = \int_0^\infty N ({x_n}|{\mathit{\boldsymbol{\mu }}_j},{u_n}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_j})\Gamma (\left. {{u_n}} \right|\\ \left. {\frac{{{\nu _j}}}{2},\frac{{{\nu _j}}}{2}} \right){\rm{d}}{u_n} \end{array} $

(2)

式中，π_j为iSMM中第j个子成分的权重，π_j≥0且$\sum\limits_{j=1}^{\infty} \pi_{j}=1$；μ_j为iSMM中第j个子成分的均值矩阵；Λ_j为iSMM中第j个子成分的协方差矩阵；ν_j为iSMM中第j个子成分的自由度；N(·)为高斯分布概率密度函数，Γ(·)为伽马分布概率密度函数^[10]；u_n为一个服从Gamma分布的隐变量。

为使iSMM能够自学习数据的类数，模型参数需满足条件如式(3)~(5)所示。

$ {{\pi _j} = {\pi _j}(\nu ) = {\nu _j}\prod\limits_{i = 1}^{j - 1} {(1 - {\nu _i})} } $

(3)

$ {{\nu _j}\backsim {\rm{B}}(1,\alpha )} $

(4)

$ {{Z_n}\backsim {\rm{Mult}} (\pi (\nu ))} $

(5)

式中，α为正实数；Z_n为隐变量，用于指示x_n所属的iSMM子成分，在利用iSMM拟合数据时，产生的互不重复的Z_n值的数量即数据的类数。

1.2 模型的参数估计

采用变分贝叶斯推断方法近似估计iSMM的参数，实现对数据统计分布的自学习。该参数估计方法求解快速，且求解精度可满足一般的工程需要。变分贝叶斯推断原理如下所述。

假设存在分布p，其内参数Φ={Φ_j}_j=1^m；p的变分近似分布为q，则式(6)恒成立。

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{ln}}(p(X)) = \int q (\varPhi ){\rm{ln}}\left( {\frac{{p(X,\varPhi )}}{{q(\varPhi )}}} \right){\rm{d}}\varPhi - \\ \int q (\varPhi ){\rm{ln}}\left( {\frac{{p(\varPhi |X)}}{{q(\varPhi )}}} \right){\rm{d}}\varPhi = F + {\rm{KL}}(\left. q \right\|p) \end{array} $

(6)

其中，

$ q(\varPhi ) = \prod\limits_{i = 1}^m {{q_i}} ({\varPhi _i}) $

(7)

式中，X为观测样本，X={x_n}_n=1^N；F为ln(p(X))的变分下界；KL(q‖p)为q(Φ)和p(X, Φ)的KL散度。

由式(6)可知，F越大，则KL(q‖p)越小，即q(Φ)和p(X, Φ)越接近。令F最大化，可得到迭代公式如式(8)所示。

$ {\rm{ln}}(q_j^*({\varPhi _j})) = {E_{i \ne j}}[{\rm{ln}}(p(X,\varPhi ))] $

(8)

其中，

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {E_{i \ne j}}[{\rm{ln}}(p(X,\varPhi ))] = \int {{ \cdots _{{\varPhi _{i \ne j}}}}} \cdots \int {\prod\limits_{i \ne j} {{q_i}} } ({\varPhi _i}) \cdot \\ {\rm{ln}}(p(X,\varPhi ))\prod\limits_{i \ne j}^m {\rm{d}} {\varPhi _i} \end{array} $

(9)

式中，q_j^*(Φ_j)表示对q_j(Φ_j)更新。

利用变分贝叶斯推断方法推断iSMM参数时，式(7)可写成

$ q(\varPhi ) = q(Z)q(u)q(\mu ,\varLambda )q(\alpha ) $

(10)

另外，由于iSMM部分参数集元素数量为无限而导致变分贝叶斯推断的计算无法收敛，因此须设定截断参数T，T满足条件：当j>T，π_j(ν)=0且$\sum\limits_{j=1}^{T} \pi_{j}(\nu)=1$。iSMM的参数估计算法流程如下。

1) 选定观测样本X、{μ_j, Λ_j, ν_j, α}的先验分布。

2) 根据式(6)定义F。

3) 根据式(8)求得分布q(Z_n)和q(u_n)并覆盖原q(Z_n)、q(u_n)。

4) 根据式(8)求分布q(μ_j, Λ_j)、q(ν_j)、q(α)并覆盖原q(μ_j, Λ_j)、q(ν_j)、q(α)。

5) 步骤3)和步骤4)不断交替进行。假设当前迭代步数为s, 若F_s－F_s-1 < ε(ε为精度参数，由人为设定)，则收敛完成。输出q(Z)、q(u)、q(μ, Λ)、q(α)。

2 基于iSMM的故障预警方法

机械振动信号与激励源、传递路径间的关系如式(11)所示。

$ Y(t) = \sum\limits_{i = 1}^N {{F_i}} (t){H_i}(t) $

(11)

式中，Y(t)为t时刻测点处的振动信号函数；N为激励源数；F_i(t)为第i个激励源在t时刻激励力函数；H_i(t)为t时刻第i个激励源到测点处的传递函数。

由式(11)知，测点处的振动信号由机械各激励源激发的激励信号叠加形成。当机械运行状态改变，即部分激励力函数、传递函数发生变化，机械振动信号的统计分布也会随之改变。因此，iSMM拟合机械振动信号所得的统计分布模型能够表征机械的运行状态。基于iSMM的机械故障预警方法流程如图 1所示，包括高维特征空间构建、统计分布模型训练、自适应阈值计算及报警3个部分。

2.1 高维特征空间构建

将数采系统在第i次采样获取的振动数据特征提取后按时间顺序构成数组，并定义其为X_i。利用X_i构建高维特征空间F_i的过程见式(12)。

本文选用基于降噪自编码器的方法进行特征提取^[12]。将所提特征数量设置为10；网络结构设定为n_t－1 000－500－10－500－1 000－n_t，其中，n_t表示单组时域信号的维数；神经元激活函数设为tanh函数；学习率设为0.01。

$ {\mathit{\boldsymbol{X}}_i}\mathop \to \limits^{{\rm{高特征空间构建}}} \left[ {\begin{array}{*{20}{c}} {{f_{1,1}}}&{{f_{1,2}}}& \cdots &{{f_{1,n}}}\\ {{f_{2,1}}}&{{f_{2,2}}}& \cdots &{{f_{2,n}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{f_{m,1}}}&{{f_{m,2}}}& \cdots &{{f_{m,n}}} \end{array}} \right] $

(12)

式中，m为特征类型；n为特征数量；f_{j, k}表示第j个特征类型中的第k个特征。

为避免各类特征的幅值差异对模型训练造成影响，对F_i作归一化处理。对F_i内的任一元素f_{j, k}归一化的计算方法见式(13)。

$ f_{j,k}^* = \frac{{{f_{j,k}} - {\rm{min}}({f_j})}}{{{\rm{max}}({f_j}) - {\rm{min}}({f_j})}} $

(13)

式中，f_{j, k}^*表示对f_{j, k}归一化；min(f_j)为F_i内所有第j类别特征中的最小值；max(f_j)为F_i内所有第j类别特征中的最大值。

2.2 统计分布模型训练

由于小样本易扭曲总体分布形态，为避免此问题，训练统计分布模型所用样本内应包含基于数采系统连续采样构建的多个高维特征空间。样本S表示如下。

$ S = [{\mathit{\boldsymbol{F}}_1},{\mathit{\boldsymbol{F}}_2}, \cdots ,{\mathit{\boldsymbol{F}}_i}, \cdots ,{\mathit{\boldsymbol{F}}_v}] $

(14)

式中，v为样本容量。

利用1.2节参数估计算法使iSMM拟合样本S，可得如式(15)所示的机械状态模型m。

$ m = \sum\limits_{i = 1}^k {{\pi _j}} t({\mathit{\boldsymbol{\mu }}_j},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_j},{\nu _j}) $

(15)

式中，k为iSMM根据样本S确定的类数即激励数；π_j表示第j个类的权重。

2.3 自适应阀值计算及报警

基于iSMM的故障预警方法通过计算机械在正常和实时运行条件下统计分布模型间的距离来定量评估机械当前运行状态，因此如何准确度量模型间距离是确保预警效果的关键。

在计算基准模型与正常状态模型、基准模型与实时状态模型这些统计分布模型间的距离时，基于匹配的KL散度近似算法首先对两模型内的子成分进行匹配，然后对每对匹配的子成分计算其KL散度、权重比，最后计算出两模型间的距离。由于统计分布模型内的子成分表征的是机械振动信号激励信号的统计分布，利用该度量方法的计算机制可间接实现对机械内各激励信号统计分布变化的识别，并且该算法具有较高的计算精度、计算效率^[13]，故采用该算法度量模型间的距离，如式(16)所示。

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{KL}}(\left. a \right\|b) = \sum\limits_{i = 1}^k {{\pi _{a,i}}} (V)\left[ {{\rm{KL}}(\left. {{a_i}} \right\|{b_{w(i)}}) + } \right.\\ \left. {{\rm{ln}}\frac{{{\pi _{a,i}}(V)}}{{{\pi _{b,w(i)}}(V)}}} \right] \end{array} $

(16)

其中，

$ w(i) = {\rm{argmin}} ({\rm{KL}}(\left. {{a_i}} \right\|{b_j}) - {\rm{ln}}({\pi _{b,j}}(V))) $

(17)

式中，KL(a_i‖b_j)表示模型a中第i个子成分与模型b中第j个子成分的KL散度，相应算法参考文献[14]。

在机械历史正常时段下训练出基准状态模型m_base和n个正常状态模型m₁, …, m_i, …, m_n，利用上述算法计算m_base与m_i间距离KL_i(m_base‖m_i), i∈[1, n]。由于机械在正常运行时状态通常较稳定，认为各正常状态模型间差异由随机误差造成，并服从高斯分布。由3σ准则可知，在高斯分布中数值分布在(μ-3σ, μ+3σ)外的置信度仅占不到0.3%^[15]，对应数值判定为非随机误差造成。

报警决策规则定为统计基准状态模型m_base与各正常状态模型间的距离KL_i(m_base‖m_i)，i∈[1, n]。计算这些距离的均值μ和标准差σ。考虑到实际中μ-3σ存在小于0的可能，并且模型间距离非负，故将阈值T设定为μ+3σ。利用机械实时振动数据来训练实时状态模型m_real-times，当KL(m_base‖m_real-times)>T时，认为机械状态处于异常，触发报警。

3 试验验证 3.1 案例数据

以某企业生产现场中的往复压缩机为研究对象，利用不同结构往复式压缩机出现的活塞组件磨损、气阀泄露和液击3种故障案例数据对iSMM方法进行测试，所选故障案例数据详情如表 1所示。

以6缸M型结构的往复式压缩机为例，其传感器布置如图 2所示。在机组十字头的中体位置上安装压电式加速度传感器，采样频率设定为10 kHz；采样间隔设定为往复式压缩机的两个运转周期。依据生产现场中SF方法的报警时间点将故障案例数据划分为正常状态数据和异常状态数据两部分，并对这两部分数据按照样本容量v为100进行划分，将划分余下样本容量不足100的样本舍去。在机组正常状态数据对应生成的样本中，按时间顺序选取前10%样本用于计算报警阈值，其中选取第一个样本用于训练基准模型，其余样本用于测试。

3.2 试验结果及分析

试验分别从报警准确率和相对SF方法报警的提前时间两方面来考察iSMM方法的预警效果。此外，为探究iSMM的鲁棒性在故障预警中的作用，试验中引入基于无限高斯混合模型的故障预警方法(简称iGMM方法)。iGMM聚类机制与iSMM相同，其区别在于拟合数据时鲁棒性相对较差。试验时对同一样本分别使用iSMM和iGMM进行拟合，并分别统计预警结果。

为探究不同样本容量v对iSMM方法报警准确率的影响，将27组故障案例全部用于测试，对iSMM方法进行8次试验，试验中样本容量v由20等间隔递增至160组，测试结果见图 3。可以看出，随着样本容量v的增大，iSMM方法报警准确率随之升高，但模型训练时耗随之增加；当样本容量v超过100时，iSMM方法报警准确率受样本容量影响较小。因此，综合考虑预警效果与模型的训练效率，iSMM方法中样本容量v应设定为100。

iSMM方法和iGMM方法的试验结果见表 2。由表 2可知，对于不同结构机组出现的活塞组件磨损、气阀泄露和液击故障，iSMM方法均实现了报警，相比之下，iGMM方法未能对3缸D型机组气阀故障实现报警。

将iSMM方法和iGMM方法的试验结果进一步统计于表 3。由表 3可知，iSMM方法整体具有较高的报警准确率，可达92.59%；而iGMM方法报警准确率为88.89%。上述结果表明iSMM较iGMM性能更优，且iSMM方法能够较准确地识别出往复式压缩机故障。

为进一步评估iSMM方法的预警效果，将iSMM方法、iGMM方法较SF方法提前的报警时长定义为预警时长，分别计算iSMM方法、iGMM方法相对于SF方法对各类型故障的平均预警时长，计算结果见表 4。

由表 4可知，iSMM方法和iGMM方法均实现了比SF方法提前报警，而iSMM方法较iGMM方法可使报警时间进一步提前；在这3类故障中，iSMM方法对活塞组件磨损故障最为敏感，平均预警时长可达78.4 h；而对液击故障预警时长相对最短，平均预警时长为7.56 h。

针对上述结果，本文对故障案例中结构相对复杂的6缸-M型机组进行具体分析。利用iSMM方法和iGMM方法对该机组3种故障的预警详情分别如图 4~6所示。

由图 4和图 5可知，对于6缸-M型机组出现的活塞组件磨损故障和气阀泄漏故障，iSMM方法计算出的状态异常指标随时间的变化趋势稳定；iGMM方法计算出的状态异常指标随时间波动相对较大，且在超过阈值后的一段时间内仍在阈值线上下间波动，导致此方法较难识别出机组的状态变化。分析认为，iGMM方法计算的状态异常指标波动相对较大是由各个样本内的随机离群点导致，表明iSMM方法可有效抵抗随机离群点的干扰，从而保证较好的预警效果。

由图 6可知，当SF方法对6缸-M型机组的液击故障报警时，iSMM方法和iGMM方法计算出的状态异常指标已远超报警阈值。液击故障指液体进入机组气缸后没有在排气过程迅速排除，在被压缩过程中产生瞬间高液压，这种瞬间高液压可在极短时间内对机组的压缩受力件造成损坏。由于该类型故障在发生前征兆不明显，iSMM方法对其故障预警时长较短。分析认为，该机组在SF方法报警时已发生液击故障，而iSMM方法较SF方法可提前7 h预警机组的液击故障。

4 结论

(1) 提出一种基于iSMM的机械故障预警方法，该方法利用iSMM拟合机械振动信号的统计分布，通过比较机械在正常状态与观测状态下的统计分布模型间距离，实现了复杂机械在无故障样本条件下的故障预警。

(2) iSMM方法可有效地对往复式压缩机进行故障预警，且受样本随机离群点干扰的影响相对较小，具有较高的实际应用价值。

目前仅实现在单一运行工况条件下的机械故障预警，而未考虑时变工况。今后将考虑利用机械的运行工况识别技术，并对iSMM方法在相关方面进行改进。