CN106153179A - 中速磨煤机振动故障诊断方法 - Google Patents

Abstract

本发明公开了一种中速磨煤机振动故障诊断方法，涉及磨煤机的故障诊断方法技术领域。所述方法首先对磨煤机振动进行故障特征分析，得出磨煤机电流是磨煤机振动故障分析的重要监测参数，随后通过历史数据的相关性从分析筛选出了磨出口风压、给煤机给煤量、磨入口一次风量和磨出口温度四个辅助变量，以这四个辅助变量来预测电流，并用实测电流与预测电流作差，取磨煤机电流残差序列构造磨振动量的方法，对磨电流残差序列进行三层小波包分解，获得8个频带的能量比，统计分析两种故障的能量比，得出磨辊磨损严重和磨内进异物两种故障的特征量区分明显，所述诊断方法具有诊断准确度高、测试效果好的特点。

Description

中速磨煤机振动故障诊断方法

技术领域

本发明涉及磨煤机的故障诊断方法技术领域，尤其涉及一种中速磨煤机振动故障诊断方法。

背景技术

磨煤机振动大是对机组安全稳定运行危害非常大的一类故障，振动大故障的发生会引发一系列的问题：1)磨煤机本体各连接部件松动；2)磨煤机基础振动大可能引发汽泵振动加大；3)拉杆密封损坏，发生漏粉甚至是导致拉杆断裂。而导致磨煤机振动大故障的发生也包含很多可能的原因：1)磨煤机断煤；2)磨盘衬瓦断裂；3)磨内进异物(三块居多)；4)磨辊磨损严重等。而经过调查发现，目前中速磨煤机均未安装磨煤机振动测点，所以对磨煤机振动大故障进行诊断则需要从其他相关参数入手。

发明内容

本发明所要解决的技术问题是提供一种中速磨煤机振动故障诊断方法，所述诊断方法具有诊断准确度高、测试效果好的特点。

为解决上述技术问题，本发明所采取的技术方案是：一种中速磨煤机振动故障诊断方法，其特征在于所述方法包括如下步骤：

(1)建立磨煤机电流预测模型，得出磨煤机电流预测值，将磨煤机电流实际值与磨煤机电流预测值做差，得到磨煤机电流残差序列；

(2)对电流残差序列进行小波包故障特征提取，对提取出的故障特征采用最小二乘支持向量机分类器进行故障分类。

进一步的技术方案在于：步骤(1)中以经过数据预处理的磨煤机进出口差压、给煤机给煤量、磨入口一次风量和磨出口温度作为磨煤机电流预测模型的输入，以BP神经网络预测得到一个电流预测信号，以磨煤机电流实际值与磨煤机电流预测值做差，以得到的磨煤机电流残差序列信号表征磨煤机振动。

进一步的技术方案在于：对构造的磨煤机振动信号，采用三层小波包分解，提取磨煤机振动8个频带的能量比为特征量，通过最小二乘支持向量机分类器法进行故障的诊断。

进一步的技术方案在于：电流预测模型的建立过程如下：

步骤(1)BP神经网络初始化

BP神经网络结构为四输入一输出，八个隐层节点，即输入节点n＝4，隐层节点l＝8，输出节点m＝1；X1为磨出口风压；X2为给煤机给煤量；X3为磨入口一次风量；X4为磨出口温度；Y为磨煤机预测电流；

输入层、隐层和输出层之间的连接权值分别为ω_ij，ω_jk，初始化隐含层阈值为a，输出层阈值为b，设定学习速率为η；

步骤1)网络初始化

inputWeights＝net.IW{1,1}；％输入层与隐含层之间的连接权值

inputbias＝net.b{1}；％隐含层阈值

layerWeights＝net.LW{2,1}；％隐含层与输出层之间的连接权值

layerbias＝net.b{2}；％输出层阈值

步骤2)设置训练参数

net.trainParam.show＝50；训练显示间隔

net.trainParam.lr＝0.01；学习步长

net.trainParam.mc＝0.9；动量项系数

net.trainParam.epochs＝3000；最大训练次数

net.trainParam.goal＝1e-5；最小均方误差

步骤(2)隐含层输出计算

根据输入向量X，输入层与隐含层之间的连接权值ω_ij以及隐含层阈值a，计算隐含层输出H；

$H_j=f(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{ij}}{x}_i-a_j),j=\mathrm{1,2},...,l$

式中，n输入节点数(n＝4)，l为隐含层节点数(l＝8)，f为隐含层激励函数，本文选用tansig激励函数；

$f\left(x\right)=\frac{2}{1+e^{-2x}}-1;$

步骤(3)输出层计算

根据隐含层输出H，连接权值ω_jk和阈值b，计算神经网络预测输出O

$O_k=\underset{j=1}{\overset{l}{\mathrm{\Σ}}}{H}_j{\mathrm{\ω}}_{\mathrm{jk}}-b_k,k=1;$

步骤(4)误差计算

根据网络预测输出O和期望输出Y，计算网络预测误差e

e_k＝O_k-Y_k k＝1

步骤(5)权值更新

根据网络预测误差e更新网络连接权值ω_ij，ω_jk

${\mathrm{\ω}}_{\mathrm{ij}}={\mathrm{\ω}}_{\mathrm{ij}}+\mathrm{\η}{H}_j(1-H_j)x\left(i\right)\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{jk}}{e}_k,i=\mathrm{1,2},...,4;j=\mathrm{1,2},...,8$

ω_jk＝ω_jk+ηH_je_k j＝1,2,…,8；k＝1

式中，η为学习速率；

步骤(6)阈值更新

根据网络预测误差e更新网络节点阈值a，b

$a_j=a_j+\mathrm{\η}{H}_j(1-H_j)\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{jk}}{e}_k,j=\mathrm{1,2},...,8$

b_k＝b_k+e_k k＝1

步骤(7)判断算法迭代次数是否达到最大训练次数或是否达到最小均方差，若未达到，则返回步骤(2)，若达到，则训练结束；

步骤(8)以磨出口风压X1、给煤机给煤量X2、磨入口一次风量X3、磨出口温度X4作为训练好的BP神经网络的输入信号，预测正常的磨煤机预测电流信号Y。

采用上述技术方案所产生的有益效果在于：所述方法首先对磨煤机振动进行故障特征分析，得出磨煤机电流是磨煤机振动故障分析的重要监测参数，随后通过历史数据的相关性从分析筛选出了磨出口风压、给煤机给煤量、磨入口一次风量和磨出口温度四个辅助变量，以这四个辅助变量来预测电流，并用实测电流与预测电流作差，取磨煤机电流残差序列构造磨振动量的方法，对磨电流残差序列进行三层小波包分解，获得8个频带的能量比，统计分析两种故障的能量比，得出磨辊磨损严重和磨内进异物两种故障的特征量区分明显，即将八个频带的能量比作为磨煤机振动大故障诊断的特征量，采用最小二乘向量机对样本进行分类训练、测试，实例验证测试效果较好，表明该方法能有效分离导致的磨煤机振动大的两种原因。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是信号s(k)三层小波分解过程图；

图2是磨煤机振动构造模型图；

图3是采用三层分解的小波包分解树结构图；

图4是线性可分情况下二维输入空间中最优超平面的几何结构图；

图5是磨煤机振动大故障诊断模型图；

图6是给煤机给煤量原始数据与预处理后数据时域图；

图7是预测磨煤机电流与实测磨煤机电流对照图；

图8是磨煤机振动大样本1磨实际电流值与磨预测电流及其差值时域图；

图9是振动大样本2磨实际电流值与磨预测电流及其差值时域图；

图10正常情况下电流残差序列各频带能量比；

图11是磨辊磨损大导致振动大情况下电流残差序列各频带能量比；

图12进异物导致振动大情况下电流残差序列各频带能量比；

图13是本发明中使用的BP神经网络的结构示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

如图2和5所示，总体的，本发明公开了一种中速磨煤机振动故障诊断方法，所述方法包括如下步骤：(1)建立磨煤机电流预测模型，得出磨煤机电流预测值，将磨煤机电流实际值与磨煤机电流预测值做差，得到磨煤机电流残差序列；(2)对电流残差序列进行小波包故障特征提取，对提取出的故障特征采用最小二乘支持向量机分类器进行故障分类。

步骤(1)中以经过数据预处理的磨煤机进出口差压、给煤机给煤量、磨入口一次风量和磨出口温度作为磨煤机电流预测模型的输入，以BP神经网络预测得到一个电流预测信号，以磨煤机电流实际值与磨煤机电流预测值做差，以得到的磨煤机电流残差序列信号表征磨煤机振动，对构造的磨煤机振动信号，采用三层小波包分解，提取磨煤机振动8个频带的能量比为特征量，通过最小二乘支持向量机分类器法进行故障的诊断。

电流预测模型的建立过程如下：步骤(1)BP神经网络初始化

BP神经网络结构为四输入一输出，八个隐层节点，即输入节点n＝4，隐层节点l＝8，输出节点m＝1；如图13所示，X1为磨出口风压；X2为给煤机给煤量；X3为磨入口一次风量；X4为磨出口温度；Y为磨煤机预测电流；

输入层、隐层和输出层之间的连接权值分别为ω_ij，ω_jk，初始化隐含层阈值为a，输出层阈值为b，设定学习速率为η；

步骤1)网络初始化

inputWeights＝net.IW{1,1}；％输入层与隐含层之间的连接权值

inputbias＝net.b{1}；％隐含层阈值

layerWeights＝net.LW{2,1}；％隐含层与输出层之间的连接权值

layerbias＝net.b{2}；％输出层阈值

步骤2)设置训练参数

net.trainParam.show＝50；训练显示间隔

net.trainParam.lr＝0.01；学习步长

net.trainParam.mc＝0.9；动量项系数

net.trainParam.epochs＝3000；最大训练次数

net.trainParam.goal＝1e-5；最小均方误差

步骤(2)隐含层输出计算

根据输入向量X，输入层与隐含层之间的连接权值ω_ij以及隐含层阈值a，计算隐含层输出H；

$H_i=f(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{ij}}{x}_i-a_i),j=\mathrm{1,2},...,8$

式中，l为隐含层节点数，f为隐含层激励函数，本文选用tansig激励函数；

$f\left(x\right)=\frac{2}{1+e^{-2x}}-1;$

步骤(3)输出层计算

根据隐含层输出H，连接权值ω_jk和阈值b，计算神经网络预测输出O

$O_k=\underset{j=1}{\overset{l}{\mathrm{\Σ}}}{H}_j{\mathrm{\ω}}_{\mathrm{jk}}-b_k,k=1;$

步骤(4)误差计算

根据网络预测输出O和期望输出Y，计算网络预测误差e

e_k＝O_k-Y_k k＝1

步骤(5)权值更新

根据网络预测误差e更新网络连接权值ω_ij，ω_jk

${\mathrm{\ω}}_{\mathrm{ij}}={\mathrm{\ω}}_{\mathrm{ij}}+\mathrm{\η}{H}_j(1-H_j)x\left(i\right)\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{jk}}{e}_k,i=\mathrm{1,2},...,4;j=\mathrm{1,2},...,8$

ω_jk＝ω_jk+ηH_je_k j＝1,2,…,8；k＝1

式中，η为学习速率；

步骤(6)阈值更新

根据网络预测误差e更新网络节点阈值a，b

$a_j=a_j+\mathrm{\η}{H}_j(1-H_j)\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{jk}}{e}_k,j=\mathrm{1,2},...,8$

b_k＝b_k+e_k k＝1

步骤(7)判断算法迭代次数是否达到最大训练次数或是否达到最小均方差，若未达到，则返回步骤(2)，若达到，则训练结束；

步骤(8)以磨出口风压X1、给煤机给煤量X2、磨入口一次风量X3、磨出口温度X4作为训练好的BP神经网络的输入信号，预测正常的磨煤机预测电流信号Y。

磨煤机电流预测模型：(1)数据小波降噪处理：一般从热工测量现场采集的数据，由于仪表本身以及环境因素的影响，都不可避免的会影响到数据精度，进而增加建模的难度、降低模型的精确度。基于上述对火电厂热工测量数据的认识，为了保证模型能够较为准确的反映实际情况，本文针对测量数据进行降噪预处理。

含噪声的一维信号s(k)的模型可以表示成如下的形式：

s(k)＝f(k)+ε·e(k),k＝0,1,…,n-1 (1)

式中s(k)—含噪信号；f(k)—有用信号；e(k)—噪声信号；ε—噪声信号偏差。

实际分析中，有价值的信号通常表现为低频信号或是一些比较平稳的信号，而噪声信号则通常表现为高频信号。所以降噪的过程可按如下方法进行处理：对信号s(k)进行小波分解得到低频系数cA₁、cA₂、cA₃和高频系数cD₁、cD₂、cD₃，通过一定方式求得阈值，然后对小波的高频系数进行处理，再对信号进行重构即可达到降噪的目的。抑制信号中噪声部分的同时并恢复出有用信号f(k)，即完成对信号s(k)的降噪，三层小波分解分解过程如图1所示。

一般而言，基于小波分析的一维信号降噪过程可按如下三个步骤进行：

①信号的小波分解。确定小波基和小波分解的层数N，对含噪信号s(k)进行N层小波分解，得到高、低频系数。

②阈值量化小波分解高频系数。对各尺度下的高频系数选择适当的阈值或阈值函数进行阈值量化处理。

③小波重构。依据小波分解系数最高分辨率层的低频系数和经过阈值量化处理的各层高频系数进行小波重构。

(2)数据归一化

归一化是一种数据的无量纲化处理方法，它使物理系统数值的绝对值变成一种相对值的关系，能有效简化计算，缩小量值。对原始数据进行归一化处理，把数据变成(0，1)之间的小数。它有线性函数转化、对数函数转化、反正切函数转化等方式，其中线性函数转化如式(2)所示：

$x_0=\frac{x-x_{\min }}{x_{\max }-x_{\min }}---\left(2\right)$

其中x，x₀分别为归一化前后的值，x_max和x_min为观测的最大值和最小值。

(3)辅助变量筛选

根据运行经验，磨煤机电流的变化能够反映磨煤机振动的状况。故障点前磨煤机电流变化平缓，故障点后磨煤机电流变化剧烈，运行异常分析记录为磨煤机振动大。假若不发生振动大故障，在负荷变化不大的前提下，可以推测磨煤机电流应该像故障点前一样平缓变化，那么实际的磨煤机电流值与正常状态下的电流值之间的残差就是因为振动大故障导致的磨煤机电流变化。

从磨煤机系统的众多参数中进行筛选出磨煤机进出口差压、给煤机给煤量、磨入口一次风量和磨出口温度来表征磨煤机电流，表1、表2为五个参数振动大故障点前12小时正常和故障点后1小时异常运行数据相关性分析。

表1运行数据相关性分析(正常)

表2运行数据相关性分析(异常)

对比表1和表2可得，磨煤机电流与其他四个参数在正常运行情况下呈较高的相关度(>0.408)，而在异常情况下磨煤机电流与其他四个参数呈较低的相关度(<0.126)，磨煤机电流与磨煤机进出口差压、给煤机给煤量、磨入口一次风量、磨出口温度四个参数在正常和异常两种情况下的相关性差异正好能够表征磨煤机电流，故选择该四种参数作为辅助变量来表征磨煤机电流。

磨煤机振动量构造：(1)建立磨煤机电流预测模型：以磨煤机进出口差压、给煤机给煤量、磨入口一次风量和磨出口温度作为磨煤机电流预测模型的输入，构建四输入一输出单隐层的BP神经网络模型。(2)磨煤机振动量构造：将实测的磨煤机电流值与预测电流值作差，以磨煤机电流残差序列表示磨煤机振动量，并对残差序列进行信号分析，实现磨煤机振动故障的诊断。图2中，磨煤机电流残差序列由磨煤机电流实际值减去磨煤机电流预测值得来。

基于小波包分析和最小二乘支持向量机的磨煤机振动故障诊断：对构造的磨煤机振动量(即磨煤机电流残差序列)进行小波包故障特征提取，对提取出的故障特征采用最小二乘支持向量机分类器进行故障分类。

小波包分析理论及其分解、重构算法：小波包分析是小波分析的拓展延伸，是将频带进行多层次划分，对多分频分析无法细分的频带进一步分解，并能够根据被分析信号的特征，自适应地选择合适的频带，使之与信号频谱相匹配，从而提高了时频的分辨率。确定小波包分解的层数需要考虑具体信号和对特征参数的要求，三层分解的小波包分解树结构如图3所示。

在多分辨率分析中，表明多分辨率分析时按照不同的尺度因子j把Hilbert空间L²(R)分解为所有子空间W_j(j∈Z)的正交和。其中W_j为小波函数ψ(t)的小波子空间。对小波子空间W_j按照二进制进行频率的进一步细分，从而实现频率分辨率的提高。将尺度子空间V_j和小波子空间W_j用一个新的子空间统一来表征，由小波分析中可用的分解统一为

$U_{j+1}^0=U_j^0+U_j^1,j\&Element;Z---\left(3\right)$

定义子空间是函数u_n(t)的闭包空间，而是u_2n(t)的闭包空间，令u_n(t)满足式(4)的双尺度方程

$\left\{\begin{array}{c}{u}_{2n}\left(t\right)=\sqrt{2}\underset{k\&Element;Z}{\mathrm{\&Sigma;}}h\left(k\right)u_n(2t-K)\\ u_{2n+1}\left(t\right)=\sqrt{2}\underset{k\&Element;Z}{\mathrm{\&Sigma;}}g\left(k\right)u_n(2t-K)\end{array}\right.---\left(4\right)$

式中，g(k)＝(-1)^kh(1-k)。

由X＝AS构造的序列{u_n(t)(n∈Z₊)}称为基函数u₀(t)＝(t)确定的正交小波包。当n＝0时，u₀(t)和u₁(t)分别为尺度函数φ(t)和小波基函数ψ(t)，因{u_n(t)(n∈Z)}是正尺度函数φ(t)的正交小波包，它构成L²(R)的规范正交基。

{u_n(t)(n∈Z)}是关于h_k的小波包族，设则可表示为

$g_j^n\left(t\right)=\underset{l}{\mathrm{\&Sigma;}}{d}_l^{j,n}{u}_n(2^{j}t-1)---\left(5\right)$

由 $U_{j+1}^0=U_j^0+U_j^1,j\&Element;Z$ 可知，小波包分解就是将分解为与

小波包分解算法是由求与

$\left\{\begin{array}{c}{d}_l^{j,2n}={\mathrm{\&Sigma;h}}_{k-2l}{d}_k^{j+1,n}\\ {d_l^{j,2n+1}=\mathrm{\&Sigma;g}}_{k-2l}{d}_k^{j+1,n}\end{array}\right.---\left(6\right)$

小波包重构算法是由与求

$\left\{d_l^{j+1,n}\right\}=\underset{k}{\mathrm{\&Sigma;}}[h_{k-2l}{d}_k^{j,2n}+g_{k-2l}{d}_k^{j,2n+1}]---\left(7\right)$

小波包分析故障特征提取：由以上分析可以看出，对于振动这一非平稳信号，小波包变换具有不可比拟的优点，它具有多维多分辨率分析的特点，能将任何信号细分解到相应的频带里，用每个频带里信号含有的能量来反映设备的状态，且每个频带能量对于总能量的占比也因工作状态的不同而不同。本文利用利用基于正交滤波器算法的小波包分解故障信号的相关频段，直接利用各频段成分能量比的变化来作为信号的特征量。

以信号的三层小波包分解为例，具体分解步骤如下：

(1)首先对给定待分析信号进行三层小波包分解，分别提取第三层从低频到高频共八个频带的信号特征。

(2)对小波分解系数重构，提取个频带范围的信号，以S_j(j＝0,1,…,7)表示第3层各节点的重构信号，则总信号S可表示为

$S=\underset{j=1}{\overset{7}{\mathrm{\&Sigma;}}}{S}_j---\left(8\right)$

(3)求信号各频带的能量总和。设各频带信号S_j(j＝0,1,…,7)。对应能量为E_j(j＝0,1,…,7)，则有

$E_j=\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{\left|x_{\mathrm{jk}}\right|}^2---\left(9\right)$

式中，x_jk(j＝0,1,…,7；k＝0,1,…,n)表示重构信号S_j的时间域离散点的幅值。

(4)构造特征向量。由于系统出现不同故障时，各频带内信号能量也发生相应变化，因此，以能量比构造故障的特征向量。特征向量T构造如下

T{E₀,E₁,…,E₇} (10)

这些在不同尺度下分解得到不同频带的序列，其相互之间不但是正交的、能量守恒的，而且序列还包含了大量的非平稳故障信息，对非平稳故障可以准确合理的进行故障特征提取。

支持向量机：支持向量机法是一种基于结构风险最小化原理，以统计学习VC维理论为基础的一种学习算法，根据有限样本信息在模型的复杂性和学习能力之间寻求最佳折衷。针对样本数量有限的问题，具备良好的推广能力，有限训练样本求得的解，在进行测试求解时仍能取得良好的效果。

线性可分:在线性可分的情况下，通过构建一个超平面作为决策面来实现分类，要求该决策面不仅能够将所有训练样本进行正确分类，还要使训练样本中离分类面最近的点到分类面距离最大。图4为线性可分情况下二维输入空间中最优超平面的几何结构图，其中：实心点和空心点分别代表两种不同样本；H为分类线；H₁，H₂分别为过各类中离分类线最近的样本点且平行于分类线的直线，它们之间的距离称为分类间隔。

假设给定的样本集为{(x_i,y_i)}i＝1,…,N，其中：x_i∈R^d,y_i∈{1,-1}。x_i是n维向量，如果x_i属于样本一类，则记为y_i＝1，否则记为y_i＝-1。经过归一化处理，分类线H与H₁，H₂的方程总可以表示为

H:wx+b＝0

H₁:wx-b＝1 (11)

H₂:wx-b＝-1

式中w为法向量，b为偏置量。H₁，H₂上的点距分类线H的距离均为D＝1/||w||²，分类间隔为使分类间隔最大相当于是||w||²最小，被误判的风险就最小，因此，求最优分类线可以描述为下面的二次规划问题。

$\min \mathrm{\&Phi;}(w,b)=\frac{1}{2}{\left|\right|w\left|\right|}^2$

s.t.y_i(wx_i+b)-1≥0 (12)

(i＝1,2,…,l)

可以将上述二维最优分类线的推导过程自然扩展到多维情况下最优超平面的推导过程，这时H，H₁，H₂由二维分类线变为多维超平面，超平面H₁，H₂上的点被称为支持向量，这些点唯一地确定了一个最优的分类面H。

利用拉格朗日优化方法可以把上述求解最优分类面问题转化为其对偶问题，即在约束条件对a_i求解下列函数的最大值

$Q\left(a\right)=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{l}{\mathrm{\&Sigma;}}}{a}_i{a}_j{y}_i{y}_j\left(x_i{x}_j\right)---\left(13\right)$

其中：a_i为原问题中与每个约束条件式(13)对应的拉格朗日乘子。这是一个不等式约束下的二次函数寻优问题，存在惟一解。

根据Kuhn-Tucker条件，这个优化问题的解需满足：

a_i(y_i(wx_i+b)-1)＝0(i＝1,2,…,l) (14)

求解上述问题后得到的最优分类函数是：

$f\left(x\right)=\mathrm{sgn}\{\left(\mathrm{wx}\right)+b\}=\mathrm{sgn}\{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{a}_i^{*}\left(x_{i}x\right)+b\}---\left(15\right)$

其中：m表示求和；b是分类阈值；sgn()为符号函数。

由泛函理论可知，高维空间的内积总可以在输入空间找到一个满足条件的核函数K(x,x')，因此并不需要知道非线性映射的具体形式。支持向量机在二维空间能够很好进行分类，但是无法直接求解高维特征空间的点积，而是用原空间的核函数来代替它。上述最优分类函数为

$f\left(x\right)=\mathrm{sgn}\{\left(\mathrm{wx}\right)+b\}=\mathrm{sgn}\{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{a}_i^{*}{y}_{i}K\left(x_{i}x\right)+b\}---\left(16\right)$

根据此函数输出的正负(一般设置为±1)来判别样本所属的类别。根据不同类别的核函数选择来构造不同类型的支持向量机。

线性不可分：在遇到非线性情况时，超平面无法解决分类问题，只能通过提高特征空间的维数来进行分类。首先非线性SVM依据先验的非线性映射条件将待分类向量x_i映射到一个高维特征空间(Hilbert空间)中，然后在此高维空间中构建最优超平面。为了避免在高维特征空间中进行复杂的运算，采用核函数(KernelFunction)在向量之间只作点积运算。其中核函数K须满足：

因此，对于线性不可分，只需要作一个非线性映射将x_i映射到高维特征空间，然后在高维特征空间进行线性分类，这样分类超平面变为：

引入松弛因子ξ_i≥0和错误惩罚常数C>0来解决某些样本不能被正确分类的问题。于是优化问题式(12)重新描述为：

$\min (\frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_i),i=1,...,N---\left(19\right)$

将式(17)转化为拉格朗日函数最小值问题：

式中α_i,γ_i≥0是拉格朗日乘子。式(20)中的极值满足如下条件：

$\frac{\&PartialD;}{\&PartialD;b}L=0\&RightArrow;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_i=0---\left(22\right)$

$\frac{\&PartialD;}{\&PartialD;{\mathrm{\&xi;}}_i}L=0\&RightArrow;C-{\mathrm{\&alpha;}}_i-{\mathrm{\&gamma;}}_i=0---\left(23\right)$

同理，求出最优解w*、b*后，可以得到分类判别函数：

f(x_i)＝sgn{[w*K(x_i,x_j)]+b*} (25)

常用的核函数有：线性核、高斯径向基核和多项式核，它们分别如下

K(x,x_j)＝x·x_j (26)

$K(x,x_j)=\exp \left(\frac{-{\left|\right|x-x_j\left|\right|}^2}{{2\mathrm{\&sigma;}}^2}\right)---\left(27\right)$

K(x,x_j)＝(x·x_j+1)^d,d＝1,…,N (28)

支持向量机多类分类器支持向量机都属于二值分类，但往往实际情况是故障类型可能有很多种，需要进行多类分类。通常利用二值分类方法构造多个二类分类器，通过求解多个分类器而获得多类分类问题的最优解。

给定训练样本集{(x_i,y_i)}i＝1,…,N由多类样本组成，x_i是n维向量，y_i＝1,2,3,…,k表示x_i属于第k类样本。寻找一个决策函数

f(y_i)＝sgn{[(w·x_i)]+b} (29)

支持向量机的多分类算法主要有：“一对多”方法、“一对一”方法、基于二叉树的多类分类方法等，由于本文只需要实现二分类，故不在此详细介绍支持向量机的多分类方法。

基于最小二乘支持向量机分类器：最小二乘支持向量机是基本SVM的一种扩展延伸，采用平方项作为优化指标，并用等式约束条件代替基本SVM的不等式约束条件，即将原二次规划问题转化为线性方程组的求解问题，降低了计算复杂性，比基本SVM更简捷、收敛速度更快，并可获得更高的精度。

样本集{(x_i,y_i)}i＝1,…,N由两类样本组成，x_i是n维向量，如果x_i属于样本一类，则记为y_i＝1，否则记为y_i＝-1。分类超平面

f(x)＝w^TΦ(x)·+b,w^T∈Rⁿ,b∈R (30)

其中Φ(x)为非线性函数，w^T为权向量，b为偏置常量。最小二乘支持向量机的目标最优化函数为：

$\underset{w,b,e}{\min }J(w,b,e)=\frac{1}{2}{w}^{T}w+\mathrm{\&gamma;}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{2}{e}_i^2---\left(31\right)$

s.t:y_i(w^TΦ(x_i)+b)＝1-e

其中e_i为松弛因子，γ为正规化参数。将原空间中的样本依照核函数Φ(x)映射到高维特征空间，使之成为高维特征空间中的一个向量，以解决线性不可分的问题，用拉格朗日法求解如式(32)：

$L(w,b,e,\mathrm{\&alpha;})-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left\{y_i\right[w^T\mathrm{\&Phi;}\left(x_i\right)+b]-1+e_i\}---\left(32\right)$

其中α是拉格朗日乘子。

根据优化条件：

$\frac{\&PartialD;L}{\&PartialD;w}=0\&RightArrow;w=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\mathrm{\&Phi;}\left(x\right)---\left(33\right)$

$\frac{\&PartialD;L}{\&PartialD;b}=0\&RightArrow;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i=0---\left(34\right)$

$\frac{\&PartialD;L}{{\&PartialD;e}_i}=0\&RightArrow;{\mathrm{\&alpha;}}_i={\mathrm{\&gamma;}}_i{e}_i---\left(35\right)$

$\frac{\&PartialD;L}{{\&PartialD;\mathrm{\&alpha;}}_i}=0\&RightArrow;w^T\mathrm{\&Phi;}\left(x\right)\&CenterDot;+b+e_i-y_i=0---\left(36\right)$

可得：

$\left[\begin{array}{cc}0& y^T\\ y& \mathrm{\&Omega;}+{\mathrm{\&gamma;}}^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\&alpha;}\end{array}\right]=\left[\begin{array}{c}0\\ r\\ 1\end{array}\right]---\left(37\right)$

其中y＝[y₁,y₂,…,y_N]^T,Ω＝ZZ^T, $Z={[y_1\mathrm{\&Phi;}\left(x_1\right)\&CenterDot;\&CenterDot;\&CenterDot;y_N\mathrm{\&Phi;}\left(x_N\right)],}_1^r={[1\&CenterDot;\&CenterDot;\&CenterDot;1]}_{N\&times;1}^T,$ α＝[α₁…α_N]。由式(30)可得最小二乘支持向量机分类器的模型：

$f\left(x\right)=\mathrm{sgn}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_i\mathrm{\&Phi;}\left(x_i\right)\mathrm{\&Phi;}\left(x\right)+b)---\left(38\right)$

其中，Φ(x)是核函数，常用的核函数有：线性核、高斯核和多项式核。磨煤机振动大故障诊断模型:基于以上分析，得出磨煤机振动大故障诊断的模型，如图5所示。

磨煤机振动故障诊断实例：

①数据预处理

对预处理采集到的现场原始数据。如图6所示为给煤机给煤量原始数据和预处理后的数据对比，经过数据预处理后，局部范围内降低了信号的噪声干扰，有利于准确建模，其中黑色为原始数据，灰色为预处理后的数据。

②建模

建立给煤机给煤量、磨入口一次风量、磨进出口差压和磨出口风压四个输入量，磨煤机电流一个输出量，单隐层的BP神经网络模型，隐层神经元个数为8，以十小时(36000组)的数据为训练样本，随后一个小时(3600组)数据作为检测样本，建模精度由均方误差MSE反映。将经过预处理的数据导入模型训练检测。图7为模型检测样本的磨煤机电流实际值和预测值之间对照图。图中预测电流与实测电流变化趋势一致。

③导入故障样本

将引起故障的两类原因的12个样本(各6个)导入模型，图8-9所示为两种原因导致的磨煤机振动大的典型电流差值时域图。

图8中，电流差值信号呈无规则的变化，磨辊在磨内经过长时间的运行，特别是在运行时长接近检修周期时，会发生磨辊磨损严重导致磨煤机振动大的故障，磨辊的磨损呈内凹型，每一段的磨损并不规则，与磨盘衬瓦的间隙也不规律，在碾磨的时候产生不规则的振动，电流差值也呈无规则变化趋势，与设备异常分析内容相符。

图9可看出，磨煤机电流差值呈一峰一谷的变化趋势，且峰谷之间时间约为12.5s，此样本的故障为磨煤机进异物导致磨振动大。运行经验可知，运行实际中，磨煤机进异物导致振动大，磨电流呈峰-谷规则变化；从产生机理方面，磨内进异物，如铁块，夹在磨辊与衬板之间，每经过一段时间会与一个磨辊产生较大的摩擦，磨电流亦出现增-降的过程，与此同时产生较为规律的振动。

由上述分析可知，文中提出的电流差值能较为准确的反映出磨煤机振动大的实际故障特征。

④磨煤机振动大特征量提取

截取故障发生后第一段512s的磨煤机电流差值信号数据，对数据进行小波包分析，采用db1小波进行三层小波包分解，生成(3，0)-(3，7)共8个小波包，因为高频部分分解时会发生频带翻转，所以发生4次频带翻转，其频带分布顺序为01326754。以下为正常、磨辊磨损大导致振动大和进异物导致振动大三种情况的各频带能量比。

对比图10、图11与图12可以明显分别出三种不同的状态，表明该方法得到的频带能量比特征明显，可以作为分类的特征量。

⑤建立样本集

定义磨辊磨损大导致的磨煤机振动大为故障类型一，定义进异物导致的磨煤机振动大为故障类型二，统计18个样本的特征量，见表3，表4。

表3故障类型一特征量样本集

表4故障类型二特征量样本集

⑥诊断分类

各取两类样本集中前6个样本为训练样本导入到最小二乘支持向量机进行学习，将剩下3个样本作为检测样本，输出结果为：

表5测试样本测试结果

表中，共有6个测试样本，有5个分类正确，正确分类率达到83.3％，故障类型二的样本8分类错误，其余的均正确分类。详细查看运行故障异常分析得出，该振动大故障检查出两种原因，进异物和磨盘衬瓦断裂，导致该样本分类错误。

所述方法首先对磨煤机振动进行故障特征分析，得出磨煤机电流是磨煤机振动故障分析的重要监测参数，随后通过历史数据的相关性从分析筛选出了磨出口风压、给煤机给煤量、磨入口一次风量和磨出口温度四个辅助变量，以这四个辅助变量来预测电流，并用实测电流与预测电流作差，取磨煤机电流残差序列构造磨振动量的方法，对磨电流残差序列进行三层小波包分解，获得8个频带的能量比，统计分析两种故障的能量比，得出磨辊磨损严重和磨内进异物两种故障的特征量区分明显，即将八个频带的能量比作为磨煤机振动大故障诊断的特征量，采用最小二乘向量机对样本进行分类训练、测试，实例验证测试效果较好，表明该方法能有效分离导致的磨煤机振动大的两种原因。

Claims (4)

1.一种中速磨煤机振动故障诊断方法，其特征在于所述方法包括如下步骤：

(1)建立磨煤机电流预测模型，得出磨煤机电流预测值，将磨煤机电流实际值与磨煤机电流预测值做差，得到磨煤机电流残差序列；

(2)对电流残差序列进行小波包故障特征提取，对提取出的故障特征采用最小二乘支持向量机分类器进行故障分类。

2.根据权利要求1所述的中速磨煤机振动故障诊断方法，其特征在于：步骤(1)中以经过数据预处理的磨煤机进出口差压、给煤机给煤量、磨入口一次风量和磨出口温度作为磨煤机电流预测模型的输入，以BP神经网络预测得到一个电流预测信号，以磨煤机电流实际值与磨煤机电流预测值做差，以得到的磨煤机电流残差序列信号表征磨煤机振动。

3.根据权利要求2所述的中速磨煤机振动故障诊断方法，其特征在于：对构造的磨煤机振动信号，采用三层小波包分解，提取磨煤机振动8个频带的能量比为特征量，通过最小二乘支持向量机分类器法进行故障的诊断。

4.根据权利要求1所述的中速磨煤机振动故障诊断方法，其特征在于电流预测模型的建立过程如下：

步骤(1)BP神经网络初始化

BP神经网络结构为四输入一输出，八个隐层节点，即输入节点n＝4，隐层节点l＝4，输出节点m＝1；X1为磨出口风压；X2为给煤机给煤量；X3为磨入口一次风量；X4为磨出口温度；Y为磨煤机预测电流；

输入层、隐层和输出层之间的连接权值分别为ω_ij，ω_jk，初始化隐含层阈值为a，输出层阈值为b，设定学习速率为η；

步骤1)网络初始化

inputWeights＝net.IW{1，1}；％输入层与隐含层之间的连接权值inputbias＝net.b{1}；％隐含层阈值

layerWeights＝net.LW{2，1}；％隐含层与输出层之间的连接权值layerbias＝net.b{2}；％输出层阈值

步骤2)设置训练参数

net.trainParam.show＝50；训练显示间隔

net.trainParam.lr＝0.01；学习步长

net.trainParam.mc＝0.9；动量项系数

net.trainParam.epochs＝3000；最大训练次数

net.trainParam.goal＝1e-5；最小均方误差

步骤(2)隐含层输出计算

根据输入向量X，输入层与隐含层之间的连接权值ω_ij以及隐含层阈值a，计算隐含层输出H；

$H_j=f(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{ij}}{x}_i-a_j),j=\mathrm{1,2},...,l$

式中，n输入节点数(n＝4)，l为隐含层节点数(l＝8)，f为隐含层激励函数，本文选用tansig激励函数；

$f\left(x\right)=\frac{2}{1+e^{-2x}}-1;$

步骤(3)输出层计算

根据隐含层输出H，连接权值ω_ik和阈值b，计算神经网络预测输出O

$O_k=\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{H}_j{\mathrm{\&omega;}}_{\mathrm{jk}}-b_k,k=1;$

步骤(4)误差计算

根据网络预测输出O和期望输出Y，计算网络预测误差e

e_k＝O_k-Y_k k＝1

步骤(5)权值更新

根据网络预测误差e更新网络连接权值ω_ij，ω_jk

${\mathrm{\&omega;}}_{\mathrm{ij}}={\mathrm{\&omega;}}_{\mathrm{ij}}+\mathrm{\&eta;}{H}_j(1-H_j)x\left(i\right)\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{jk}}{e}_k,i=\mathrm{1,2},...,4;j=\mathrm{1,2},...,8$

ω_jk＝ω_jk+ηH_je_k j＝1，2，…，8；k＝1

式中，η为学习速率；

步骤(6)阈值更新

根据网络预测误差e更新网络节点阈值a，b

$a_j=a_j+\mathrm{\&eta;}{H}_j(1-H_j)\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{jk}}{e}_k,j=\mathrm{1,2},...,8$

b_k＝b_k+e_k k＝1

步骤(7)判断算法迭代次数是否达到最大训练次数或是否达到最小均方差，若未达到，则返回步骤(2)，若达到，则训练结束；

步骤(8)以磨出口风压X1、给煤机给煤量X2、磨入口一次风量X3、磨出口温度X4作为训练好的BP神经网络的输入信号，预测正常的磨煤机预测电流信号Y。