改进的变分模态分解在锚杆和锚固质量检测中的应用

闫志勋，刘怡明2，林 童，孙晓云，王明明

(1. 石家庄铁道大学 电气与电子工程学院，河北 石家庄 050043；2. 苏州大学 机电工程学院，江苏 苏州 215137)

作为一种新兴的超声导波检测手段，磁致伸缩导波技术在锚杆和锚固质量无损检测领域中有检测快速、可靠性高、检测距离长等优点，在各工程结构的健康检测和监测中有良好的应用前景[1-4]。由于电磁超声信号幅值较微弱，常被噪声所淹没，使得采集到的信号失真[5]，因此，在磁致伸缩锚杆检测中，对信号进行降噪十分必要。

在磁致伸缩锚杆检测中，检测信号因噪声的叠加而被视为非稳态信号。针对此类信号，所采用的降噪方法主要有经验模态分解(EMD)、小波阈值法等，但是小波阈值法的小波基和阈值选取较困难，EMD方法存在模态混叠与端点效应等问题，并且抗噪性较差。文献[6]中提出变分模态分解(VMD)方法，该方法相对于递归式模态分解方法有完善的数学理论基础，能更准确地分解信号，运算效率高，抗噪性好，并且不易发生模态混叠和端点效应。该方法需要预先设置惩罚因子和模态个数，分解结果在较大程度上受这2个参数的影响，然而磁致伸缩锚杆检测信号是复杂的，模态个数和惩罚因子难以确定。

为了实现信号的最优分解，针对模态个数和惩罚因子难以确定的问题，本文中利用布谷鸟搜索(CS)算法对VMD的最优参数组合进行搜寻[7]，从而确定VMD的最优参数，实现信号的最优分解。单独使用VMD算法对信号进行处理可能导致不能完全剥离反射波信号中的回波信息。独立分量分析(ICA)是一种基于样本高阶统计信息的新的信号处理方法[8]，目前的研究大都是基于多观测通道的，而通常采集到的锚杆磁致伸缩信号均为一维数据，不能满足ICA对于多维数据的要求。

针对以上问题，本文中提出CS算法优化VMD与ICA结合的方法，以达到更好的降噪效果。

1 CS算法优化VMD参数

1.1 CS算法与VMD基本理论

1.1.1 CS算法

CS算法是Yang等[7]基于布谷鸟的巢寄生特性及莱维飞行而提出的，首先假设以下3条理想化规则：

1)每只布谷鸟一次只产一个卵，孵化时的鸟巢是随机选择；

2)最好的巢保留到下一代；

3)寄生巢数量是固定的，寄生宿主发现外来鸟蛋的概率也是固定的，设为Pa∈[0,1]。

基于以上理想化规则，布谷鸟对鸟巢的位置和路径更新是按照莱维飞行模式进行的，路径和位置的更新公式为

(1)

1.1.2 VMD算法

VMD算法是对变分问题的求解，分解步骤如下。

1)构建目标函数。

①进行Hilbert变换，寻求每个模态函数uk(t)解析信号，并得到单边频谱

(2)

式中：δ(t)为冲击函数；k为分解的本征模态分量(IMF)的层数。

②通过加入指数项，估计每个模态的中心频率e-jωkt，其中ωk为中心频率，把每个模态的频谱调制到相应的基频带，即

(3)

③估计每个模态信号带宽，从而获得由变分问题组成的目标函数，即

min{uk},{ωk}∑Kk=1∂∂tδ(t)+jπt[]uk(t){}e-jωkt22{},∑Kk=1uk=f(t),ìîíïïïï(4)

2)引入拉格朗日算子λ(t)和二次惩罚因子α，得到增广拉格朗日函数

L({uk},{ωk},λ)=

(5)

式中：〈·〉为内积；K为当前分层数。

(6)

式中arg min表示式子达到最小值时的变量的取值。

利用Parseval-Plancherel傅里叶等距变换后，将频率ω替代为ω-ωk，并将得到的结果转换为负频率区间积分的形式，则优化问题的解为

(7)

根据同样的过程，求得中心频率更新方式

(8)

(9)

综上所述，VMD的实现过程如下。

2)根据式(7)—(9)，在频域内更新uk、ωk、λ；

1.1.3 峭度准则

峭度是描述波形尖峰度的一个无量纲参数，对冲击信号特别敏感。峭度值C的定义为

(10)

式中：μ、σ分别为信号x的均值和标准差[9]。当峭度值C=3时，定义分布曲线具有正常峰度，即零峭度，如果信号中叠加了大量冲击脉冲信号，则信号的峭度值会明显增大。利用VMD算法对锚杆的磁致伸缩信号进行分解后，如果IMF中更多的是反射波信息，波形中出现规律性冲击脉冲，则分量信号峭度值较大。

1.2 CS-VMD参数优化

现实中磁致伸缩信号是复杂的，在进行VMD运算时，核心问题是如何选取合适的K及α。如果仅改变一个参数，另一个参数设置为常数，这样得到的只是局部最优结果[10]。为了得到一组最佳参数，本文中利用CS算法对VMD的输入参数K和α进行寻优，在进行参数优化时，选用峭度值作为适应度函数。

在VMD运算后，计算所有IMF的峭度值，将峭度值的最大值称为局部极大值，用maxEIMF表示，与该值所对应的分量称为局部最佳分量。

VMD算法的α和K的优化步骤如下：

1)初始化鸟巢个数、更新概率、迭代次数等参数，初始化群体，随机产生J个鸟巢的初始位置yi(i=1,2,…,J)，确定适应度函数为峭度值；

2)在不同参数条件下对信号进行VMD运算，计算每组参数相应的适应度值maxEIMF并记录当前的最优解；

3)保留上一代最佳IMF的鸟巢，利用式(1)对鸟巢位置进行更新；

4)比较新的鸟巢适应度值与上一代鸟巢适应度值，选取更好的鸟巢位置并保留；

5)比较随机数r和发现概率Pa，若r>Pa，则更新鸟巢位置，反之保留位置不变；

6)若未满足终止条件，则返回步骤2)；

7)输出全局最优位置。

算法流程如图1所示。

图1 布谷鸟搜索(CS)算法优化变分模态分解(VMD)算法流程图

2 VMD与ICA联合降噪方法

2.1 ICA

ICA是在源信号的先验信息和信道传输参数均未知的情况下，利用随机矩阵对信号进行混合，对混合后的信号进行线性变换，使得变换后的信号尽量逼近源信号。设Y(t)=(Y1(t),Y2(t),…,YM(t))T为M个观测信号，设源信号S(t)=(S1(t),S2(t),…,SN(t))T为N个未知独立源信号，各分量之间互相独立，观测到的线性混合模型为

Y(t)=AS(t)，

(11)

式中A为满秩的M×N型混合矩阵(M≥N)。ICA算法的目的是找到分离矩阵W，使得

(12)

在混合矩阵A及源信号S都未知的情况下，估计出的各独立分量的排列顺序及幅值具有不确定性，但是不影响信号整体的判断。为了提高计算的收敛速率，通常采用快速ICA(FastICA)算法。

2.2 排列熵算法

排列熵(permutation entropy,PE)是由Bandt等[11]提出的一种衡量序列复杂度和动力学突变的方法。计算步骤如下。

给定时间序列{X(i),i=1,2,…,n}进行相空间重构，可得到重构矩阵[12]为

j=1,2,…,Q,

(13)

式中：m为嵌入维数；T为延迟时间；Q为重构分量个数。对X(i)中第j个向量升序排列，得

x[i+(j1-1)T]≤x[i+(j2-1)T]≤…≤

x[i+(jm-1)T]。

(14)

如果重构向量存在相等的情况，则按原顺序排列，因此，对于任意一个时间序列X(j)，符号序列为S(l)=(j1,j2,…,jm)其中l=1,2,…,q，且q≤m!。每种符号序列出现概率记为P1,P2,…,Pq，此时，时间序列X(i)的q种不同符号序列的PE可定义为

(15)

当Pj=1/m!时，Hp(m)最大，对应值为ln(m!)。对Hp(m)进行归一化处理，即

0≤Hp=Hp/ln(m!)≤1。

(16)

Hp值表示时间序列的复杂度，其值越小，信号越规则，反之，随机性越强。

2.3 VMD-ICA联合降噪

本文中通过CS-VMD运算，得到若干IMF，选取保留较多振动信号的分量重构成一组信号，利用FastICA算法对其与源信号组成的2路输入通道进行解混，得到降噪后的信号。对于由噪声引起的毛刺，可以进一步进行小波阈值处理。

选取IMF重构信号时，选用PE作为性能指标，PE对突变的信号有有效的检测效果[13]，因此，根据PE值进行信号的重构。

综上所述，提出基于CS-VMD参数优化与ICA联合降噪方法，实现步骤如下：

1)利用CS算法优化VMD参数;

2)利用VMD方法分解回波信号;

3)依据PE值对分解的各模态进行重构，对单通道信号进行升维;

4)利用FastICA算法解混升维后的信号;

5)对解混后的信号进行小波阈值处理，滤除由噪声引起的毛刺。

算法流程如图2所示。

3 裸锚杆长度实验分析

测试中采用的锚杆为长度为3 m，直径为22 mm。锚杆磁致伸缩信号采集示意图如图3所示。检测方法如下：通过激励线圈产生磁致伸缩导波，当导波到达端部时会发生反射，对应的反射波会有相位改变。由安放在锚杆另外一端的感应线圈接收反射回来的信号，再通过采集模块传送至上位机。假设锚杆的长度为L0，感应线圈接收到端面反射信号的时间为t，超声导波在锚杆中的传播波速为c，通常取值为5 000 m/s。计算锚杆长度的公式为

(17)

CS—布谷鸟搜索；IMF—本征模态分量；FastICA—快速独立分量分析。图2 变分模态分解-独立分量分析(VMD-ICA)联合降噪算法流程图

为了验证本文中提出方法的降噪效果，以试验研究中实际测得的裸锚杆磁致伸缩信号为对象进行降噪结果对比。裸锚杆原始信号波形如图4所示。

图4 裸锚杆原始信号波形

由图可知，信号中存在大量噪声信号，端面反射信号淹没在噪声中，很难判断反射时刻，这使得锚杆的计算长度与实际长度相差很大。

3.1 CS算法优化VMD参数分析与验证

为了验证CS算法与传统频率观察法相比在优化VMD参数方面的优势，选取试验研究中的实际测得的裸锚杆磁致伸缩信号进行以下分析：对实测信号选取不同K值进行VMD运算，惩罚因子α使用默认值，即α=2 000，观察每个模态分量的中心频率，如表1所示。由表可知，当K=8时，出现了模态分量中心频率相近的现象(过分解现象)，因此K的取值为7。

表1 不同当前分层数K值时各模态分量的中心频率

利用CS算法优化VMD参数组合[K,α]。CS算法参数设置如下：更新概率Pa=0.25，鸟巢数n=20，分解层数K的搜索范围为(2,15)，惩罚因子α的搜索范围为(100,3 000)，最大迭代次数为50。寻优结果为K=7，α=1 525，分解结果如图5所示。

计算传统频率观察法与本文中方法VMD运算后每层IMF的PE值，如表2所示。

将传统频率观察法与优化后得到的参数组合应用到本文中提出的VMD-ICA降噪方法进行对比。选取PE值较大的IMF构成虚拟通道，因此传统频率观察法取IMF3—IMF6，本文中方法取IMF2—IMF6，分别重构信号。2种方法的降噪效果对比如图6所示，锚杆长度的计算结果如表3所示。由于ICA没有关于原始信号的任何先验知识，幅值具有不确定性，因此，降噪后的信号幅值可能会有改变，但对信号分析没有影响。由表3可知，2种方法都能有效地降低噪声影响，但是本文中提出的方法能更大限度地恢复原始信号，计算锚杆长度误差较小。

综上所述，针对模态个数K和惩罚因子α很难确定的问题，CS算法能有效地搜寻到一组最佳输入参数。

图5 模态分量波形及频谱

表2 各层本征模态分量(IMF)的排列熵(PE)值

(a)传统频率观察法(b)本文中提出的方法图6 传统频率观察法与本文中提出方法的降噪效果

表3 传统频率观察法与本文中提出方法的锚杆长度计算结果

3.2 锚杆长度实验分析

为了验证本文中提出的方法相较于其他降噪方法的优势，分别采用小波阈值去噪、EMD-ICA联合降噪加小波阈值处理[13]、VMD-ICA降噪加小波阈值处理[14]及本文中提出的方法进行裸锚杆降噪处理，结果如图7所示。原始信号和上述4种降噪方法得到的锚杆长度值和相对误差率如表4所示。

(a)小波阈值去噪法(b)文献[13]中方法(c)文献[14]中方法(d)本文中提出的方法图7 裸锚杆降噪效果对比

表4 锚杆长度计算结果

文献[13]中各IMF与源信号的互相关系数如表5所示，文献[14]中参数组合为K=7，α=2 000，各层IMF的峭度值如表6所示。本文中提出方法的参数在3.1节中已给出，此处不再赘述。

由图7和表4可知，小波阈值去噪法和文献[13]中方法降噪效果不明显，波形不平滑，含噪声较多，锚杆的计算长度与实际长度误差较大；文献[14]中方法波形平滑，但是误差率最大；而本文中提出的降噪算法的降噪效果好，含噪声较少且波形平滑，最大限度地恢复了原始信号，计算锚杆长度误差最小，仅为0.83%。

表5 各层本征模态分量(IMF)与原信号的互相关系数

表6 各层本征模态分量(IMF)的峭度

4 锚杆和锚固实验分析

测试中采用的锚杆长度为2 m，直径为22 mm，锚固长度为1 m，锚固层材质为混凝土，脱锚段为沙土。锚固锚杆信号采集示意图如图8所示。假设锚杆、锚固长度分别为L1、L2，超声导波在锚杆、锚固中的传播波速分别为c1、c2，感应线圈接收到锚固底端面反射信号的时间为t1，裸锚杆波速c1一般取5 000 m/s，经验证，锚固中传播速度c2为3 000 m/s。计算锚固长度的公式为

(18)

图8 锚杆和锚固信号采集示意图

锚杆和锚固去除电磁脉冲后的原始信号的波形如图9所示。由图可知，信号被噪声所污染，缺陷回波和底端回波被噪声淹没。VMD参数优化参数组合为分解层数K=9、惩罚因子α=1 424、每层IMF的PE值，如表7所示。

图9 原始波形

表7 各层本征模态分量(IMF)的排列熵(PE)值

分别应用4种方法得到的降噪后波形，如图10所示。原始信号和4种降噪方法所得锚固长度值和相对误差率如表8所示。

(a)小波阈值去噪法

(b)文献[13]中方法

(c)文献[14]中方法

(d)本文中提出的方法图10 采用不同方法的锚杆和锚固系统降噪效果

表8 锚固长度计算结果

由图10和表8可知，VMD-ICA降噪方法降噪效果不明显，缺陷回波和底端回波不明显，小波阈值和EMD-ICA去噪波形不平滑，缺陷回波信号和底端面回波有失真现象，而本文中提出的降噪算法的降噪效果好，含噪声较少，能明显看出缺陷回波及底端回波信号，计算锚固长度相对误差最小，仅为0.5%。

5 结论

本文中采取CS算法、VMD与ICA相结合的方式，利用CS算法对VMD的最优参数组合进行搜寻，借助VMD方法较好的噪声稳健性和自适应分解的特性，以及ICA通过信号的不同特征，能够分离出未知原始信号能力，更大限度地去除噪声。通过实测信号结果分析来验证本文中所提方法的有效性。本文中提出的方法应用在锚杆磁致伸缩导波检测上，能更有效地保留有用的反射波信号。结果表明，相对于传统小波阈值去噪法、VMD-ICA降噪以及EMD-ICA降噪方法，本文中提出方法的降噪效果有较大改善，计算锚杆长度的误差更小；在锚固的质量检测中能更有效地提取缺陷信号及回波信号。