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

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


背景技术：

2.水泵是输送液体或使液体增压的机械，它将原动力机的机械能或其他外部能量传送给液体，使液体能量增加，主要用来输送液体包括水、油、酸碱液、乳化液、悬乳液和液态金属等，也可输送液体、气体混合物以及含悬浮固体物的液体。水泵机组主要包括水泵、电动机、控制系统等部分，因其在工业系统中扮演的角色，导致机组运行工况变化频繁，会易产生机组振动现象，且振动异常可能导致机械故障、损失和安全事故。现有的关于水泵机组振动的研究越来越多，而研究的第一步就是需要对振动信号进行处理。
3.对水泵机组振动信号处理的一般步骤是将其进行特征提取以得到其中有用的特征信息，再使用分类器对其状态进行分类识别，以便于后续的故障诊断、预测和控制工作。对于水泵机组振动信号处理，目前有一部分采用传统的特征提取算法和分类器进行处理，时域分析、频域分析、小波变换、奇异值分解、支持向量机等，但这些方法特征提取的特征不够准确，在分类器的训练中容易出现过拟合等问题，导致结果的可靠性和稳定性差；另一部分采用深度学习方法，如卷积神经网络，递归神经网络等，虽然这些方法在特征提取和分类器构建上相对传统方法更具优势，但在样本计算和分类器训练时需要较大的数据量，而目前水泵机组振动信号数据较为有限，因此这种方法在实际应用中仍存在一定的局限性。因此，亟需开发一种适用性更广、可靠性更高的水泵机组振动信号特征提取与状态识别方法。


技术实现要素：

4.本发明的目的在于提供一种水泵机组振动信号特征提取与状态识别方法，解决现有技术方法中效果不理想、可靠性差且局限性大的问题。
5.为解决上述的技术问题，本发明采用以下技术方案：一种水泵机组振动信号特征提取与状态识别方法，其特征在于包括如下步骤：
6.s1. 分别采集水泵机组健康状态和故障状态下的振动信号，将采集到的水泵机组的振动信号采用modwt小波变换进行降噪；
7.s2. 对降噪后的状态信号进行ceemdan分解，得到 7个imf分量和1个最终余量r(n)；
8.s3.通过相关系数法计算每一个imf分量的相关系数，并选取相关系数最高的4个imf分量；
9.s4.计算选取后的imf分量的样本熵，并整合成特征样本集，并按状态分为健康样本和故障样本；
10.s5. 将特征样本集按照1：9划分为测试样本和训练样本，采用lcgsa-rbf方法对训练样本进行状态识别，得到训练好的模型；将测试样本输入到训练好的模型中进行识别准
确率验证；所述lcgsa-rbf方法具体为将训练样本输入到rbf神经网络对数据进行拟合，得到初始的rbf神经网络模型，再将rbf神经网络中的参数转化为搜索空间中的个体，根据其质量和位置计算出相互之间的作用力和速度，更新其位置和速度，在每次迭代中，利用gsa算法重新计算每个个体的质量和位置，更新搜索空间中的个体，在更新位置和速度时，引入莱维飞行和混沌理论的随机因素，根据搜索空间中的最优个体，得到最优的rbf神经网络模型。
11.更进一步的技术方案是所述步骤s1对水泵机组振动冲击信号采用最大重叠离散小波变换（modwt）进行去噪处理。modwt是dwt的升级版本，
12.s101.设x为n维向量，元素为实值时间序列，其中n是2的整数倍；再设是一个低通滤波器，其中，l表示滤波器的长度，该滤波器对于所有非零整数n均满足如下条件：
13.，
14.是与对应的高通滤波器；
15.s102. 设尺度
ꢀꢀ
下的尺度变换系数为 , , 其中 ，且有 ；
16.根据 mallat 算法，计算出尺度
ꢀꢀ
下的尺度变换系数和小波变换系数分别如下：
17.，
18.，
19.其中，mod表示两数相除取余数；其中
[0020][0021][0022]
s103.dwt中的尺度滤波器和小波滤波器关系如下：
[0023]
满足：
[0024][0025]
其中，
[0026][0027][0028]
s104.在尺度j下的插入 个零：
[0029][0030]
s105.根据mallat 算法，计算出尺度j下的尺度变换系数和小波变换系数分别如下：
[0031][0032]
。
[0033]
更进一步的技术方案是所述步骤s2对于降噪后的信号进行ceemdan分解，ceemdan是基于eemd所改进，加入经emd分解后含辅助噪声的 imf 分量，而不是将高斯白噪声信号直接添加在原始信号中，抑制了eemd的模态混叠和计算复杂度高以及虚假imf多的问题。
[0034]
ceemdan分解具体步骤如下：
[0035]
定义为经emd所产生的第个分量，为ceemdan产生的第个模态函数：
[0036]
s201.采用以eemd相同的方式获得各阶固有模态分量，首先设原始信号中加入标准差为的次高斯白噪声，即
[0037]
，
[0038]
s202.对于每个加入高斯白噪声的原始信号进行emd分解进而获得个固有模态分量并进行总体平均计算如下：
[0039]
，
[0040]
s203.根据上述公式计算得到第一阶固有模态分量
[0041]
，
[0042]
计算获得第一阶段固有模态分量后余量为：
[0043]
，
[0044]
s204.向第一阶段分解的余量加入高斯白噪声，并对处理后的信号余量进行emd分解，对第一阶段余量产生的第一个固有模态分量进行总体平均计算即获得第二阶段的固有模态分量：
[0045]
，
[0046]
计算获得第二阶段固有模态分量后余量为：
[0047]
，
[0048]
s205.依此，对振动信号进行多次emd分解直至第次，得到第个余量为：
[0049]
，
[0050]
得到ceemdan的第阶模态函数：
[0051]
，
[0052]
s206.直至分解后的余量无法继续分解，完成ceemdan分解流程，最终经ceemdan分解后的信号为：
[0053]
，
[0054]
式中：为信号经ceemdan分解得到的模态函数个数；为最终余量，表达式为：
[0055]
。
[0056]
更进一步的技术方案是所述步骤s3中所有imf分量的相关系数计算，并选取相关系数最好的4个分量，其具体步骤包括：
[0057]
在本发明的筛选流程中，经 ceemdan分解后的固有模态分量可依据各分量的相关系数进行筛选，相关系数的数学表达式为： [0058]
式中：x，y 代表两个样本变量，r(x,y)表 示 两 样 本 变 量 之 间 的 相 关 系 数 ， 表示变量 x 的方差，
ꢀꢀ
表示 y 的方差。
[0059]
更进一步的技术方案是所述步骤s4中对于相关系数最高的4个分量进行样本熵计算其具体步骤包括：
[0060]
对于由n个数据组成的时间序列：
[0061]
s401.按序号组成一组维数为m的向量序列，，则有：
[0062]
，
[0063]
式中：，这些向量代表从第i点开始的m个连续的x的值；
[0064]
s402.定义向量与之间的距离为两者对应元素中最大差值的绝对值，即：
[0065]
，
[0066]
式中：；
[0067]
s403.对于给定的 , 统计
ꢀꢀ
与
ꢀꢀ
之间距离小于等于的
ꢀꢀ
的数目，并记作；对于 , 则：
[0068]
，
[0069]
s404.对所有的取平均，得到：
[0070][0071]
s405.增加维数到m+1，计算
ꢀꢀ
与 距离小于等于的个数，得到和；这样，是两个序列在相似容限下匹配个点的概率，而是两个序列匹配m+1的点的概率；
[0072]
得到序列样本熵的理论值为：
[0073]
，
[0074]
但是在实际中n不可能为无穷大，而为一有限值，则样本熵的估计值为：
[0075]
，
[0076]
s406.根据上式可知，嵌入维数m和相似容限r的大小与样本熵的结果相关；设定
ꢀꢀ
，为离散的时间序列 的标准差。
[0077]
更进一步的技术方案是所述步骤s5中对将样本熵作为样本集输入到lcgsa-rbf算法中进行状态识别与分类，其具体步骤包括：
[0078]
设置初始参数如下：粒子维度为2，方式为8 ，种群数量为20，迭代次数为15，具体地：
[0079]
s501.是通过随机分布的大质量个体吸引其他粒子来完成gsa的初始化，两个粒子间的作用力为：
[0080]
，
[0081]
式中：和分别为粒子j和粒子i的质量；为一个极小的常量；粒子i与粒子j之间的欧式距离；为k时刻的万有引力常数；分别表示粒子i、粒子j的位置；
[0082]
s502.其具体定义如下：
[0083][0084]
式中：为进化参数；为初始时刻的引力常数；指数递减系数；ci 初始模拟量；
[0085]
s503.当个体在搜索空间中移动以寻找最优解时，其受到其它粒子的合力为：
[0086][0087]
式中：为变化区间（0，1）之间的随机数；
[0088]
s504.该质量的加速度为：
[0089][0090]
s505.在进化过程中，粒子的速度和位置的更新方式为：
[0091][0092]
s506.粒子的质量与适应度值有关，质量越大的粒子越接近最优，并且它对其他粒子的作用力相应地会更大，但移动速度较慢，粒子质量的计算方式为：
[0093][0094]
式中表示适应度参数、表示适应度参数最差值、表示适应度参数最优值、表示粒子i的质量、n代表粒子数量；
[0095]
s507.levy的飞行公式如下所示：
[0096][0097]
s508.levy飞行的步长公式为：
[0098][0099]
式中：； ；，都表示gamma函数；
[0100]
参数的方差为：
[0101][0102]
s509.混沌归一化的表达式为：
[0103][0104]
式中：和为混沌映射的范围；为混沌索引；为混沌归一化区间。
[0105]
的表达式为：
[0106][0107]
式中：和为自适应间隔；
[0108]
s510.将levy和混沌序列与gsa的引力常数结合，得到lcgsa的引力常数的表达式为：
[0109][0110]
s511.rbf使用欧氏距离及高斯函数，其表达式为：
[0111][0112]
式中：为欧式范数；为高斯函数的中心，
[0113]
s512.最终的输出的数学表达式为：
[0114][0115]
式中：为的输入样本；为网络隐藏层的中心；为从隐藏层到输出层的权重；为隐藏层中节点的数量；为的实际输出；为高斯函数的方差，
[0116]
样本的期望输出值，其表达式为：
[0117]
。
[0118]
工作原理：本发明提出了基于ceemdan样本熵和lcgsa-rbf的水泵机组振动信号特征提取与状态分类识别方法，首先对原始信号进行wodwt小波降噪，再对降噪后的信号进行ceemdan分解，得到一系列imf分量，同时计算各个分量的相关系数，选取相关系数最高的4个分量进行样本熵计算，将计算好的样本熵作为样本集输入到rbf神经网络对数据进行拟合，得到初始的rbf神经网络模型，之后，将rbf神经网络中的参数转化为搜索空间中的个体，并根据其质量和位置计算出相互之间的作用力和速度，更新其位置和速度，在每次迭代中，根据gsa算法的原理，重新计算每个个体的质量和位置，更新搜索空间中的个体，最后在更新位置和速度时，引入莱维飞行和混沌理论的随机因素，增加算法的随机性和全局搜索能力，最终，根据搜索空间中的最优个体，得到最优的rbf神经网络模型，从而实现对数据的优化拟合，使得模型能更加准确地将振动信号所对应的水泵机组状态进行识别。
[0119]
与现有技术相比，本发明的有益效果是：
[0120]
（1）相较于传统的原始信号分解方法，采用本发明提出的ceemdan分解方法所提取的样本熵更适用于模型的读取与辨识。
[0121]
（2）基于ceemdan样本熵和lcgsa-rbf的水泵机组振动信号特征提取与状态分类识别系统对于旋转机械的声音振动信号的识别优于传统系统，利用噪声传感器可以在和机组非接触的条件下，精准获取分析识别机组的突变信号，可以有效地监测机组的工作状态。
[0122]
（3）引入莱维飞行和混沌映射方法对gsa进行优化，增加算法的收敛速度和全局搜索能力，同时建立lcgsa-rbf的全新模型应用于水泵机组振动信号的识别应用中。
[0123]
（4）本发明设计全新模型解决水泵机组故障诊断效率低、反应速度慢等问题，在水泵机组特征提取领域具有很强的优越性和实用价值。
[0124]
（5）本发明的方法对于促进水泵机组的安全可靠运行具有重要而积极的意义，通过提供自动化的振动信号分类和识别，可以有效地监测机组的工作状态，及时发现异常情况并采取相应的措施，从而提高机组的安全性和可靠性，延长其使用寿命。
附图说明
[0125]
图1为本发明的流程图。
[0126]
图2为健康状态时水泵机组的实测信号波形图。
[0127]
图3为去噪后的健康状态下水泵机组的实测信号波形图。
[0128]
图4为健康状态时经ceemdan分解的imf分量图。
[0129]
图5为健康状态时和振动冲击时两组imf分量相关系数图。
[0130]
图6为rbf神经网络基本架构图。
[0131]
图7为振动冲击时水泵机组的实测信号波形图。
[0132]
图8为去噪后的振动冲击时水泵机组的实测信号波形图。
[0133]
图9为振动冲击时经ceemdan分解的imf分量图。
[0134]
图10为状态分类识别准确率图。
[0135]
图11为本发明与ceemdan样本熵和rbf、emd样本熵和lcgsa-rbf、ceemd样本熵和lcgsa-rbf、vmd样本熵和lcgsa-rbf四组模型对比图。
[0136]
图12为5种模型20次运行后的递增平均值图。
具体实施方式
[0137]
为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述，但以下实施例仅是说明性的，本发明的保护范围并不受这些实施例的限制。
[0138]
该水泵机组额定功率55kw，额定电压400v，额定水头10m，转速为81.6r/min，因实际机组中很难发生故障，本发明采取机组运行状态下对顶盖进行连续敲击模拟故障条件下振动信号的产生，并采用非接触试噪声传感器crysoud进行实时采集，该传感器内置先进的dsp处理器。集传声器、前置放大器和数据采集卡于一体，在结构布置紧密的前提下，又保证了数据处理高校准确。用噪声传感器共采集得到机组2种不同状态下的振动信号，分别为故障和健康两组样本，每种状态的数据集为1000个样本，每个样本的点数为4096。传感器参数
如表1所示，数据采集卡的型号为研华idaq801。
[0139]
表1噪声传感器cry2301参数
[0140]
参数名称规格单位采样速率48khz测量频率范围10-20000hz标准测量范围25-130dba测量动态范围≥110dba通讯接口usbaudio+usbhid尺寸25
×
115mm
[0141]
下面结合实例描述本发明所述的一种基于ceemdan样本熵和lcgsa-rbf的水泵机组振动信号特征提取与状态识别方法，针对实例中的水泵机组的振动冲击信号进行状态分类识别，并将结果与ceemdan样本熵和rbf、emd样本熵和lcgsa-rbf、ceemd样本熵和lcgsa-rbf、vmd样本熵和lcgsa-rbf四种系统处理同一种信号的结果进行对比对比，以验证本发明的优越性，包括如下所述具体过程：
[0142]
1.水泵机组健康状态时振动信号的特征提取与状态分类识别。
[0143]
如图1所示，本发明实例提供的基于ceemdan样本熵和lcgsa-rbf的水泵机组振动信号特征提取与状态识别方法，包括如下步骤：
[0144] s1.将采集到的水泵机组原始振动冲击信号采用modwt小波变换进行降噪；
[0145]
s101.modwt是dwt的升级版本，其具体原理如下：
[0146]
s102.设x为n维向量，它的元素为实值时间序列，其中n是2的整数倍；再设是一个低通滤波器，其中，l表示滤波器的长度，该滤波器对于所有非零整数n均满足如下条件：
[0147][0148] 是与对应的高通滤波器。
[0149]
s103.设尺度下的尺度变换系数为,,其中,且有。
[0150]
则根据mallat算法,可以计算出尺度下的尺度变换系数(近似)和小波变换系数(细节)分别如下:
[0151][0152][0153]
其中，mod表示两数相除取余数。
[0154]
s104.且：
[0155][0156]
s105.modwt与dwt中的尺度滤波器和小波滤波器关系如下：
[0157]
[0158][0159]
s106.且同样要满足：
[0160][0161][0162]
或者
[0163]
s107.为了避免由于抽取采样导致变换后样本点数减少的情况，即在尺度j下的插入 个零：
[0164][0165] s108.根据mallat算法，可以计算出尺度j下的尺度变换系数（近似）和小波变换系数（细节）分别如下：
[0166][0167][0168]
去噪前后的波形图分别如图2和图3所示。
[0169]
对降噪后的状态信号进行ceemdan分解，可以得到 7个imf分量和1个趋势分量r，具体步骤包括：
[0170] s201.ceemdan的分解过程采用以eemd相同的方式获得各阶固有模态分量，首先设原始信号中加入标准差为=的次高斯白噪声。
[0171][0172] s202.对于每个加入高斯白噪声的原始信号进行emd分解进而获得个固有模态分量并进行总体平均计算如下：
[0173][0174]
s203.根据这一公式可得第一阶固有模态分量
[0175][0176]
s204.计算获得第一阶段固有模态分量后余量为：
[0177][0178] s205.向第一阶段分解的余量加入高斯白噪声，并对处理后的信号余量进行emd分解，对第一阶段余量产生的第一个固有模态分量进行总体平均计算即可获得第二阶段的固有模态分量：
[0179]
计算获得第二阶段固有模态分量后余量为：
[0180][0181] s206.对振动信号进行多次emd分解直至第次，得到第个余量信号为：
[0182][0183] s207.得到ceemdan的第阶模态函数：
[0184][0185] s208.直至分解后的余量无法继续分解，完成ceemdan分解流程，最终经ceemdan分解后的信号为：
[0186][0187] 式中：为信号经ceemdan分解得到的模态函数个数；为最终余量信号，表达式为：
[0188][0189]
健康状态时经ceemdan分解的imf分量如图4所示。
[0190]
s3.通过相关系数法计算每一个imf分量的相关系数，并选取相关系数最高的4个imf分量；
[0191]
s301. 在本发明的筛选流程中，经 ceemdan分解后的 固有模态分量可依据各分量的相关系数进行筛选，相关系数的数学表达式为：
[0192][0193]
式中：x，y 代表两个样本变量，r(x,y)表 示 两 样 本 变 量 之 间 的 相 关 系 数 ， 表示变量 x 的方差，
ꢀꢀ
表示 y 的方差。
[0194]
imf分量计算得到的相关系数如图5所示。
[0195] s4.计算选取后的imf分量的样本熵，并整合成样本集；
[0196] s401.一般地，对于由n个数据组成的时间序列，样本熵的计算方法为：
[0197]
s402.按序号组成一组维数为m的向量序列，。则有：
[0198][0199]
式中：。这些向量代表从第i点开始的m个连续的x的值。
[0200] s403.定义向量与之间的距离为两者对应元素中最大差值的绝对值。即：
[0201][0202]
式中：。
[0203] s404.对于给定的 , 统计
ꢀꢀ
与
ꢀꢀ
之间距离小于等于的
ꢀꢀ
的数目，并记作。对于 , 则：
[0204][0205]
s405.对所有的取平均，得到：
[0206][0207] s406.增加维数到m+1，计算
ꢀꢀ
与 距离小于等于的个数，得到和。这样，是两个序列在相似容限下匹配个点的概率，而是两个序列匹配m+1的点的概率。
[0208]
s407.得到序列样本熵的理论值为：
[0209] s408.但是在实际中n不可能为无穷大，而为一有限值，则样本熵的估计值为：
[0210]
s409.根据上式可知，嵌入维数m和相似容限r的大小与样本熵的结果相关。当或或(为离散的时间序列为离散的时间序列的标准差)时,所得样本熵才有统计意义。本发明选择。
[0211]
计算的样本熵值如表2所示。
[0212]
表2机组不同运行状态振动信号样本熵特征值
[0213][0214] s5.将特征样本集按照1：9划分为测试样本和训练样本，健康样本设置为样本1。并采用lsgsa-rbf对样本进行状态识别，设置初始参数如下：粒子维度为2，方式为8，种群数量为20，迭代次数为15，权重、初始速度、初始位置和加速度等参数根据样本复杂程度自适应调整。将测试样本输入到训练好的模型中进行识别率验证。lcgsa算法的全称为基于莱维飞行的混沌引力搜索算法，主要目的是解决gsa快速陷入局部最优和收敛时间长的问题。其基本原理如下：
[0215] s501.gsa的初始化是通过随机分布的大质量个体吸引其他粒子来完成的，两个粒子间的作用力为：
[0216][0217]
式中：和分别为粒子j和粒子i的质量；为一个极小的常量；粒子i与粒子j之间的欧式距离；为k时刻的万有引力常数。
[0218] s502.其具体定义如下：
[0219][0220]
式中：为进化参数；为初始时刻的引力常数；
[0221] s503.当个体在搜索空间中移动以寻找最优解时，其受到其它粒子的合力为：
[0222][0223]
式中：为变化区间【0，1】之间的随机数。
[0224] s504.依据牛顿第二定律，该质量的加速度为：
[0225][0226] s505.在进化过程中，粒子的速度和位置的更新方式为：
[0227][0228] s506.粒子的质量与适应度值有关，质量越大的粒子越接近最优，并且它对其他粒子的作用力相应地会更大，但移动速度较慢。粒子质量的计算方式为：
[0229][0230]
s507.从gsa的基本原理中不难看出，gsa的参数选取对算法的搜索性能起着至关重要的作用。一般在gsa算法中，万有引力常数会被设定一个固定值，很大程度上限制了算法的搜索效率，且无法保证可以获得全局最优解。
[0231] s508.levy飞行是一种特殊的的随机步长的方法，指的是步长的概率分布为重尾分布的随机行走，在随机行走的过程中有相对较大的概率出现大跨步，这种方法可以很好的绘制出莱维分布的方式。levy的飞行公式如下所示：
[0232][0233] s509.levy飞行的步长公式为：
[0234][0235]
式中：； ；，都表示gamma函数。
[0236] s510.参数的方差为：
[0237][0238]
s511.莱维飞行的作用主要是帮助gsa克服陷入局部最优的问题，它可以控制gsa中的万有引力常数，以解决gsa的缺陷。将levy-gsa算法与混沌优化算法结合，能够提高levy-gsa算法的收敛能力和局部搜索能力。在levy-gsa算法中，惯性权重在（0-1）之间随机选取，直接应用上式来计算粒子总力。而本发明采用混沌映射方法来替换随机值选取过程，其中每个粒子的混沌映射都遵循相应的顺序。具体的混沌序列随机行为分布方式如图6所示，嵌入gsa的混沌序列函数如表3所示。
[0239]
表3嵌入gsa的混沌序列函数
[0240][0241] s512.混沌归一化的表达式为：
[0242][0243]
式中：和为混沌映射的范围；为混沌索引；为混沌归一化区间。
[0244] s513. 的表达式为：
[0245][0246]
式中：和为自适应间隔。
[0247] s514.将levy和混沌序列与gsa的引力常数结合，得到lcgsa的引力常数的表达式为：
[0248][0249] s515.rbfnn构成的隐层空间可以将输入向量直接映射到隐层空间，因此不要按照权重连接，输入层和隐藏层之间的连接权值都为1。隐藏层（激活函数为rbf）实现输入向量的非线性投影，输出层（激活函数为线性函数）负责最终的线性化加权求和。径向基函数rbf是中心点径向对称、取值仅依赖于距中心点距离的非负实值函数。其基本架构如图6所示。
[0250] s516.常用的rbf使用欧氏距离及高斯函数，其表达式为：
[0251][0252]
式中：为欧式范数；为高斯函数的中心。
[0253] s517.最终的输出的数学表达式为：
[0254][0255]
式中：为的输入样本；为网络隐藏层的中心；为从隐藏层到输出层的权重；为隐藏层中节点的数量；为的实际输出；为高斯函数的方差.
[0256] s518.设为样本的期望输出值，其表达式为：
[0257][0258]
2.水泵机组故障状态时振动信号的特征提取与状态分类识别：
[0259]
s1.将采集到的故障状态的水泵机组原始振动冲击信号采用modwt小波变换进行降噪；去噪前后的振动冲击时水泵机组的实测信号波形图，如图7、图8所示。
[0260]
s2. 对降噪后的状态信号进行ceemdan分解，可以得到 7个imf分量和1个趋势分量r；
[0261]
振动冲击时经ceemdan分解的imf分量如图9所示。
[0262] s3.通过相关系数法计算每一个imf分量的相关系数，并选取相关系数最高的4个imf分量；
[0263]
s4.计算选取后的imf分量的样本熵，并整合成样本集，如表2所示；
[0264]
s5. 将特征样本集按照1：9划分为测试样本和训练样本，振动冲击样本设置为样本2。并采用lsgsa-rbf对样本进行状态识别，设置初始参数如下：粒子维度为2，方式为8 ，种群数量为20，迭代次数为15，权重、初始速度、初始位置和加速度等参数根据样本复杂程度自适应调整。将测试样本输入到训练好的模型中进行识别率验证。健康状态与故障状态分类识别准确率如图10所示。
[0265]
为体现本发明的优越性，将本发明与ceemdan样本熵和rbf、emd样本熵和lcgsa-rbf、ceemd样本熵和lcgsa-rbf、vmd样本熵和lcgsa-rbf四组处理实例一与实例二的相同信号的结果作对比，对比图如图11所示。
[0266]
为防止一次实验具有偶然性和随机性，将本发明的系统和四组对比系统分别运行20次，准确率如表4所示。 为进一步直观的表述不同模型的差异，采用递增平均值表示每一个模型的平均准确率变化情况，如图12所示。
[0267]
表4 不同系统运行20次状态识别准确率及平均准确率
[0268]
ceemdan-lcgsa-rbfemd-lcgsa-rbfceemd-lcgsa-rbfvmd-lcgsa-rbfceemdan-rbf1100%94.5%98.5%99.5%87.5%2100%96%95.5%99.5%85.5%399.5%92.5%95%99%84.5%4100%91.5%96%95.5%85%
599.5%97%96%99%86%699%95%94%94%89.5%7100%95.5%98%97.5%84%899.5%92%92%98%86%999%94.5%90%99.5%87.5%1099%95%90%97.5%90.5%1199%94%89.5%98%85.5%1299.5%91.5%86%94.5%83%1399%97%92%98%89%14100%98%92.5%96%86.5%15100%95.5%96%92%84%16100%93.5%96.5%99.5%91.5%1799.5%95%93%95%92%1899%94%94.5%96%87%1999%90%95%95.5%86%20100%92%91%97%89.5%平均准确率99.5%94.2%93.55%97.025%87%
[0269]
图12是递进平均值，也是对表14数据的另一种表述方式，递进平均值意思就是第二个点是第一个数据加第二个数据的平均值，第三个点是前前三组数据的平均值，第四个点是前四组数据的平均值，以此类推，这样曲线越接近于直线，说明数据波动性越小，进一步证明本发明所提出的机组诊断系统的优越性。
[0270]
以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。本领域的技术人员容易理解，以上所述仅为本发明专利的较佳实例而已，并不用以限制本发明专利，凡在本发明专利的精神和原则之内所作的任何修改，等同替代和改进等，均应包含在本发明专利的保护范围之内。本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

技术特征：
1.一种水泵机组振动信号特征提取与状态识别方法，其特征在于包括如下步骤：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神经网络模型。2.根据权利1所述的一种水泵机组振动信号特征提取与状态识别方法，其特征在于：所述步骤s1的modwt小波变换进行降噪具体步骤如下：s101.设x为n维向量，元素为实值时间序列，其中n是2的整数倍；再设是一个低通滤波器，其中，l表示滤波器的长度，该滤波器对于所有非零整数n均满足如下条件：，是与对应的高通滤波器；s102. 设尺度
ꢀꢀ
下的尺度变换系数为 , , 其中 ，且有 ；根据 mallat 算法，计算出尺度
ꢀꢀ
下的尺度变换系数和小波变换系数分别如下：，，其中，mod表示两数相除取余数；其中，，，s103.dwt中的尺度滤波器和小波滤波器关系如下：，，满足：
，其中，，，s104.在尺度j下的插入 个零：，s105.根据mallat 算法，计算出尺度j下的尺度变换系数和小波变换系数分别如下：，。3.根据权利1所述的一种水泵机组振动信号特征提取与状态识别方法，其特征在于：所述步骤s2的ceemdan分解具体步骤如下：定义为经emd所产生的第个分量，为ceemdan产生的第个模态函数：s201.采用以eemd集合经验模态分解的方式获得各阶固有模态分量，首先设原始信号中加入标准差为的次高斯白噪声，即，s202.对于每个加入高斯白噪声的原始信号进行emd分解进而获得个固有模态分量并进行总体平均计算如下：，s203.根据上述公式计算得到第一阶固有模态分量，计算获得第一阶段固有模态分量后余量为：，s204.向第一阶段分解的余量加入高斯白噪声，并对处理后的信号余量进行emd经验模态分解，对第一阶段余量产生的第一个固有模态分量进行总体平均计算即获得第二阶段的固有模态分量：，其中，为经emd所产生的第个分量；表示加入高斯白噪声的标准差；计算获得第二阶段固有模态分量后余量为：，
s205.依此，对振动信号进行多次emd分解直至第次，得到第个余量为：，得到ceemdan的第阶模态函数：，其中，表示加入高斯白噪声的标准差s206.直至分解后的余量无法继续分解，完成ceemdan分解流程，最终经ceemdan分解后的信号为：，式中：为信号经ceemdan分解得到的模态函数个数；为最终余量，表达式为：。4.根据权利1所述的一种水泵机组振动信号特征提取与状态识别方法，其特征在于：所述步骤s3中的相关系数法为皮尔逊相关系数的计算方式，其产生的相关系数范围在[-1,1]之间，相关系数的绝对值越大，变量之间的线性相关程度越大；反之，相关系数越小，变量之间线性相关程度越小，相关系数的数学表达式为：，式中：x，y 代表两个样本变量，r(x,y)表 示 两 样 本 变 量 之 间 的 相 关 系 数 ， 表示变量 x 的方差，
ꢀꢀ
表示 y 的方差，cov表示协方差。5.根据权利1所述的一种水泵机组振动信号特征提取与状态识别方法，其特征在于：所述步骤s4中样本熵的具体步骤如下：对于由n个数据组成的时间序列：s401.按序号组成一组维数为m的向量序列，，则有：，式中：，这些向量代表从第i点开始的m个连续的x的值；s402.定义向量与之间的距离为两者对应元素中最大差值的绝对值，即：，式中：；s403.对于给定的 , 统计
ꢀꢀ
与
ꢀꢀ
之间距离小于等于的
ꢀꢀ
的数目，并记作；对于 , 则：，s404.对所有的取平均，得到：，
s405.增加维数到m+1，计算
ꢀꢀ
与 距离小于等于的个数，得到和；这样，是两个序列在相似容限下匹配个点的概率，而是两个序列匹配m+1的点的概率；s406.得到序列样本熵的理论值为：，s407.但是在实际中n不可能为无穷大，而为一有限值，则样本熵的估计值为：，s408.根据上式可知，嵌入维数m和相似容限r的大小与样本熵的结果相关；设定
ꢀꢀ
，为离散的时间序列为离散的时间序列的标准差。6.根据权利1所述的一种水泵机组振动信号特征提取与状态识别方法，其特征在于：所述步骤5中采用lcgsa-rbf对样本进行状态识别，设置初始参数如下：粒子维度为2，方式为8 ，种群数量为20，迭代次数为15，具体地：s501.是通过随机分布的大质量个体吸引其他粒子来完成gsa的初始化，两个粒子间的作用力为：，式中：和分别为粒子j的主动引力质量和粒子i的被动引力质量；为一个极小的常量；粒子i与粒子j之间的欧式距离；为k时刻的万有引力常数；分别代表粒子i、粒子j的位置；s502.其具体定义如下：，式中：为进化参数；为初始时刻的引力常数；为指数递减系数；ci 初始模拟量；s503.当个体在搜索空间中移动以寻找最优解时，其受到其它粒子的合力为：，式中：为变化区间（0，1）之间的随机数；初始值为m，终值为1；s504.该质量的加速度为：，s505.在进化过程中，粒子的速度和位置的更新方式为：，式中的表示粒子的速度；
s506.粒子的质量与适应度值有关，质量越大的粒子越接近最优，并且它对其他粒子的作用力相应地会更大，但移动速度较慢，粒子质量的计算方式为：，式中表示适应度参数、表示适应度参数最差值、表示适应度参数最优值、表示粒子i的质量、n代表粒子数量；s507.levy的飞行公式如下所示：，式中表示飞行步长；s508.levy飞行的步长公式为：，式中：； ；，都表示正态分布函数；参数的方差为：，s509.混沌归一化的表达式为：，式中：和为混沌映射的范围；为混沌索引；为k时刻的混沌指数范围；为混沌归一化区间；的表达式为：，式中：和为自适应间隔；ci 初始模拟量;s510.将levy和混沌序列与gsa的引力常数结合，得到lcgsa的引力常数的表达式为：，其中，为搜索空间维度为 d、候选解数量n为的levy飞行函数；s511.rbf使用欧氏距离及高斯函数，其表达式为：，
式中：为欧式范数；为高斯函数的中心，s512.最终的输出的数学表达式为：，式中：为的输入样本；为网络隐藏层的中心；为从隐藏层到输出层的权重；为隐藏层中节点的数量；为的实际输出；为高斯函数的方差，。

技术总结
本发明公开了一种水泵机组振动信号特征提取与状态识别方法，涉及振动冲击信号处理技术领域。将采集到的水泵机组原始振动冲击信号采用MODWT小波变换进行降噪；对降噪后的状态信号进行CEEMDAN分解，得到7个IMF分量和1个最终余量R(n)；通过相关系数法计算每一个IMF分量的相关系数，并选取相关系数最高的4个IMF分量；计算选取后的IMF分量的样本熵，并整合成特征样本集；采用LCGSA-RBF对训练样本进行状态识别，得到训练好的模型。采用CEEMDAN分解方法提取的样本熵更适用于模型的读取与辨识，引入莱维飞行和混沌映射方法对GSA进行优化，增加算法的收敛速度和全局搜索能力，解决现有方法效果不理想、可靠性差且局限性大的问题。可靠性差且局限性大的问题。可靠性差且局限性大的问题。


技术研发人员：曾云 李想 钱晶 刀方 王芳芳 张晓旭
受保护的技术使用者：昆明理工大学
技术研发日：2023.08.02
技术公布日：2023/9/6