一种水泵机组振动信号特征提取与状态识别方法

技术领域

本发明涉及振动冲击信号处理技术领域，具体是涉及一种基于CEEMDAN样本熵和LCGSA-RBF的水泵机组振动信号特征提取与状态分类识别方法。

背景技术

水泵是输送液体或使液体增压的机械，它将原动力机的机械能或其他外部能量传送给液体，使液体能量增加，主要用来输送液体包括水、油、酸碱液、乳化液、悬乳液和液态金属等，也可输送液体、气体混合物以及含悬浮固体物的液体。水泵机组主要包括水泵、电动机、控制系统等部分，因其在工业系统中扮演的角色，导致机组运行工况变化频繁，会易产生机组振动现象，且振动异常可能导致机械故障、损失和安全事故。现有的关于水泵机组振动的研究越来越多，而研究的第一步就是需要对振动信号进行处理。

对水泵机组振动信号处理的一般步骤是将其进行特征提取以得到其中有用的特征信息，再使用分类器对其状态进行分类识别，以便于后续的故障诊断、预测和控制工作。对于水泵机组振动信号处理，目前有一部分采用传统的特征提取算法和分类器进行处理，时域分析、频域分析、小波变换、奇异值分解、支持向量机等，但这些方法特征提取的特征不够准确，在分类器的训练中容易出现过拟合等问题，导致结果的可靠性和稳定性差；另一部分采用深度学习方法，如卷积神经网络，递归神经网络等，虽然这些方法在特征提取和分类器构建上相对传统方法更具优势，但在样本计算和分类器训练时需要较大的数据量，而目前水泵机组振动信号数据较为有限，因此这种方法在实际应用中仍存在一定的局限性。因此，亟需开发一种适用性更广、可靠性更高的水泵机组振动信号特征提取与状态识别方法。

发明内容

本发明的目的在于提供一种水泵机组振动信号特征提取与状态识别方法，解决现有技术方法中效果不理想、可靠性差且局限性大的问题。

为解决上述的技术问题，本发明采用以下技术方案：一种水泵机组振动信号特征提取与状态识别方法，其特征在于包括如下步骤：

S1. 分别采集水泵机组健康状态和故障状态下的振动信号，将采集到的水泵机组的振动信号采用MODWT小波变换进行降噪；

S2. 对降噪后的状态信号进行CEEMDAN分解，得到 7个IMF分量和1个最终余量R(n)；

S3.通过相关系数法计算每一个IMF分量的相关系数，并选取相关系数最高的4个IMF分量；

S4.计算选取后的IMF分量的样本熵，并整合成特征样本集，并按状态分为健康样本和故障样本；

S5. 将特征样本集按照1：9划分为测试样本和训练样本，采用LCGSA-RBF方法对训练样本进行状态识别，得到训练好的模型；将测试样本输入到训练好的模型中进行识别准确率验证；所述LCGSA-RBF方法具体为将训练样本输入到RBF神经网络对数据进行拟合，得到初始的RBF神经网络模型，再将RBF神经网络中的参数转化为搜索空间中的个体，根据其质量和位置计算出相互之间的作用力和速度，更新其位置和速度，在每次迭代中，利用GSA算法重新计算每个个体的质量和位置，更新搜索空间中的个体，在更新位置和速度时，引入莱维飞行和混沌理论的随机因素，根据搜索空间中的最优个体，得到最优的RBF神经网络模型。

更进一步的技术方案是所述步骤S1对水泵机组振动冲击信号采用最大重叠离散小波变换（MODWT）进行去噪处理。MODWT是DWT的升级版本，

S101.设X为N维向量，元素为实值时间序列

S102. 设尺度 

根据 Mallat 算法，计算出尺度 

其中，mod表示两数相除取余数；其中

S103.DWT中的尺度滤波器

满足：

其中，

S104.在尺度j下的

S105.根据Mallat 算法，计算出尺度j下的尺度变换系数和小波变换系数分别如下：

更进一步的技术方案是所述步骤S2对于降噪后的信号进行CEEMDAN分解，CEEMDAN是基于EEMD所改进，加入经EMD分解后含辅助噪声的 IMF 分量，而不是将高斯白噪声信号直接添加在原始信号中，抑制了EEMD的模态混叠和计算复杂度高以及虚假IMF多的问题。

CEEMDAN分解具体步骤如下：

定义

S201.采用以EEMD相同的方式获得各阶固有模态分量，首先设原始信号

S202.对于每个加入高斯白噪声的原始信号

S203.根据上述公式计算得到第一阶固有模态分量

计算获得第一阶段固有模态分量

S204.向第一阶段分解的余量加入高斯白噪声，并对处理后的信号余量进行EMD分解，对第一阶段余量产生的第一个固有模态分量进行总体平均计算即获得第二阶段的固有模态分量

计算获得第二阶段固有模态分量

S205.依此，对振动信号进行多次EMD分解直至第

得到CEEMDAN的第

S206.直至分解后的余量无法继续分解，完成CEEMDAN分解流程，最终经CEEMDAN分解后的信号为：

式中：

更进一步的技术方案是所述步骤S3中所有IMF分量的相关系数计算，并选取相关系数最好的4个分量，其具体步骤包括：

在本发明的筛选流程中，经 CEEMDAN分解后的固有模态分量可依据各分量的相关系数进行筛选，相关系数的数学表达式为： 

式中：X，Y 代表两个样本变量，r(X,Y)表 示 两 样 本 变 量 之 间 的 相 关系 数 ，

更进一步的技术方案是所述步骤S4中对于相关系数最高的4个分量进行样本熵计算其具体步骤包括：

对于由N个数据组成的时间序列

S401.按序号组成一组维数为m的向量序列，

式中：

S402.定义向量

式中：

S403.对于给定的 

S404.对所有的

S405.增加维数到m+1，计算 

得到序列样本熵的理论值为：

但是在实际中N不可能为无穷大，而为一有限值，则样本熵的估计值为：

S406.根据上式可知，嵌入维数m和相似容限r的大小与样本熵的结果相关；设定

更进一步的技术方案是所述步骤S5中对将样本熵作为样本集输入到LCGSA-RBF算法中进行状态识别与分类，其具体步骤包括：

设置初始参数如下：粒子维度为2，方式为8 ，种群数量为20，迭代次数为15，具体地：

S501.是通过随机分布的大质量个体吸引其他粒子来完成GSA的初始化，两个粒子间的作用力为：

式中：

S502.其具体定义如下：

式中：

S503.当个体在搜索空间中移动以寻找最优解时，其受到其它粒子的合力为：

式中：

S504.该质量的加速度为：

S505.在进化过程中，粒子的速度和位置的更新方式为：

S506.粒子的质量与适应度值有关，质量越大的粒子越接近最优，并且它对其他粒子的作用力相应地会更大，但移动速度较慢，粒子质量的计算方式为：

式中

S507.Levy的飞行公式如下所示：

S508.Levy飞行的步长公式为：

式中：

参数的方差为：

S509.混沌归一化的表达式为：

式中：

式中：

S510.将Levy和混沌序列与GSA的引力常数结合，得到LCGSA的引力常数

S511.RBF使用欧氏距离及高斯函数，其表达式为：

式中：

S512.最终的输出

式中：

样本的期望输出值，其表达式为：

工作原理：本发明提出了基于CEEMDAN样本熵和LCGSA-RBF的水泵机组振动信号特征提取与状态分类识别方法，首先对原始信号进行WODWT小波降噪，再对降噪后的信号进行CEEMDAN分解，得到一系列IMF分量，同时计算各个分量的相关系数，选取相关系数最高的4个分量进行样本熵计算，将计算好的样本熵作为样本集输入到RBF神经网络对数据进行拟合，得到初始的RBF神经网络模型，之后，将RBF神经网络中的参数转化为搜索空间中的个体，并根据其质量和位置计算出相互之间的作用力和速度，更新其位置和速度，在每次迭代中，根据GSA算法的原理，重新计算每个个体的质量和位置，更新搜索空间中的个体，最后在更新位置和速度时，引入莱维飞行和混沌理论的随机因素，增加算法的随机性和全局搜索能力，最终，根据搜索空间中的最优个体，得到最优的RBF神经网络模型，从而实现对数据的优化拟合，使得模型能更加准确地将振动信号所对应的水泵机组状态进行识别。

与现有技术相比，本发明的有益效果是：

（1）相较于传统的原始信号分解方法，采用本发明提出的CEEMDAN分解方法所提取的样本熵更适用于模型的读取与辨识。

（2）基于CEEMDAN样本熵和LCGSA-RBF的水泵机组振动信号特征提取与状态分类识别系统对于旋转机械的声音振动信号的识别优于传统系统，利用噪声传感器可以在和机组非接触的条件下，精准获取分析识别机组的突变信号，可以有效地监测机组的工作状态。

（3）引入莱维飞行和混沌映射方法对GSA进行优化，增加算法的收敛速度和全局搜索能力，同时建立LCGSA-RBF的全新模型应用于水泵机组振动信号的识别应用中。

（4）本发明设计全新模型解决水泵机组故障诊断效率低、反应速度慢等问题，在水泵机组特征提取领域具有很强的优越性和实用价值。

（5）本发明的方法对于促进水泵机组的安全可靠运行具有重要而积极的意义，通过提供自动化的振动信号分类和识别，可以有效地监测机组的工作状态，及时发现异常情况并采取相应的措施，从而提高机组的安全性和可靠性，延长其使用寿命。

附图说明

图1为本发明的流程图。

图2为健康状态时水泵机组的实测信号波形图。

图3为去噪后的健康状态下水泵机组的实测信号波形图。

图4为健康状态时经CEEMDAN分解的IMF分量图。

图5为健康状态时和振动冲击时两组IMF分量相关系数图。

图6为RBF神经网络基本架构图。

图7为振动冲击时水泵机组的实测信号波形图。

图8为去噪后的振动冲击时水泵机组的实测信号波形图。

图9为振动冲击时经CEEMDAN分解的IMF分量图。

图10为状态分类识别准确率图。

图11为本发明与CEEMDAN样本熵和RBF、EMD样本熵和LCGSA-RBF、CEEMD样本熵和LCGSA-RBF、VMD样本熵和LCGSA-RBF四组模型对比图。

图12为5种模型20次运行后的递增平均值图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述，但以下实施例仅是说明性的，本发明的保护范围并不受这些实施例的限制。

该水泵机组额定功率55KW，额定电压400V，额定水头10m，转速为81.6r/min，因实际机组中很难发生故障，本发明采取机组运行状态下对顶盖进行连续敲击模拟故障条件下振动信号的产生，并采用非接触试噪声传感器CRYSOUD进行实时采集，该传感器内置先进的DSP处理器。集传声器、前置放大器和数据采集卡于一体，在结构布置紧密的前提下，又保证了数据处理高校准确。用噪声传感器共采集得到机组2种不同状态下的振动信号，分别为故障和健康两组样本，每种状态的数据集为1000个样本，每个样本的点数为4096。传感器参数如表1所示，数据采集卡的型号为研华IDAQ801。

表1 噪声传感器CRY2301 参数

下面结合实例描述本发明所述的一种基于CEEMDAN样本熵和LCGSA-RBF的水泵机组振动信号特征提取与状态识别方法，针对实例中的水泵机组的振动冲击信号进行状态分类识别，并将结果与CEEMDAN样本熵和RBF、EMD样本熵和LCGSA-RBF、CEEMD样本熵和LCGSA-RBF、VMD样本熵和LCGSA-RBF四种系统处理同一种信号的结果进行对比对比，以验证本发明的优越性，包括如下所述具体过程：

1.水泵机组健康状态时振动信号的特征提取与状态分类识别。

如图1所示，本发明实例提供的基于CEEMDAN样本熵和LCGSA-RBF的水泵机组振动信号特征提取与状态识别方法，包括如下步骤：

 S1. 将采集到的水泵机组原始振动冲击信号采用MODWT小波变换进行降噪；

S101.MODWT是DWT的升级版本，其具体原理如下：

S102.设X为N维向量，它的元素为实值时间序列

S103. 设尺度 

则根据 Mallat 算法, 可以计算出尺度 

其中，mod表示两数相除取余数。

S104.且：

S105.MODWT与DWT中的尺度滤波器

S106.且同样要满足：

或者

S107.为了避免由于抽取采样导致变换后样本点数减少的情况，即在尺度j下的

 S108.根据MALLAT算法，可以计算出尺度j下的尺度变换系数（近似）和小波变换系数（细节）分别如下：

去噪前后的波形图分别如图2和图3所示。

对降噪后的状态信号进行CEEMDAN分解，可以得到 7个IMF分量和1个趋势分量R，具体步骤包括：

 S201.CEEMDAN的分解过程采用以EEMD相同的方式获得各阶固有模态分量，首先设原始信号

 S202.对于每个加入高斯白噪声的原始信号

S203.根据这一公式可得第一阶固有模态分量

S204.计算获得第一阶段固有模态分量

 S205.向第一阶段分解的余量加入高斯白噪声，并对处理后的信号余量进行EMD分解，对第一阶段余量产生的第一个固有模态分量进行总体平均计算即可获得第二阶段的固有模态分量

计算获得第二阶段固有模态分量

 S206.对振动信号进行多次EMD分解直至第

 S207.得到CEEMDAN的第

 S208.直至分解后的余量无法继续分解，完成CEEMDAN分解流程，最终经CEEMDAN分解后的信号为：

 式中：

健康状态时经CEEMDAN分解的IMF分量如图4所示。

S3.通过相关系数法计算每一个IMF分量的相关系数，并选取相关系数最高的4个IMF分量；

S301. 在本发明的筛选流程中，经 CEEMDAN分解后的 固有模态分量可依据各分量的相关系数进行筛选，相关系数的数学表达式为：

式中：X，Y 代表两个样本变量，r(X,Y)表 示 两 样 本 变 量 之 间 的 相 关系 数 ，

IMF分量计算得到的相关系数如图5所示。

 S4.计算选取后的IMF分量的样本熵，并整合成样本集；

 S401.一般地，对于由N个数据组成的时间序列

S402.按序号组成一组维数为m的向量序列，

式中：

 S403.定义向量

式中：

 S404.对于给定的 

S405.对所有的

 S406.增加维数到m+1，计算 

S407.得到序列样本熵的理论值为：

 S408.但是在实际中N不可能为无穷大，而为一有限值，则样本熵的估计值为：

S409.根据上式可知，嵌入维数m和相似容限r的大小与样本熵的结果相关。当

计算的样本熵值如表2所示。

表2 机组不同运行状态振动信号样本熵特征值

 S5. 将特征样本集按照1：9划分为测试样本和训练样本，健康样本设置为样本1。并采用LSGSA-RBF对样本进行状态识别，设置初始参数如下：粒子维度为2，方式为8 ，种群数量为20，迭代次数为15，权重、初始速度、初始位置和加速度等参数根据样本复杂程度自适应调整。将测试样本输入到训练好的模型中进行识别率验证。LCGSA算法的全称为基于莱维飞行的混沌引力搜索算法，主要目的是解决GSA快速陷入局部最优和收敛时间长的问题。其基本原理如下：

 S501.GSA的初始化是通过随机分布的大质量个体吸引其他粒子来完成的，两个粒子间的作用力为：

式中：

 S502.其具体定义如下：

式中：

 S503.当个体在搜索空间中移动以寻找最优解时，其受到其它粒子的合力为：

式中：

 S504.依据牛顿第二定律，该质量的加速度为：

 S505.在进化过程中，粒子的速度和位置的更新方式为：

 S506.粒子的质量与适应度值有关，质量越大的粒子越接近最优，并且它对其他粒子的作用力相应地会更大，但移动速度较慢。粒子质量的计算方式为：

S507.从GSA的基本原理中不难看出，GSA的参数选取对算法的搜索性能起着至关重要的作用。一般在GSA算法中，万有引力常数

 S508.Levy飞行是一种特殊的的随机步长的方法，指的是步长的概率分布为重尾分布的随机行走，在随机行走的过程中有相对较大的概率出现大跨步，这种方法可以很好的绘制出莱维分布的方式。Levy的飞行公式如下所示：

 S509.Levy飞行的步长公式为：

式中：

 S510.参数的方差为：

S511.莱维飞行的作用主要是帮助GSA克服陷入局部最优的问题，它可以控制GSA中的万有引力常数

表3嵌入GSA的混沌序列函数

 S512.混沌归一化的表达式为：

式中：

 S513. 

式中：

 S514.将Levy和混沌序列与GSA的引力常数结合，得到LCGSA的引力常数

 S515.RBFNN构成的隐层空间可以将输入向量直接映射到隐层空间，因此不要按照权重连接，输入层和隐藏层之间的连接权值都为1。隐藏层（激活函数为RBF）实现输入向量的非线性投影，输出层（激活函数为线性函数）负责最终的线性化加权求和。径向基函数RBF是中心点径向对称、取值仅依赖于距中心点距离的非负实值函数。其基本架构如图6所示。

 S516.常用的RBF使用欧氏距离及高斯函数，其表达式为：

式中：

 S517.最终的输出

式中：

 S518.设

2.水泵机组故障状态时振动信号的特征提取与状态分类识别：

S1.将采集到的故障状态的水泵机组原始振动冲击信号采用MODWT小波变换进行降噪；去噪前后的振动冲击时水泵机组的实测信号波形图，如图7、图8所示。

S2. 对降噪后的状态信号进行CEEMDAN分解，可以得到 7个IMF分量和1个趋势分量R；

振动冲击时经CEEMDAN分解的IMF分量如图9所示。

 S3.通过相关系数法计算每一个IMF分量的相关系数，并选取相关系数最高的4个IMF分量；

S4.计算选取后的IMF分量的样本熵，并整合成样本集，如表2所示；

S5. 将特征样本集按照1：9划分为测试样本和训练样本，振动冲击样本设置为样本2。并采用LSGSA-RBF对样本进行状态识别，设置初始参数如下：粒子维度为2，方式为8 ，种群数量为20，迭代次数为15，权重、初始速度、初始位置和加速度等参数根据样本复杂程度自适应调整。将测试样本输入到训练好的模型中进行识别率验证。健康状态与故障状态分类识别准确率如图10所示。

为体现本发明的优越性，将本发明与CEEMDAN样本熵和RBF、EMD样本熵和LCGSA-RBF、CEEMD样本熵和LCGSA-RBF、VMD样本熵和LCGSA-RBF四组处理实例一与实例二的相同信号的结果作对比，对比图如图11所示。

为防止一次实验具有偶然性和随机性，将本发明的系统和四组对比系统分别运行20次，准确率如表4所示。 为进一步直观的表述不同模型的差异，采用递增平均值表示每一个模型的平均准确率变化情况，如图12所示。

表4 不同系统运行20次状态识别准确率及平均准确率

图12是递进平均值，也是对表14数据的另一种表述方式，递进平均值意思就是第二个点是第一个数据加第二个数据的平均值，第三个点是前前三组数据的平均值，第四个点是前四组数据的平均值，以此类推，这样曲线越接近于直线，说明数据波动性越小，进一步证明本发明所提出的机组诊断系统的优越性。

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。本领域的技术人员容易理解，以上所述仅为本发明专利的较佳实例而已，并不用以限制本发明专利，凡在本发明专利的精神和原则之内所作的任何修改，等同替代和改进等，均应包含在本发明专利的保护范围之内。本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。