基于优化VMD组合降噪和LMD的水轮机空化声发射信号特征提取

刘忠，潘宜桦，邹淑云，陈星宇，李志鹏

（长沙理工大学能源与动力工程学院，湖南 长沙 410114）

空化是流体中局部压力小于其饱和蒸气压而产生的局部流体的复杂流动现象[1]，是造成水轮机出力和效率降低、过流部件表面材料破坏并产生剧烈振动和噪声的重要原因[2].空化发生时伴有强烈的冲击波和微射流[3]，学者们对水轮机空化声发射（acoustic emission，AE）信号进行了诸多研究.RUS等[4]探究了AE信号的叶片流道调制水平相对值随空化系数变化的规律.刘忠等[5]以灰狼与布谷鸟混合算法进行参数寻优，建立本征模态函数（intrinsic mode function，IMF）能量与空化系数的变化关系.刘忠等[6]采用改进经验模态分解方法（empirical mode decomposition，EMD）提取各IMF的关联维数,分析关联维数随水轮机空化系数的变化关系.

局部均值分解（local mean decomposition，LMD）[7]是一种新的时频分析方法，它能将信号自适应地分解成多个有物理意义的乘积函数（product function，PF）分量之和.较之EMD方法，LMD迭代次数少，能更好地抑制端点效应，保留信号信息，故广泛应用于机械设备的故障诊断领域.由于AE信号易受外界因素影响[8]，采集到的水轮机空化AE信号中，常带有噪声成分.为提高特征提取准确度，须对信号进行降噪处理.刘忠等[9]采用EMD阈值对空化声发射信号进行降噪，证明了该方法降噪效果优于小波及小波包阈值降噪.EMD能有效分解信号，但存在模态混叠、端点效应.

变分模态分解（variational mode decomposition，VMD）[10]是一种基于经典维纳滤波的自适应方法，且精度高、噪声鲁棒性强.避免了EMD递归的束缚，可将其运用于信号降噪.VMD中存在分解层数K、惩罚因子α不确定的问题，影响分解效果.唐贵基等[11]采用中心频率观察法对K进行选取，但α人为设定.冉茂霞等[12]将粒子群算法优化VMD的方法运用到磁瓦声振信号分析中.但粒子群算法易陷入局部最优.哈里斯鹰优化算法（HHO）是灵感来源于鹰的捕食行为的一种仿生算法[13]，具有收敛性强、设置参数少和易实现等优点.文中拟采用HHO对VMD参数进行全局搜寻，得到最优的VMD参数组合，形成优化VMD方法，并结合Birge-Massart策略[14]应用于水轮机空化AE信号的降噪中.

1 基本理论

1.1 VMD算法

VMD算法的实质是构造变分问题，并求解变分模型.VMD算法将信号分解为K个模态函数uk（t），uk（t）的中心频率ωk固定，使得每个uk（t）的估计带宽之和最小.约束变分模型为

（1）

式中：δ（t）为冲击函数.

由于拉格朗日乘子λ（t）、二次惩罚因子α的引入，约束变分问题变为非约束变分问题.α作为保证重构信号保真度的一个经典方法，具有很好的收敛性.即使存在噪声，也可使信号的重构精度得到保障.λ（t）使约束条件严格性得到保障.拉格朗日扩展表达式为

L（{uk},{ωk},λ）=

（2）

（3）

VMD 算法在频域内不断更新并利用傅立叶逆变换得到模态分量uk和ω.具体步骤如下：

2） 根据式（3）在频域内更新uk和ωk.

3） 更新拉格朗日乘子λ：

1.2 LMD算法

LMD可将非线性非平稳信号分解为一组单分量调幅-调频信号.算法的实现过程如下：

1） 给定x（t），找出x（t）的极值点ni，求相邻的局部极值点ni，ni+1的均值mi，求其包络估计值ai.mi和ai经平滑操作，得到局部均值函数m11（t）、包络估计函数a11（t）.

（4）

2） 将m11（t）从x（t）中分离，即

h11（t）=x（t）-m11（t）.

（5）

3） 将h11（t）解调，即

（6）

4） 重复上述迭代过程：

（7）

其迭代终止条件为

（8）

5） 将步骤（4）中得到的全部a1q（t）做乘积，得到包络信号为

（9）

6） 将纯调频信号s1n（t）与a1（t）相乘，得到的即为x（t）的第1个PF分量：

PF1=a1（t）s1n（t）.

（10）

7） 从x（t）中分离PF1，余下的成分构建成u1（t）.再将u1（t）作为初始信号，重复上述过程，u（k）为单调函数则停止.

（11）

经过LMD处理，x（t）被分解为一系列PF分量和残余分量u（k），即：

（12）

1.3 HHO理论

HHO算法主要由以下3个阶段组成.

1）探索阶段：哈里斯鹰随机栖息在某地.该阶段如式（13）所示.

（13）

式中：X（t）为鹰的位置；Xr（t）为兔子的位置；q及r1～r4∈[0,1]，且为随机数，由q来决定需采用什么策略；Ub为搜索范围的上界；Lb为下界；Xm（t）为鹰的平均位置；G为鹰的总数.

2） 过渡阶段：兔子在逃跑时，其能量会大大地降低.定义兔子的逃逸能量E为

（14）

3） 开发阶段：该阶段根据不同的条件，分为4个不同的围攻方式.

① 当0.5<|E|≤1.0和r≥0.5时，采用软围攻策略，其位置更新表达式为

X（t+1）=ΔX（t）-E|JXrabbit（t）-X（t）|,

（15）

式中：ΔX（t）=Xrabbit（t）-X（t）,为兔子的位置向量和当前位置的差值；J为0～2的随机数.

② 当|E|<0.5且r>0.5时，用硬围攻策略更新位置：

X（t+1）=Xrabbit（t）-E|ΔX（t）|.

（16）

③ 当0.5≤|E|<1.0和r<0.5时，采用渐进式快速俯冲的软包围更新位置为

（17）

式中：Y=Xrabbit（t）-E|JXrabbit（t）-X（t）|；Z=Y+S×LF（D），其中D为求解问题的维度，S为随机向量，LF为莱维飞行的表达式.

④ 当|E|<0.5且r<0.5时，用渐进式快速俯冲的硬包围更新位置为

（18）

式中：Y=Xrabbit（t）-E|JXrabbit（t）-Xm（t）|.

2 HHO优化VMD组合降噪

采用HHO对VMD进行组合参数寻优.具体步骤如下：

1） 初始化参数.迭代次数T，群体个数G，优化参数个数dim，搜索范围上下界Ub，Lb.

2） 设置适应度函数.计算经VMD分解后，各IMF的散布熵[15][E1,E2,…,EK]，及各IMF的平均散布熵Ea.定义散布熵差异系数C为各IMF的散布熵与平均散布熵的差方和，即

（19）

相关系数可表述变量间的近似程度，其表达式为

（20）

式中：imfi为经VMD分解后，得到的第i个分量;r*为imfi与原始信号的相关系数.

C越大，则分解效果越好；r*越大，则相关程度越高.为得到最大化差异的IMF，且使原始信号尽可能多地保留信息，定义散布熵差异相关系数最小值Q，作为适应度函数，其表达式为

Q=-max （C·r）,

（21）

利用HHO对VMD进行寻优.

3） 循环迭代至终止，输出野兔的位置，其对应的K和α即为VMD最佳参数组合.

4） 通过HHO优化VMD算法，得到一系列IMF分量，计算各分量的相关系数，对相关系数小于0.1的IMF进行剔除，大于0.5的保留，0.1～0.5的IMF采用Birge-Massart准则进行降噪，并与保留的分量重构.降噪方法的流程图如图1所示.

3 仿真信号分析

构造高频仿真信号以检验HHO优化VMD与Birge-Massart策略组合降噪的可行性.构造纯净信号x0（t），x0（t）如式（22）所示.向纯净信号x0（t）加入10 dB的高斯白噪声，得到含噪信号x（t），2种信号的时域图如图2所示，图中A为振幅，N为点数.

（22）

HHO优化VMD算法中，分解层数K设为3～15，惩罚因子设为200～3 000.经迭代寻优，获得适应度函数的收敛曲线，如图3所示，图中C为迭代次数.得到全局最优的适应度值为-0.44，最佳K和α组合为[10,232].

计算各IMF的相关系数r，如表1所示.对大于0.500 0的分量保留，小于0.1的剔除.0.100 0～0.500 0的进行Birge-Massart策略降噪，并与保留的分量重构，获得降噪后的信号.从表1中可看出IMF1—IMF6的相关系数极小，可以直接剔除，IMF7,IMF8,IMF10的相关系数小于0.500 0，则对其进行Birge-Massart策略降噪处理.降噪后，与相关系数大于0.5的IMF9进行重构，得到经优化VMD组合降噪后的信号.

表1 IMF的相关系数

为证明该方法的有效性，与传统VMD降噪及Birge-Massart策略降噪进行对比，降噪后，时域对比图如图4所示.从图中可以看出，每种降噪方法都能在一定程度上降低部分噪声，具有一定的降噪效果.通过3种方法进行处理后，得到的波形图区别不是很大，但通过比较信号幅值与整体趋势，可观察到文中所提方法的效果最佳，且与原始信号最接近.

为了对降噪效果做进一步的评定，除了通过观察波形图，还可以对降噪效果做定量分析，即可采用常用的降噪效果评定指标对其进行评定，信噪比SNR与均方差MSE是常用的指标，其公式如式（23）所示.SNR越大，MSE越小，表明降噪效果越好.3种方法的降噪效果如表2所示.

（23）

式中：x0（n）为初始信号；s（n）为降噪信号.

表2 3种方法的降噪效果对比

从图4及表2中可以看出，文中所提到的降噪方法得到的波形更加平滑，同时更接近原始信号，且信噪比最大，均方根最小，从而证明该方法的可行性.

4 试验信号分析

混流式水轮机空化试验信号来自一座处于国内领先水平、综合精度<±0.2%的水轮机模型试验台.该试验台实时监测模型水轮机的相关物理量.水轮机尾水管段为透明有机玻璃，可用闪频仪记录水轮机流态的变化.采用PCI-9846H多功能数据采集卡进行信号采集.文中使用的AE信号数据截取自水轮机导叶拐臂上不同空化状态对应的一段.

采用优化VMD组合降噪方法分别对不同空化状态的信号进行降噪.限于篇幅，以无空化为例，对比原始信号和降噪后的信号的时域图及频谱图，如图5所示.

从图5可看出，试验信号的频率范围为0～500 kHz.经优化VMD组合降噪后主要集中分布于200 kHz以下，去除了200 kHz以上的高频噪声，并保存了有用成分，为后续的特征提取研究提供了比较好的试验信号.

不同空化状态下所采集到的AE信号，经过降噪处理，并通过LMD分解所得到的PF分量不同，数量也可能不一样，对应各频带能量也有所差异.因此，将能量作为特征参数来描述空化的变化情况.设能量Ei为

（24）

式中：xi为LMD分解所得到的第i个PF分量；k为LMD的分解层数.

计算各PF分量的能量，保留相关程度较大的分量，绘制能量与空化系数的关系变化曲线.图6为PF能量与空化系数图，图中Ea为能量值;θ为空化系数.

从图6可以看出，PF1和PF2的能量随空化系数呈现明显的负相关.这是由于在转轮叶片背面，压力逐渐降低，且未达到饱和汽化压力，未产生空化.当压力低至饱和汽化压力，则开始有空泡产生，此时为空化初生状态，AE信号能量有所增多.伴随空泡数的变多，其溃灭程度加剧，AE信号能量变得更多.因此，AE信号经LMD分解后，求得的主要PF的能量随空化系数变化的情况，可以反映空化状态的变化.

5 结 论

1） 引入了优化VMD组合降噪方法对空化AE信号进行处理，与传统VMD降噪、Birge-Massart降噪进行对比，表明优化VMD组合降噪方法更具优势.

2） 利用LMD算法对降噪后的水轮机空化AE信号进行分解，分析结果表明，PF分量的能量与空化系数之间呈现负相关，可用于水轮机空化状态的动态识别.