CN111611533B - 基于自适应vmd和改进功率谱的磨机负荷特征提取方法 - Google Patents

Abstract

本发明公开了一种基于自适应VMD和改进功率谱的磨机负荷特征提取方法，步骤1，对球磨机筒体振动信号x(i)进行EMD分解、筛选、重构以及VMD，得到模态分量{u(i)}；步骤2：基于模态分量{u(i)}得到改进自相关函数{r(τ)}；采用序列延拓，将模态分量{u(i)}进行改进自相关得到改进自相关函数{r(τ)}；步骤3：对{r(τ)}功率谱用Nuttall自卷积窗能量重心法进行频率校正，将校正后的功率谱最大值对应频率作为表征球磨机负荷变化的特征量。本方法能够有效提取磨矿过程中球磨机内部负荷特征，实现球磨机负荷状态的准确识别，为磨矿优化控制和效率提高提供准确、可靠的依据。

Description

基于自适应VMD和改进功率谱的磨机负荷特征提取方法

技术领域

本发明涉及一种基于自适应VMD和改进功率谱的磨机负荷特征提取方法。

背景技术

球磨机作为磨矿过程中最重要的核心生产设备，其效率的高低直接影响整个选矿厂的生产水平。球磨机负荷是指球磨机筒体内矿、研磨介质、水、矿浆等物料的总和，准确获取球磨机负荷状态，能为球磨机的优化控制提供重要依据。而实际生产中，球磨机具有时变性、非线性、随机干扰大等特性，使得球磨机负荷特征难以准确获取。

为准确获取球磨机负荷状态，近年来国内外领域研究者提出多种方法与先进技术。澳大利亚联邦科学和工业研究组织(CSIRO)提出基于球磨机筒体振动信号分析方法，相比传统功率法、压差法、噪声法，测量方式更直接，具有灵敏度高、抗干扰性强的优点。

球磨机筒体振动信号最主要和最本质的原因是筒体中钢球抛落过程中对筒壁的冲击而形成的，而振动的大小受到筒内负荷多少的影响。实际测量中的球磨机筒体振动信号，除了钢球冲击筒体的特征信号外，还有其他固有的一些振动信号，比如钢球之间的撞击振动、球磨机的机械振动、故障振动信号等，这些信号不随球磨机负荷的改变而改变。

现有经验模态分解(Empirical Mode Decomposition，EMD)和集成EMD(EnsembleEMD，EEMD)能对球磨机的振动信号进行自适应分解，进而提取负荷特征量，但EMD和EEMD均属于递归分解，存在模态混叠和端点效应的问题，导致从模态分量中提取的负荷特征不明显，负荷状态识别率低。

因此，有必要设计一种新的磨机负荷特征提取方法。

发明内容

本发明所要解决的技术问题是提供一种基于自适应变分模态分解(VariationalMode Decomposition，VMD)和改进功率谱的磨机负荷特征提取方法能够有效提取磨矿过程中球磨机内部负荷特征，实现球磨机负荷状态的准确识别，为磨矿优化控制和效率提高提供准确、可靠的依据。

发明的技术解决方案如下：

一种基于自适应VMD和改进功率谱的磨机负荷特征提取方法，包括以下步骤：

步骤1，对球磨机筒体振动信号x(i)进行EMD分解、筛选、重构以及VMD，得到模态分量{u(i)}；

步骤2：基于模态分量{u(i)}得到改进自相关函数{r(τ)}；

采用序列延拓，将模态分量{u(i)}进行改进自相关得到改进自相关函数{r(τ)}；

步骤3：对改进自相关函数{r(τ)}功率谱用Nuttall自卷积窗能量重心法进行频率校正，将校正后的功率谱最大值对应频率{f_k}作为表征球磨机负荷变化的特征量。

步骤1的过程为：

利用EMD对球磨机筒体振动信号x(i)进行分解得到本征模态分量(IntrinsicMode Function，IMF)，依据球磨机负荷变化与筒体振动信号冲击成分强相关原理结合峰度参数ρ滤除与球磨机负荷变化不相关的IMF，由峰度变化趋势选取峰度参数ρ大于85％所对应的K个IMF进行信号重构并同时确定VMD的分解层数为K，对重构信号进行VMD(变分模态分解)，得到K个模态分量{u(i)}。

球磨机筒体振动信号为x(i)，i＝1，2，...n，n为信号序列长度；

利用EMD的快速自适应性能对球磨机筒体振动信号x(i)进行分解得到一系列模态函数{u_e(i)}，求取每个IMF的峰度{Ku}：

式中，

为u_e(i)的平均值。

将每个IMF的峰度{Ku}按从大到小的顺序排序，并计算峰度参数ρ：

式中，Kus为所有峰度的和，k_max为信号经EMD后得到的最大模态数；

选取峰度参数ρ大于85％所对应的K个IMF进行信号重构得到x_re(i)：

对重构信号x_re(i)进行VMD，设定分解的模态数为K，得到K个模态分量{u(i)}。VMD的变分分解是现有的技术；VMD分解的过程对应公式1-公式6。

步骤2中，对于离散序列u(i)，改进的自相关函数

为：

式中，τ为时延，其取值范围为[1，n]，U(i)为

原来的自相关函数r(τ)式为：

改进过程可理解为：对离散序列进行序列延拓，避免自相关函数r(τ)时延τ增加，使得求和项减少，导致波形峰值幅度衰减问题。

步骤3中，根据球磨机筒体振动信号模态分量频谱特征，综合考虑窗函数旁瓣峰值电平与衰减速率，采用Nuttall自卷积窗能量重心法，对改进自相关函数{r(τ)}功率谱进行频率校正，将校正后的功率谱最大值对应频率{f_k}作为表征球磨机负荷变化的特征量。

步骤3中，Nuttall自卷积窗能量重心法的校正式为：

式中，f_k为前文所提的功率谱最大值对应频率，f_s为采样频率，n为采样点数，k₀为谐波信号经离散傅里叶变换得到的峰值谱线，此处谐波信号为改进自相关函数

P_N-2(k₀+i)为加窗功率谱第i条谱线功率值，P_N-2(ω)为Nuttall自卷积窗的功率谱函数，表达式如下：

式中，M为Nuttall窗函数的项数；b_m为Nuttall的窗系数，且满足

的约束条件，综合考虑旁瓣峰值电平和衰减速率，选用4项3阶Nuttall窗，即M＝4，b₀＝0.338946，b₁＝0.481973，b₂＝0.161054，b₃＝0.018027，W_R(·)为矩形窗的频谱函数，表达式为

VMD能有效对振动信号进行自适应分解，并能有效克服模态混叠和端点效应，其本质是构造和解决变分问题，对于原始信号f(t)，其对应的约束变分模型为：

式中，δ(t)表示冲击函数，‘*’为卷积符号，{u_k}＝{u₁，u₂，...，u_k}为k个本征模态分量，{ω_k}＝{ω₁，ω₂，...，ω_k}为各模态中心频率。

引入拉格朗日乘数算子λ和惩罚因子α将有约束性条件变为非约束性的变分问题。扩展的拉格朗日表达式如下：

利用乘法算子的交替方向法迭代求解式公式(12)的“鞍点”，得到模态分量u_k与中心频率ω_k分别为：

每次解得模态函数和中心频率后，利用下式对拉格朗日乘数运算符λ进行更新：

式中，σ为更新因子。

按照上述更新迭代，直到满足如下收敛条件：

式中，ε是收敛标准的容差，使用标准VMD的默认值。经过上述计算，VMD的变分分解过程完成，得到k个本征模态分量{u_k}。

利用VMD算法对信号分解，需要预先设定本征模态分量分解的数量，分解模态数量太少，会导致信息缺失，取值太多则会产生频率混叠现象，并且当数据序列过大同时分解模态数过多时，会降低其分解速度，因此需要选取合适的分解模态数K。

将EMD的快速自适应性能对球磨机筒体振动信号x(i)进行分解得到一系列IMF，每个IMF的峰度为Ku：

式中，u_e(i)为EMD分解得到IMF，n为序列长度，

为u_e(i)的平均值。

将每个IMF的峰度{Ku}按从大到小的顺序排序，得到峰度参数ρ：

式中，Kus为所有峰度的和，k_max为信号经EMD后得到的最大模态数。

依据球磨机负荷变化与筒体振动信号冲击成分强相关原理结合峰度参数ρ滤除与球磨机负荷变化不相关的IMF，由峰度变化趋势选取峰度参数ρ大于85％所对应的K个IMF进行信号重构得到x_re(i)：

对重构x_re(i)信号进行VMD，设定分解的模态数为K，得到K个模态分量{u(i)}。

对于离散序列u(i)，其自相关函数{r(τ)}式为：

式中，τ为时延。

对于任一τ＞0，可利用数据只有n-τ个，因此实际运算时，公式(10)变为

由公式(11)可看出，r(τ)随着τ的增加，其求和项逐渐减少，导致波形的峰值幅度逐渐衰减，为检测峰值点和峰值点清晰度带来困难。对于非周期函数，其自相关函数也不具周期性，峰值点会随着τ的增加迅速下降。为有效改正这一缺点，从充分利用采样数据的角度，对自相关函数进行改进。

改进自相关函数将序列长度从n延拓到2n，当i＞n时，令u(i)＝u(i-n)。再以公式(10)计算自相关函数，保证τ在[1，n]间取值时，做相关求和的数据长度为n。

将得到的K个模态分量{u(i)}进行改进自相关，得到自相关函数{r(τ)}。

能量重心法依据窗离散频谱的能量重心无限靠近坐标原点的特性，对频率进行校正，能量重心法校正式为：

式中，f_k为校正频率，f_s为采样频率，n为采样点数，k₀为谐波信号经离散傅里叶变换得到的峰值谱线，P(k₀+i)为加窗功率谱第i条谱线功率值，l为窗函数的旁瓣谱线数。

能量重心法的频谱分析精度易受到窗函数的影响，依据球磨机筒体振动信号模态分量频谱特征，比较现有窗函数性能，构建具有更低的旁瓣峰值电平、更快衰减速率Nuttall自卷积窗，改进能量重心法，以提高球磨机筒体振动信号模态分量频谱分析的精度。

Nuttall窗的离散表达式为：

式中，M为Nuttall窗函数的项数；b_m为Nuttall的窗系数，且满足

的约束条件，综合考虑，选用4项3阶Nuttall窗。

令Nuttall窗的频谱模函数W_N(ω)为：

式中，W_R()为矩形窗的频谱函数。

由余弦自卷积窗定义，可得Nuttall自卷积窗的时域表达式为：

w_N-2(n)＝[w(n)]²   …………………(公式15)

根据卷积定理，函数的时域卷积相当于频域的乘积，可得Nuttall自卷积窗的频谱模函数：

W_N-2(ω)＝[W_N(ω)]²   …………………(公式16)

令P_N-2(ω)为Nuttall自卷积窗的功率谱函数，则Nuttall自卷积窗窗谱能量表达式为：

由公式(12)和公式(17)可得Nuttall自卷积窗能量重心法频率校正式为：

为保证主瓣内较大的谱线均被利用，又减少旁瓣对结果的影响，采用l＝7谱线进行频率校正。

所述球磨机负荷特征提取方法还包括：

依据球磨机负荷变化与筒体振动信号冲击成分强相关原理，确定模态分量{u(i）}进行改进自相关得到改进自相关函数{r(τ)}的功率谱最大值谱线{k₀}，以公式(18)计算Nuttall自卷积窗能量重心法校正后的功率谱最大值对应频率{f_k}作为表征球磨机负荷变化的特征量。

有益效果：

本发明的基于自适应VMD和改进功率谱的磨机负荷特征提取方法，步骤1，利用EMD对球磨机筒体振动信号x(i)进行分解得到IMF，依据球磨机负荷变化与筒体振动信号冲击成分强相关原理结合峰度参数ρ滤除与球磨机负荷变化不相关的IMF，由峰度变化趋势选取峰度参数ρ大于85％所对应的K个IMF进行信号重构，对重构信号进行VMD，得到K个模态分量{u(i)}；步骤2，采用序列延拓，对自相关函数进行改进，将模态分量{u(i)}进行改进自相关得到改进自相关函数{r(τ)}；步骤3，依据球磨机筒体振动信号模态分量频谱特征，综合考虑窗函数旁瓣峰值电平与衰减速率，构建Nuttall自卷积窗能量重心法，对改进自相关函数{r(τ）}功率谱进行频谱校正，将校正后的功率谱最大值对应频率{f_k}作为表征球磨机负荷变化的特征量。本方法能够有效提取磨矿过程中球磨机内部负荷特征，实现球磨机负荷状态的准确识别，为磨矿优化控制和效率提高提供准确、可靠的依据。

本发明基于现有EMD与EEMD方法，依据球磨机负荷变化与筒体振动信号冲击成分强相关原理结合峰度参数ρ，构建的新型自适应VMD信号分解方法，该方法能够自适应确定VMD分解模态数，避免现有EMD和EEMD方法的模态混叠和端点效应。应用序列延拓对自相关函数进行改进，依据球磨机筒体振动信号模态分量频谱特征，综合考虑窗函数旁瓣峰值电平与衰减速率，构建Nuttall自卷积窗能量重心法，对球磨机筒体振动信号的IMF的进行频率校正，选取更加准确表征球磨机负荷变化的功率谱最大值对应频率作为特征量，能够有效提取磨矿过程中球磨机内部负荷特征，实现球磨机负荷状态的准确识别，为磨矿优化控制和效率提高提供准确、可靠的依据。

本发明的创新点体现在：自适应VMD，改进自相关函数以及构建新的功率谱估计方法，实现频率校正。

附图说明

图1为本发明的实施流程图；

图2为在不同频率偏差下，采用不同窗能量重心法对单频正弦信号校正误差率比较图；

图3为工业实测3种负荷状态球磨机筒体振动信号时域图；

图4为球磨机筒体振动信号进行EMD得到各IMF的峰度，并按从大到小顺序排序及其对应的峰度参数图；

图5为球磨机筒体重构信号VMD分解得到的IMF及其功率谱图；

图6为采用EMD方法提取球磨机筒体振动信号负荷特征量，导入支持向量机(Support Vector Machines，SVM)进行训练，得到的测试集预测结果图；

图7为EEMD方法提取球磨机筒体振动信号负荷特征量，导入SVM进行训练，得到的测试集预测结果图；

图8为采用自适应VMD方法提取球磨机筒体振动信号负荷特征量，导入SVM进行训练，得到的测试集预测结果图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：

为比较构建的Nuttall自卷积窗能量重心法的频率校正精度，采用采样频率f_s＝1024Hz，数据长度N＝4096的单频正弦信号，频率f在分辨率区间[49.5Hz，50.5Hz]等间隔取51个频率点，分别计算各常见余弦组合窗与Nuttall自卷积窗能量重心法在不同频偏下的归一化校正频率误差，仿真结果如图2所示。

为验证本发明提出方法的提取球磨机负荷特征的准确性和有效性，基于山东黄金集团焦家金矿选矿厂磨浮车间的溢流型球磨机(MQYΦ5.5×8.5m)进行实验验证。在球磨机筒体上安装加速度传感器，设置采样频率f_s为18750Hz，分别采集球磨机处于欠负荷、正常负荷、过负荷3种球磨机负荷状态的筒体振动信号，每个振动信号的序列长度为N＝164000，从3种负荷状态的振动信号中各提取50组样本共150组实验样本，3种负荷状态的时域波形如图3所示。

由图3可知，不同负荷状态的振动信号波形十分相似，均表现为中间大两头小的纺锤形状，且波峰位置也基本相同，纺锤形状中间的波峰较大是因为此时加速度传感器正处于球磨机最下方衬板上，球磨机内部钢球等物体抛落到该衬板，产生更大的加速度造成，总体规律性一致，但振动信号时域波形受到噪声干扰严重，只通过振动信号时域特征难以判断球磨机的内部负荷状态。

依据球磨机筒体振动信号生成原理结合峰度参数p滤除与球磨机负荷变化不相关的IMF，由峰度变化趋势选取峰度参数p大于85％所对应的K个IMF进行信号重构，图4为球磨机筒体振动信号进行EMI)得到各IMF的峰度按从大到小顺序排序及其对应的峰度参数图，由图4可确定前4个IMF的峰度参数以达到86.72％，因此选取进行信号重构的IMF以及数量K为4。

采用EMD方法、EEMD方法和本发明所提的自适应VMD和改进功率谱方法对3种负荷状态下的球磨机筒体振动信号提取特征，得到的功率谱最大值对应的频率{f_k}如表1所示，图5为球磨机筒体重构信号VMD分解得到的IMF及其功率谱图。

表1不同方法提取的球磨机负荷特征

对EMD方法、EEMD方法与本发明方法提取的特征量(同取各IMF功率谱最大值频率)进行比较研究，将各方法提取的特征向量分别导入SVM进行测试。每种负荷状态的50个样本数据随机分为两组，一组为训练集30个样本，另一组为测试集20个样本，再将每种状态的训练集与测试集分别构成总的训练集与测试集，即总训练集为90个样本，总测试集为60个样本。EMD方法、EEMD方法和本发明所提方法测试样本的识别结果分别如图6-8。

图6-8中测试样本类别1代表欠负荷、2代表正常负荷、3代表过负荷，基于发明方法提出的基于自适应VMD和改进功率谱估计的球磨机负荷特征提取方法测试集预测结果准确率达到91.6％，相比于EMD、EEMD方法具有更高的准确性。

虽然本发明已以较佳实施例公开如上，但实施例并不限定本发明。在不脱离本发明之精神和范围内，所做的任何等效变化或润饰，同样属于本发明之保护范围。因此本发明的保护范围应当以本申请的权利要求所界定的内容为标准。

Claims (1)

1.一种基于自适应VMD和改进功率谱的磨机负荷特征提取方法，其特征在于，包括以下步骤：

步骤1，对球磨机筒体振动信号x(i)进行EMD分解、筛选、重构以及VMD，得到模态分量{u(i)}；

步骤2：基于模态分量{u(i)}得到改进自相关函数{r(τ)}；

采用序列延拓，将模态分量{u(i)}进行改进自相关得到改进自相关函数{r(τ)}；

步骤3：对改进自相关函数{r(τ)}功率谱用Nuttall自卷积窗能量重心法进行频率校正，将校正后的功率谱最大值对应频率{f_k}作为表征球磨机负荷变化的特征量；

步骤1的过程为：

利用EMD对球磨机筒体振动信号x(i)进行分解得到IMF，依据球磨机负荷变化与筒体振动信号冲击成分强相关原理结合峰度参数ρ滤除与球磨机负荷变化不相关的IMF，由峰度变化趋势选取峰度参数ρ大于85％所对应的K个IMF进行信号重构并同时确定VMD的分解层数为K，对重构信号进行VMD，得到K个模态分量{u(i)}；

步骤1中，球磨机筒体振动信号为x(i)，i＝1,2,…n，n为信号序列长度；

利用EMD的快速自适应性能对球磨机筒体振动信号x(i)进行分解得到一系列模态函数{u_e(i)}，求取每个IMF的峰度{Ku}：

式中，

为u_e(i)的平均值；

将每个IMF的峰度{Ku}按从大到小的顺序排序，并计算峰度参数ρ：

式中，Kus为所有峰度的和，k_max为信号经EMD后得到的最大模态数；

选取峰度参数ρ大于85％所对应的K个IMF进行信号重构得到x_re(i)：

对重构信号x_re(i)进行VMD，得到K个模态分量{u(i)}；

步骤2中，对于离散序列u(i)，改进的自相关函数

为：

式中，τ为时延，其取值范围为[1,n]，U(i)的表达式为：

步骤3中，根据球磨机筒体振动信号模态分量频谱特征，综合考虑窗函数旁瓣峰值电平与衰减速率，采用Nuttall自卷积窗能量重心法，对改进自相关函数{r(τ)}功率谱进行频率校正，将校正后的功率谱最大值对应频率{f_k}作为表征球磨机负荷变化的特征量；

步骤3中，Nuttall自卷积窗能量重心法的校正式为：

式中，f_s为采样频率，n为采样点数，k₀为谐波信号经离散傅里叶变换得到的峰值谱线，此处谐波信号为改进自相关函数

P_N-2(k₀+i)为加窗功率谱第i条谱线功率值，P_N-2(ω)为Nuttall自卷积窗的功率谱函数，表达式如下：

式中，M为Nuttall窗函数的项数；b_m为Nuttall的窗系数，且满足

的约束条件，综合考虑旁瓣峰值电平和衰减速率，选用4项3阶Nuttall窗，即M＝4，b₀＝0.338946，b₁＝0.481973，b₂＝0.161054，b₃＝0.018027，W_R(·)为矩形窗的频谱函数，表达式为：

Images