CN113049684A - 基于变分模态分解的超声兰姆波缺陷信号识别方法 - Google Patents

Abstract

本发明涉及信息处理技术领域，提供一种基于变分模态分解的超声兰姆波缺陷信号识别方法，包括：步骤1，利用电磁超声换能器产生激发信号探测铝板，在铝板另一侧接收换能器接收超声兰姆波的原始信号；步骤2，确定变分模态分解的输入参数，所述输入参数有包括分解层数、二次惩罚因子和终止条件；步骤3，进行变分模态分解；步骤4，确定最终分解层数，得到分离后处在中心频率的子信号；步骤5，采用粒子群算法对变分模态分解输入参数‑惩罚因子进行选取和优化；步骤6，利用最终分解层数及优化后的惩罚因子进行变分模态，保证完整提取缺陷信号，并通过得到中心频率的缺陷信号进行缺陷识别。本发明能够提高超声兰姆波缺陷信号检测的准确率与识别率。

Description

基于变分模态分解的超声兰姆波缺陷信号识别方法

技术领域

本发明涉及信息处理技术领域，尤其涉及一种基于变分模态分解的超声兰姆波缺陷信号识别方法。

背景技术

由于超声检测信号中存在很多突发变量，在对其进行信号处理时需要根据超声信号自身特点选择适合的处理方法。实验中采集到的超声信号中往往存在着大量的噪声，这就会在对信号处理及特征分析时造成较大影响。在对各种缺陷进行评估时，并没有一套统一的标准，这就会导致很多检测信号即使含有缺陷也无法与其它相似的信号相互区别开来。

Dragomiretskiy提出了变分模态分解方法(VMD)，该方法与其他分解方法相比，不仅具有较好的理论基础，计算效率高，而且不易发生模态混叠与端点效应。目前，变分模态分解方法常被用来处理非线性和非平稳信号，如振动信号、超声检测信号等。

变分模态分解效果受分解个数K和变分函数中的惩罚参数α的影响。如何在分解前选择合适的参数组合是VMD是否能够有效分解的关键。Mohanty S等提出根据在选择不同分解层数K值时，分解后的固有模态分量的频率中心是否相互接近，来选定K-1分解个数的观测方法。但是，现有技术的变分模态分解中分解效果受分解层数k和惩罚因子α的影响，会存在模态混叠产生虚假分量现象以及无法准确检测到缺陷信号。

发明内容

本发明主要解决现有技术的变分模态分解中分解效果受分解层数k和惩罚因子α的影响，会存在模态混叠产生虚假分量现象以及无法准确检测到缺陷信号的技术问题，提出一种基于变分模态分解的超声兰姆波缺陷信号识别方法，以达到提高超声兰姆波缺陷信号检测的准确率与识别率的目的。

本发明提供了一种基于变分模态分解的超声兰姆波缺陷信号识别方法，包括以下过程：

步骤1，利用电磁超声换能器产生激发信号探测铝板，在铝板另一侧接收换能器接收超声兰姆波的原始信号；

步骤2，根据获得的超声兰姆波的原始信号，确定变分模态分解的输入参数，所述输入参数有包括分解层数、二次惩罚因子和终止条件；

步骤3，根据确定变分模态分解的输入参数，进行变分模态分解；

步骤4，利用固有模态函数计算K_临界-1个IMF分量与原始信号之间的互信息值I(X,Y)，并将该互信息值与互信息值阈值进行比较，将互信息值小于阈值β的IMF分量作为虚假分量，剔除虚假分量，确定最终的K值，进而确定最终分解层数为K_最终-1，其中，K_最终表示VMD分解的最终K值；同时，确定最终分解层数后，得到分离后处在中心频率的子信号；

步骤5，采用粒子群算法对变分模态分解输入参数-惩罚因子进行选取和优化，包括步骤501至步骤503：

步骤501，根据适应度值或适应度函数，将群中有限个随机粒子进行初始化，记录每个粒子的速度特性极值和位置特性极值；在每次迭代中，粒子会通过以下公式来更新的速度特性和位置特性：

v_i＝v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i) (10)

x_i＝x_i+v_i (11)

其中，i＝1,2,3,.....,N，N表示此群中粒子的总数；v_i表示粒子目前的速度特性；rand()表示介于(0,1)之间的随机数；x_i表示粒子当前的位置特性；c₁和c₂分别表示学习因子，通常情况下c₁＝c₂＝2；vi的最大值记为V_max，如果v_i大于V_max，则将vi＝V_max；

步骤502，找到速度特性极值和位置特性极值，引入动态惯性因子ω，达到最大迭代次数结束优化；引入惯性因子w后粒子的速度特性和位置特性公式变为：

v_i＝w×v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i) (12)

其中，w表示惯性因子，w≥0；

步骤503，根据步骤501、步骤502对变分模态分解方法中惩罚因子进行优化时，选取信噪比、平滑度、均方根误差三个评价指标为粒子群算法的适应度函数；通过对信号中这三个评价指标的计算，选取最适合变分模态分解方法的输入参数；

信噪比计算公式：

SNR＝20×log(norm(x)/norm(f-x)) (14)

平滑度计算公式：

均方根误差计算公式：

其中，f为原始信号，x为经处理后信号，N为信号样本长度；

步骤6，利用最终分解层数及优化后的惩罚因子进行变分模态，保证完整提取缺陷信号，并通过得到中心频率的缺陷信号进行缺陷识别。

进一步的，步骤2，包括步骤201-203：

步骤201，根据获得的超声兰姆波的原始信号，获得固有模态函数的瞬时频率；

步骤202，计算不同K值各瞬时频率均值的曲率，将曲率首次出现局部极大值处的K作为分解临界值，其中，K表示变分模态分解中的模态数；

步骤203，根据预设分解层数，确定变分模态分解的输入参数。

进一步的，步骤3，包括步骤301-305：

步骤301，通过希尔伯特变换得到每个模态函数的解析信号，得到单边频谱：

其中，t表示时间，u_k(t)表示模态函数，δ(t)表示冲激函数，j表示复数的虚部；

步骤302，对各模态解析信号预估中心频率，将每个模态的频谱调制到相应的基频带：

其中，

表示模态解析信号预估中心频率；

步骤303，计算以上解调信号梯度的平方范数，估计出各模态信号带宽，受约束的变分问题如下：

其中，{u_k}＝{u₁,...,u_k},{w_k}＝{w₁,...w_k}，u_k表示中心频率，w_k表示瞬时频率，

步骤304，根据二次惩罚因子α和拉格朗日乘法算子λ，将约束性变分问题变为非约束性变分问题，拉格朗日算子使得约束条件保持严格性，扩展的拉格朗日表达式为：

其中，α表示二次惩罚因子，λ表示拉格朗日乘法算子；

步骤305，变分模态分解中采用了乘法算子交替方向法通过不断交替更新

和λⁿ⁺¹寻求变分方程(5)中的极值点。

进一步的，步骤4，包括步骤401-403：

步骤401，计算各模态分量与原始信号的互信息值I(X,Y)：

其中，x表示信源发出消息；y表示接收到信号；p(x)表示先验概率；p(x|y)表示后验概率；

步骤402，对各模态分量与原始信号的互信息I(X,Y)做归一化处理：

β＝I(X,Y)/max(I(X,Y))

当分解后的IMF分量与原始信号的互信息值小于阈值时，分解出现了虚假分量，剔除该分量；

步骤403，剔除虚假分量后，确定最终的K值，进而确定最终分解层数为K_最终-1。

本发明提供的一种基于变分模态分解的超声兰姆波缺陷信号识别方法，引入K-瞬时频率均值方法、互信息值法与粒子群优化算法对变分模态分解过程中的主要参数分解层数、惩罚因子进行优化，避免了模态混叠与虚假分量的产生；对包含噪声的超声兰姆波信号进行分解，通过对超声信号进行处理去除超声信号中的噪声成分，得到准确的频率缺陷信号，继而可以对含缺陷的信号进行特征参数提取和有效识别，提高超声兰姆波缺陷信号检测的准确率与识别率。

附图说明

图1是本发明提供的基于变分模态分解的超声兰姆波缺陷信号识别方法的实现流程图。

图2为实施例一中的分解层数K选取流程图。

图3为实施例一的粒子群优化算法流程图。

图4为实施例二的实验设置流程图。

图5为实施例二的电磁超声采集实际信号。

图6为实施例二的K-瞬时频率均值关系图。

图7为实施例二的瞬时关系均值曲率图。

图8为实施例二的惩罚因子取不同值对应首波幅值。

图9为实施例二的惩罚因子取不同值时对应缺陷回波幅值。

具体实施方式

为使本发明解决的技术问题、采用的技术方案和达到的技术效果更加清楚，下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部内容。

实施例一

如图1所示，本发明实施例提供的基于变分模态分解的超声兰姆波缺陷信号识别方法，包括以下过程：

步骤1，利用电磁超声换能器产生激发信号探测铝板，在铝板另一侧接收换能器接收超声兰姆波的原始信号。

信号发生器产生激发信号，激发信号经功率放大器进行放大后，激发待探测铝板，待探测铝板中产生和传播超声兰姆波，超声兰姆波在经过铝板上缺陷时会发生反射和散射现象，通过铝板另一侧的接收换能器来接收回波信号f(t)，作为变分模态分解(VMD)需要处理超声兰姆波的原始信号f(t)。

其中，激发信号可以是调频信号，也可以是固定频率的正弦波信号。超声兰姆波的原始信号f(t)是经铝板传播后接收到的回波信号，包含了铝板中的缺陷信号、低频高频噪声及超声兰姆波在遇到铝板边界时产生的多次反射信号成分。

本发明利用变分模态分解方法(VMD)，变分模态分解方法是一种自适应、完全非递归的模态变分和信号处理的方法，其自适应性表现在根据实际情况确定所给序列的模态分解个数，随后的搜索和求解过程中可以自适应地匹配每种模态的最佳中心频率和有限带宽，并且可以实现固有模态分量(IMF)的有效分离、信号的频域划分、进而得到给定信号的有效分解成分，最终获得变分问题的最优解。与其他分解方法相比，不仅克服了EMD方法存在端点效应和模态分量混叠的问题，并且具有更坚实的数学理论基础。

步骤2，根据获得的超声兰姆波的原始信号f(t)，确定变分模态分解的输入参数。

在变分模态分解中，所述输入参数有：分解层数、二次惩罚因子α、终止条件ε；分解层数由变分模态分解中的模态数K的值决定，本发明中K值由瞬时频率均值来确定；二次惩罚因子α，当信号中含有强噪声时，为了达到良好的去噪效果，需使用二次惩罚因子；终止条件ε保证信号不会被VMD无休止分解下去，当ε＝1e-7时，则停止迭代。VMD的K值作为一个影响到最终分解效果的重要因素，K值过小，信号出现欠分解，发生模态混叠现象；K值过大，信号被过分解了，导致虚假分量的产生。因此本发明根据原始信号的瞬时频率均值与K之间的关系并对瞬时频率均值曲率进行量化分析来确定K值。输入参数的确定，如图2所示，具体过程为：

步骤201，根据获得的超声兰姆波的原始信号f(t)，获得固有模态函数(IMF)的瞬时频率w(t)。

超声兰姆波的原始信号受噪声、结构变化等一些因素的影响，具有非线性、非平稳特征，其频率具有时变性，因此，瞬时频率w(t)可以表征原始信号的局部特征。瞬时频率w(t)的计算公式如下：

其中，θ(t)表示瞬时相位，

f(t)表示超声兰姆波的原始信号，f(t)为时域内的连续信号，H[f(t)]表示f(t)解析信号的虚部。

步骤202，计算不同K值各瞬时频率均值的曲率，将曲率首次出现局部极大值处的K作为分解临界值。

在本步骤中，各瞬时频率均值的曲率首次出现极大值时，此时认为信号出现过分解。该极大值处的K值就是分解临界值，确定预设分解层数为K_临界-1，K_临界表示各瞬时频率均值的曲率首次出现极大值时的K值。

步骤203，根据预设分解层数，确定变分模态分解的输入参数。

步骤3，根据确定变分模态分解的输入参数，进行变分模态分解。

根据上述方法确定预设分解层数，对原始信号进行变分模态分解，得到固有模态函数(IMF)。其中，所述模态函数IMF都围绕各自的中心频率u_k扰动。变分模态分解的详细过程如下：

步骤301，通过希尔伯特变换(Hilbertransform)得到每个模态函数u_k(t)的解析信号，得到单边频谱：

其中，t表示时间，u_k(t)表示模态函数，δ(t)表示冲激函数，j表示复数的虚部.

步骤302，对各模态解析信号预估中心频率

将每个模态的频谱调制到相应的基频带：

步骤303，计算以上解调信号梯度的平方(L²)范数，估计出各模态信号带宽，受约束的变分问题如下：

其中，{u_k}＝{u₁,...,u_k},{w_k}＝{w₁,...w_k}，u_k表示中心频率，w_k表示瞬时频率，

步骤304，根据二次惩罚因子α和拉格朗日乘法算子λ，将约束性变分问题变为非约束性变分问题，其中二次惩罚因子可在高斯噪声存在的情况下保证信号的重构精度，拉格朗日算子使得约束条件保持严格性，扩展的拉格朗日表达式为：

步骤305，变分模态分解中采用了乘法算子交替方向法(ADMM)解决以上变分问题，通过交替更新

和λⁿ⁺¹寻求扩展式(5)中的极值点。

u_k最小化为：

w_k最小化为：

步骤306，首先初始化

和迭代次数n；之后，根据式(6)和式(7)更新u_k和w_k；再利用式子

更新λ；之后，对于给定的判别精度ε>0且ε<1e^-7，若

则停止迭代，否则再根据式(6)和式(7)输出最终的u_k和w_k。

步骤4，利用固有模态函数(IMF)计算K_临界-1个IMF分量与原始信号之间的互信息值I(X,Y)，并将该互信息值与互信息值阈值进行比较，将互信息值小于阈值β的IMF分量作为虚假分量，剔除虚假分量，确定最终的K值，进而确定最终分解层数为K_最终-1，其中，K_最终表示VMD分解的最终K值；同时，确定最终分解层数后，得到分离后处在中心频率的子信号，即选取所需频率的缺陷信号以备后续处理。剔除虚假分量的具体过程包括以下子步骤：

步骤401，计算各模态分量与原始信号的互信息值I(X,Y)：

其中，x表示信源发出消息；y表示接收到信号；p(x)表示先验概率；p(x|y)表示后验概率；

步骤402，对各模态分量与原始信号的互信息I(X,Y)做归一化处理：

β＝I(X,Y)/max(I(X,Y)) (9)

当分解后的IMF分量与原始信号的互信息值小于阈值时，即分解出现了虚假分量，剔除该分量。β表示互信息值阈值，可取为0.02。

步骤403，剔除虚假分量后，确定最终的K值，进而确定最终分解层数为K_最终-1。

步骤5，采用粒子群算法(Partical Swarm Optimization，PSO)对变分模态分解输入参数-惩罚因子α进行选取和优化。如图3所示，步骤5包括步骤501至步骤503：

步骤501，根据适应度值或适应度函数，将群中有限个随机粒子进行初始化，记录每个粒子的速度特性极值和位置特性极值，在每次迭代中，更新极值。

群中的所有粒子仅具有两个特性：即速度特性和位置特性。速度特性代表粒子寻优的快慢，位置特性代表粒子下一步寻优的方向。首先根据适应度值或适应度函数，将群中有限个随机粒子进行初始化，记录每个粒子的个体极值，找到第一个全局最优值。在接下来每一次的迭代中，粒子都会通过跟踪这两个速度特性极值和位置特性极值(pbest_i，gbest_i)来不断进行更新。在找到这两个最优值后，粒子会通过以下公式来更新自己的速度特性和位置特性：

v_i＝v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i) (10)

x_i＝x_i+v_i (11)

其中，i＝1,2,3,.....,N，N表示此群中粒子的总数；v_i表示粒子目前的速度特性；rand()表示介于(0,1)之间的随机数；x_i表示粒子当前的位置特性；c₁和c₂分别表示学习因子，通常情况下c₁＝c₂＝2；vi的最大值记为V_max，如果v_i大于V_max，则将vi＝V_max；

速度特性中的第一项表示上一次粒子的速度特性大小和方向；第二项是从当前粒子指向粒子个体极值的一个矢量值，表示粒子个体极值的求取过程；第三项是一个从当前粒子指向种群最优粒子的矢量值，表示粒子向全局最优解的运动过程。每个粒子都是通过这两部分来决定下一步的运动。

步骤502，找到速度特性极值和位置特性极值，引入动态惯性因子ω，达到最大迭代次数结束优化。

为了调节粒子在全局或局部上的寻优能力，将惯性因子w引入粒子群算法中。加入惯性因子w后粒子的速度特性和位置特性公式变为：

v_i＝w×v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i) (12)

其中，w表示惯性因子，w≥0。与固定值相比，动态的惯性因子w可以在问题求解过程中获得更好的寻优结果。w可以在粒子群算法搜寻最优解的过程中按照线性进行变化，也可以根据算法性能的某个测度函数来进行动态调整。

对惯性因子的权值进行线性递减是目前采用较多的方法之一，具体如下：

w^(t)＝(w_ini-w_end)(G_k-g)/G_k+w_end (13)

G_k表示设置的最大迭代次数；w_ini表示设置的初始惯性权值；w_end表示设置的迭代至最大次数时的惯性权值；

本实施例中惯性因子权值设置值：w_ini＝0.9，w_end＝0.4；

动态惯性因子w的引入，使得粒子群算法在对变分模态分解输入参数的优化过程中表现更加灵活，使分解后的效果有很大的提升。

步骤503，根据步骤501、步骤502对变分模态分解方法中惩罚因子进行优化时，本实施例选取信噪比、平滑度、均方根误差三个评价指标为粒子群算法的适应度函数。通过对信号中这三个参数的计算与评估来选取最适合变分模态分解方法的输入参数。

信噪比计算公式：

SNR＝20×log(norm(x)/norm(f-x)) (14)

平滑度计算公式：

均方根误差计算公式：

其中，f为原始信号，x为经处理后信号，N为信号样本长度；

将联合特征参数L作为粒子群算法中的适应度函数来对群中的随机粒子进行初始化操作。通过分别对粒子的两个极值不断迭代更新来寻找群中适应度函数在最优解时的惩罚因子。

联合特征参数L：

L＝a×SNR+b×R+c×RMSE (17)

a,b,c分别为信号的信噪比，平滑度与均方根误差参数的权值系数，取值范围(0，1)。

设β＝[a,b,c]，X＝[SNR,R,RSME]^T

故有：

L＝βX (18)

通过求取不同惩罚因子时，联合特征参数L的最大值来确定惩罚因子的最优选取值。

引入残差平方和函数S：

S(β)＝||Xβ-L||² (19)

当

时，S(β)取最小值，记作：

对S(β)进行微分求最值，有：

如果矩阵X^TX非奇异，则β有唯一解：

β求解出来即确定L的各项参数系数矩阵，将L作为适应度函数求取个体与整体最优值并带入公式(10)，(11)进行迭代求解。

步骤6，利用最终分解层数及优化后的惩罚因子进行变分模态，保证完整提取缺陷信号，并通过得到中心频率的缺陷信号进行缺陷识别。

将与原始信号互信息值大于阈值的模态分量提取出来，得到处在中心频率的缺陷信号，分析分解后得到的缺陷信号时域图，可以对缺陷信号进行识别。通过计算各兰姆波的对称模态或反对称模态在铝板中的传播速度与时间的乘积得到缺陷与发射探头之间的距离，从而可以确定缺陷的具体位置。

实施例二

本发明实施例提供的基于变分模态分解的超声兰姆波缺陷信号识别方法，包括：

步骤1，利用信号发生器产生激发信号探测铝板，在铝板另一侧接收换能器接收超声兰姆波的原始信号。

图4为检测板材为铝板材料，具体参数为：1000mm×900mm×0.5mm。激发频率为130KHz，A0为主要激发模态，群速度为1507m/s，相速度为784.8m/s。实验采用的信号由超声换能器激发，信号在1mm铝板表面及内部传播并且遇到边界和缺陷时发生反射与散射现象，最终形成兰姆波。实际信号如图5所示，信号由接收换能器采集到的超声兰姆信号与实验噪声构成，采样点为1400个。

步骤2，根据获得的超声兰姆波的原始信号f(t)，确定变分模态分解的输入参数。

VMD的分解层数K作为一个影响到最终分解效果的重要因素，K值过小，信号出现欠分解，发生模态混叠现象；K值过大，信号被过分解了，导致虚假分量的产生。

步骤3，根据确定变分模态分解的输入参数，进行变分模态分解。

绘制分解层数与信号分解后各IMF分量的瞬时频率均值之间的关系如图6。从实际信号瞬时频率均值的曲率图可以看出，当K＝4时，瞬时频率均值开始出现骤减，此时信号被过度分解，临界分解层为4层，因此分解层数K值为3。

步骤4，利用固有模态函数(IMF)计算K_临界-1个IMF分量与原始信号之间的互信息值I(X,Y)，并将该互信息值与互信息值阈值进行比较，将互信息值小于阈值β的IMF分量作为虚假分量，剔除虚假分量，确定最终的K值，进而确定最终分解层数为K_最终-1，其中，K_最终表示VMD分解的最终K值。不同K值下各固有模态函数(IMF)与原始信号互信息值如下表1所示：

表1各IMF与实际信号互信息值表

从表中可以看出，当K＝3时，各IMF分量与实际信号之间的互信息值均高于阈值，不存在虚假分量；当K＝4时，IMF4分量与实际信号之间的互信息值低于0.02，此时认为IMF4分量为虚假分量,即出现了过度分解现象，IMF4分量需要被剔除；同样的情况，当K＝5时，IMF分量3与IMF4分量的互信息值均小于阈值，因此被剔除。从IMF分量与实际信号之间的互信息值方法可以得知VMD对实际信号的分解准确层数为3层，与上述K-瞬时频率均值法确定的分解层数一致。

步骤5，采用粒子群算法对变分模态分解输入参数-惩罚因子α进行选取和优化。

将联合特征参数L作为粒子群算法中的适应度函数来对群中的随机粒子进行初始化操作。通过分别对粒子的两个极值不断迭代更新来寻找群中适应度函数在最优解时的惩罚因子。

通过求解计算得到分解层数为3时，惩罚因子α＝2000左右时，适应度(联合特征参数)L最大。此时，信号信噪比与平滑度较高，均方根误差也相对较小，且α＝2000左右的信号更加稳定，能够获得最优的分解效果。

惩罚因子α取不同值时候，根据变分模态分解的信号首波幅值与缺陷回波幅值大小，绘制表2，将惩罚因子α＝2000时的首幅值与缺陷回波幅值作为基准并与其他惩罚因子做比较，计算误差值，具体如下：

表2惩罚因子α对应信号首波与缺陷回波幅值表

在不同惩罚因子α下，变分模态分解后信号首波幅值如图8。在不同惩罚因子α下，变分模态分解后信号缺陷回波幅值如图9。从图8和图9可以看出，当惩罚因子α在200-1200之间，变分模态分解之后的信号首波幅值和缺陷回波幅值都处于不稳定状态，这是由于惩罚因子α过小，信号还未完全分解；惩罚因子α在1400-2400之间，信号首波幅值与缺陷回波幅值较为平稳；惩罚因子α取值大于2400时，信号的首波幅值与缺陷回波幅值出现大范围波动且没有规律(α值大于3000时信号首波幅值与缺陷回波幅值未在表格与图中标出，结果与2400-3000结论一致)，这是由于惩罚因子α越大，信号在中心频率u_k两侧的衰减就越快，容易导致模态混叠现象。

因此，在根据实验中的超声信号选择变分模态分解参数时，将分解层数K设置为3，惩罚因子α设置为2000。经上述优化方法得到VMD分解层数为3，惩罚因子大小为2000。

步骤6，利用最终分解层数及优化后的惩罚因子进行变分模态，保证完整提取缺陷信号，并通过得到中心频率的缺陷信号进行缺陷识别。

从分解结果可以看出，VMD不仅可以有效地去除虚假分量，而且每个IMF都在固有频率周围一定带宽扰动，没有模态混叠现象的产生。该方法实现了对超声无损检测信号的多尺度表征，具有较好的噪声鲁棒性。

本发明实施例提供的基于变分模态分解的超声兰姆波缺陷信号识别方法，在变分模态分解中通过瞬时频率均值来确定分解层数，利用粒子群算法对输入参数-惩罚因子优化，计算各固有模态函数与原始信号之间的互信息值来剔除小于阈值的虚假分量，最终得到在特定频率的缺陷信号，实现对缺陷信号的准确识别。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (4)

1.一种基于变分模态分解的超声兰姆波缺陷信号识别方法，其特征在于，包括以下过程：

步骤1，利用电磁超声换能器产生激发信号探测铝板，在铝板另一侧接收换能器接收超声兰姆波的原始信号；

步骤2，根据获得的超声兰姆波的原始信号，确定变分模态分解的输入参数，所述输入参数有包括分解层数、二次惩罚因子和终止条件；

步骤3，根据确定变分模态分解的输入参数，进行变分模态分解；

步骤4，利用固有模态函数计算K_临界-1个IMF分量与原始信号之间的互信息值I(X,Y)，并将该互信息值与互信息值阈值进行比较，将互信息值小于阈值β的IMF分量作为虚假分量，剔除虚假分量，确定最终的K值，进而确定最终分解层数为K_最终-1，其中，K_最终表示VMD分解的最终K值；同时，确定最终分解层数后，得到分离后处在中心频率的子信号；

步骤5，采用粒子群算法对变分模态分解输入参数-惩罚因子进行选取和优化，包括步骤501至步骤503：

步骤501，根据适应度值或适应度函数，将群中有限个随机粒子进行初始化，记录每个粒子的速度特性极值和位置特性极值；在每次迭代中，粒子会通过以下公式来更新的速度特性和位置特性：

v_i＝v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i) (10)

x_i＝x_i+v_i (11)

其中，i＝1,2,3,.....,N，N表示此群中粒子的总数；v_i表示粒子目前的速度特性；rand()表示介于(0,1)之间的随机数；x_i表示粒子当前的位置特性；c₁和c₂分别表示学习因子，通常情况下c₁＝c₂＝2；v_i的最大值记为V_max，如果v_i大于V_max，则将v_i＝V_max；

步骤502，找到速度特性极值和位置特性极值，引入动态惯性因子ω，达到最大迭代次数结束优化；引入惯性因子w后粒子的速度特性和位置特性公式变为：

v_i＝w×v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i) (12)

其中，w表示惯性因子，w≥0；

步骤503，根据步骤501、步骤502对变分模态分解方法中惩罚因子进行优化时，选取信噪比、平滑度、均方根误差三个评价指标为粒子群算法的适应度函数；通过对信号中这三个评价指标的计算，选取最适合变分模态分解方法的输入参数；

信噪比计算公式：

SNR＝20×log(norm(x)/norm(f-x)) (14)

平滑度计算公式：

均方根误差计算公式：

其中，f为原始信号，x为经处理后信号，N为信号样本长度；

步骤6，利用最终分解层数及优化后的惩罚因子进行变分模态，保证完整提取缺陷信号，并通过得到中心频率的缺陷信号进行缺陷识别。

2.根据权利要求1所述的基于变分模态分解的超声兰姆波缺陷信号识别方法，其特征在于，步骤2，包括步骤201-203：

步骤201，根据获得的超声兰姆波的原始信号，获得固有模态函数的瞬时频率；

步骤202，计算不同K值各瞬时频率均值的曲率，将曲率首次出现局部极大值处的K作为分解临界值，其中，K表示变分模态分解中的模态数；

步骤203，根据预设分解层数，确定变分模态分解的输入参数。

3.根据权利要求2所述的基于变分模态分解的超声兰姆波缺陷信号识别方法，其特征在于，步骤3，包括步骤301-305：

步骤301，通过希尔伯特变换得到每个模态函数的解析信号，得到单边频谱：

其中，t表示时间，u_k(t)表示模态函数，δ(t)表示冲激函数，j表示复数的虚部；

步骤302，对各模态解析信号预估中心频率，将每个模态的频谱调制到相应的基频带：

其中，

表示模态解析信号预估中心频率；

步骤303，计算以上解调信号梯度的平方范数，估计出各模态信号带宽，受约束的变分问题如下：

其中，{u_k}＝{u₁,...,u_k},{w_k}＝{w₁,...w_k}，u_k表示中心频率，w_k表示瞬时频率，

步骤304，根据二次惩罚因子α和拉格朗日乘法算子λ，将约束性变分问题变为非约束性变分问题，拉格朗日算子使得约束条件保持严格性，扩展的拉格朗日表达式为：

其中，α表示二次惩罚因子，λ表示拉格朗日乘法算子；

步骤305，变分模态分解中采用了乘法算子交替方向法，通过不断交替更新

和λⁿ⁺¹寻求变分方程(5)中的极值点。

4.根据权利要求3所述的基于变分模态分解的超声兰姆波缺陷信号识别方法，其特征在于，步骤4，包括步骤401-403：

步骤401，计算各模态分量与原始信号的互信息值I(X,Y)：

其中，x表示信源发出消息；y表示接收到信号；p(x)表示先验概率；p(x|y)表示后验概率；

步骤402，对各模态分量与原始信号的互信息I(X,Y)做归一化处理：

β＝I(X,Y)/max(I(X,Y))

当分解后的IMF分量与原始信号的互信息值小于阈值时，分解出现了虚假分量，剔除该分量；

步骤403，剔除虚假分量后，确定最终的K值，进而确定最终分解层数为K_最终-1。

Images