基于EEMD和低秩稀疏分解的超声缺陷回波检测方法

周航锐, 孙 坚, 徐红伟, 缪存坚, 宋 鑫

(1.中国计量大学 机电工程学院,浙江 杭州 310018；2.浙江省特种设备检验研究院,浙江 杭州 310012)

1 引 言

超声技术在无损检测中得到了广泛应用,但在实际工程使用超声对金属材料中的微小缺陷检测时受随机噪声和金属材料晶粒散射引起的相干噪声影响,缺陷回波检测仍具有挑战性。国内外学者提出了许多从含噪超声信号中识别缺陷回波的方法,如小波变换[1]、裂谱分析[2]、稀疏表示[3]等。

经验模态分解(empirical mode decomposition, EMD)可以自适应地将信号分解为有限个本征模态函数(intrinsic mode functions, IMF),广泛应用于非平稳非线性信号的处理[4～6]。在EMD的基础上,通过多次添加高斯白噪声的集合经验模态分解(ensemble empirical mode decomposition,EEMD)被提出,克服了EMD分解信号发生模态混叠的问题。

选择信号主导的IMF分量进行局部重构是一种经典的EMD去噪方法,文献[1]应用EMD滤除聚焦超声谐波的噪声分量；文献[2]将EEMD应用于SO2荧光信号去噪。然而,由于EMD算法分解信号的不间歇特性导致大量噪声会被分解到信号主导的IMF分量中,所以需要对重构后的信号作进一步处理。进一步降噪的方法可以大体分为两类,第一类是采用诸如小波阈值去噪的方法对信号主导的IMF分量进行阈值处理后再重构信号。文献[3]研究了基于改进的EMD和局部熵的超声信号阈值降噪方法,通过相关系数法选取信号主导的IMF分量进行局部熵阈值降噪处理进一步抑制结构噪声；文献[4]应用EMD间隔阈值对信号主导的IMF分量阈值处理后加窗提取缺陷回波信号再重构。第二类是先重构信号主导的IMF分量后再结合其他方法实现二次降噪。文献[5]提出对超声信号EMD局部重构初步去噪后结合S变换二次去噪的缺陷回波检测方法。

近年来,低秩稀疏分解(low-rank sparse decomposition, LRSD)作为一种新的矩阵分解理论被提出并广泛用于视频前背景分离,图像去噪、语音增强和故障诊断等领域。LRSD算法在充分考虑数据矩阵特征的基础上,可以鲁棒的将矩阵分解为低秩部分和稀疏部分,又称为鲁棒主成分分析(robust principal component analysis, RPCA)[6]。

本文提出了一种基于EEMD和低秩稀疏分解的超声缺陷回波检测方法,克服以往去噪方法未充分考虑超声缺陷回波在时域和时频域具有稀疏特性的缺点。

该方法首先对含噪信号进行EEMD得到一系列IMF分量,计算原始信号和每个IMF分量的概率密度函数(probability density function, pdf)并根据相似度测量准测选取信号主导的IMF；然后基于短时傅里叶变换(short time Fourier transform, STFT)计算幅度谱和相位谱,执行低秩稀疏分解算法将幅度谱分解为稀疏、低秩和噪声3部分并舍弃噪声部分；最后通过时频掩蔽分离出缺陷信号幅度谱并结合相位谱逆STFT变换得到回波信号。实验结果表明该方法能够在保留缺陷回波完整性的同时有效减少无用信息并提高缺陷识别率。

2 基于EEMD和低秩稀疏分解的去噪方法

2.1 基于EEMD的高频成分去除

原始信号x(t)经EEMD后可以得到一系列IMF分量和余项r(t)：

(1)

式中：h(i)(t)(i≤N)为第i阶IMF分量。在本研究中,使用基于pdf和l2-norm的相似性测量方法来区分噪声主导分量和信号主导分量[7]。pdf是表示随机变量取值分布规律的数学表示,不同的分布可以反映两个变量之间的差异。用L(i)表示x(t)与第i个IMF分量之间的相似性,其定义如下：

(2)

式中：dist是衡量两个pdf之间差异的度量模型。本文采用l2-norm,定义为：

(3)

式中：P和Q分别表示2个信号的pdf。

IMF分量的频率分布随着阶数增加逐渐向低频集中,根据相似性测量准则,噪声主导的IMF分量中信号含量较少,且随阶数递增高频的成份减少导致相似度是逐渐降低的,即L(i)呈现递增趋势,当IMF分量开始由信号主导时相似度较上一阶分量会变高,所以通常选择L(i)达到第一个极大值后点作为分界点,记为kth=argmax1≤i≤L{L(i)}+1,其以后的分量作为信号主导分量。考虑到超声回波是瞬时信号,理想情况下只能在有限的时间间隔内出现振荡,而高阶IMF分量即使不包含高频噪声,仍会出现类似调制正弦波的低频虚假振荡,即与原始信号的相似度越来越低,导致即使分解到首阶信号主导的IMF分量也不会有明显的相似度提高,但经大量实验发现,在首阶信号IMF出现时,相似度L(i)的变化量L=L(i)-L(i-1)会出现极小值,故重新定义kth：

kth=first arg min2≤i≤L{L}

(4)

(5)

2.2 基于低秩稀疏分解的信噪混叠抑制

2.2.1 低秩稀疏分解算法

传统的低秩稀疏分解模型是非凸的NP难问题,一系列求解该模型的算法被众多学者提出。半软去分解(semi-soft godec, SSGD)是在RPCA的基础上优化了矩阵的分解模型而提出的一种低秩稀疏分解算法[8]。

假设受到噪声干扰的矩阵为Y,在其满足低秩和稀疏的优化准则下,SSGD算法可以有效估计出Y中的低秩部分L和稀疏部分S,如式(6)所示：

Y=L+S+N, rank(L)≤r, card(S)≤k

(6)

式中：N代表噪声部分；rank(L)表示矩阵L的秩；card(S)表示矩阵S的基数,即非零元素的个数。对于式(6)可以通过最小化分解误差求解：

s.t. rank(L)≤r

(7)

(8)

2.2.2 低秩矩阵构建

(9)

式中：N为信号长度；W为帧长；h(n)(n=0,…,W-1)为归一化窗函数；R=ti-ti-1为相邻两帧交叠的点数,计算得到的结果即为描述原始信号的时频谱。

2.2.3 基于LRSD的二次去噪

低秩稀疏分解算法处理的对象是含有潜在低秩和稀疏结构特征的矩阵,超声检测中缺陷反射信号为不平稳的信号,其在时频域的稀疏性是一个广泛接受的观点,而缺陷检测中的非声学噪声信号多以自相似的背景混响为主,在时频域具有相似、重复的结构特点。时频谱幅值矩阵X可以视作一帧像素灰度值可为任意非负值的灰度矩阵,从而可以使用低秩稀疏分解算法估计幅度谱X中的低秩和稀疏成份以进一步抑制信噪混叠。

(10)

式中：l和n代表幅度谱中时频点索引。将计算出来二值时频掩蔽值作用于估计出的低秩和稀疏部分的幅度谱,得到缺陷检测信号的幅度谱：

(11)

2.3 基于EEMD-LRSD的回波检测方法

基于EEMD和低秩矩阵分解的超声缺陷回波检测步骤如下。

步骤1：对原始信号x(t)进行EEMD得到N个IMF分量。

步骤4：使用SSGD算法分解幅度谱X得到低秩部分L、稀疏部分S和噪声部分N,舍弃噪声部分。

算法流程如图1所示。

图1 缺陷回波信号检测方法总体设计框图

3 实验与分析

3.1 缺陷检测信号建模

金属缺陷超声检测采集的缺陷回波信号通常还包含两类干扰信号：由测量系统、环境波动等因素引起的非相干噪声,通常可以建模为高斯白噪声；由材料晶粒散射形成的相干噪声,也称结构噪声。在脉冲反射式超声检测中,单个缺陷回波s(t)可以表示为[11]：

s(t)=βe-α(1-rtanh(m(t-τ)))(t-τ)2cos(ωc(t-τ)+φ)

(12)

式中：β为幅值；α为带宽因子；m为双曲正切函数阶次；r为不对称因子；τ为到达时间；a(t)为单位峰值的包络函数；ωc为回波信号的中心频率；φ为初相位。

入射超声波经过金属材料晶粒散射形成的结构噪声可以建模为[10]：

(13)

式中：b为大于零的常数；γ为材料衰减因子；K为超声探测声束范围内的晶粒数；σk为第k个晶粒散射回波的强度系数；ωk为第k个晶粒散射回波的频率漂移；τk为第k个晶粒散射回波延时；σk服从瑞利分布,ωk和τk服从均匀分布。

非声学噪声n(t)通常可以简化为方差σ2的高斯白噪声。综上,缺陷回波检测信号x(t)可以表示为：

(14)

式中：n为缺陷回波数量。为简化表述缺陷回波信号si(t)的各个参数,记Λi=(βi,αi,ri,mi,τi,ωci,φi)为si(t)的参数向量。

利用以上模型生成超声无损检测信号,如图2所示,采样频率取fs=100 MHz,采样时间4 s,参数向量设置为：缺陷回波数量n=3,Λ1=(1,35,0.5,10,0.8,5,0),Λ2=(0.9,30,0.2,12,2,5,0),Λ3=(0.8,35,0,10,3.2,5,0),构造s(t),设晶粒数量K为2500,衰减系数γ为10-28构造u(t),再添加能量大小为-3 dBW的随机噪声n(t)得到信噪比为-10.5 dB的含噪信号x(t)。

图2 仿真含噪信号

3.2 仿真去噪与结果分析

对仿真构造的原始信号x(t)进行EEMD得到如图3所示的8阶IMF分量和一个余项R。

图3 原始信号的EEMD分解结果

3.2.1 信号IMF分量的筛选

图4绘制了每阶IMF分量和原始信号x(t)的pdf分布,并给出了基于l2-norm计算的pdf相似度L(i)。可以看出首阶IMF分量和原始信号相似程度最大,随着IMF的阶数增加,IMF分量和x(t)的L值逐渐增大,即相似程度逐渐降低,但L的变化量L是逐渐减少的,说明IMF分量中依旧是高频成份占主导,但是含量逐渐减少。当阶数增加到4时,L达到最小,故可以确定此后的IMF分量逐渐由信号主导,并且由EEMD的特性也可得知低阶IMF分量不含高频成份。

图4 原始信号及其每阶IMF分量的pdf分布

为了进一步验证信号IMF分量筛选方法的稳定性,改变信号x(t)中的噪声强度,在不同信噪比下分别生成80组信号进行EEMD分解并统计kth的分布,如图5所示,可见筛选方法较为稳定。根据图6(a)中L变化曲线可以确定噪声IMF与信号IMF的分界点kth为4,根据式(5)对IMF4～IMF8分量以及余项R进行重构,得到的结果如图6(b)所示,可以看出EEMD重构的信号相较于原始信号缺陷回波位置显著,高频噪声得到有效消除,但在时域依旧有明显的信噪混叠现象。

图5 不同信噪比下kth统计直方图

图6 相似性测量值及EEMD重构的信号

3.2.2 低秩稀疏分解

图7 低秩稀疏分解二次去噪

3.2.3 去噪效果评价

为了定量描述上述方法的性能,引入相关系数(COR)和均方根误差(RMSE)作为评价标准。COR的计算公式为：

(15)

RMSE的计算公式为:

3.2.4 不同方法去噪结果对比

为了评估本文所提方法对缺陷检测信号处理性能,引入小波阈值(wavelet thresholding denoise, WTD),经验模态分解直接阈值(EMD direct threshold, EMD-DT),EMD离散小波变换(EMD discrete wavelet transform, EMD-DWT)与本文方法进行比较,处理结果如图8所示。

图8 各方法去噪结果对比

由时域波形图可以看出,4种方法均可以从含噪信号中恢复出缺陷信号的波形,但是前3种方法处理后的信号均产生了不同程度的畸变,波形中仍含有较多噪声,起振位置难以识别。本文所提方法处理后的信号整体效果最好,波形畸变程度最小,在缺陷回波之外几乎不含有干扰杂波,信号失真更小,3个缺陷回波的起振位置和振幅以及其他细节信息更易于识别。

为了进一步比较4种方法的性能,计算了在不同信噪比下所得去噪信号的COR和RMSE并绘制曲线如图9所示。

图9 不同SNR下降噪效果比较

可以看出,随着含噪信号的SNR增加,4种方法的降噪性能均有不同程度的提升,但相较于其他3种方法,本文方法即使在信噪比较低时也可以保持较稳定的降噪性能。

3.3 实际检测回波去噪

3.3.1 缺陷回波检测系统

应用本文方法处理实测信号,采用如图10(a)所示的超声纵波脉冲反射法检测缺陷,收发一体式超声探头发射脉冲信号后会多次接收到表面波、缺陷回波、底面波和干扰信号。如图10(b)所示,检测对象是侧面依次钻有2个直径为2 mm孔的铝合金,实验使用RPR-4000型脉冲发生器激励中心频率为5 MHz的多浦乐N5P6L型防水直探头产生超声信号,对其中一个缺陷进行检测,回波信号由RIGOL公司的DS2102型示波器采集并存储。

图10 缺陷检测实验系统

实验采集的超声信号如图11(a)所示,回波信号中受到了较多的噪声干扰。为了简化分析,本文仅截取超声检测信号表面波和底面波之间信号部分数据段进行处理,如图11(b)所示。

图11 金属缺陷超声检测信号

3.3.2 信号去噪及评价

分别使用和上节相同的4种方法对信号进行去噪处理,结果如图12所示。

图12 去噪结果对比

可以看出,4种方法处理后的信号都能较大程度抑制干扰信号,但前3种方法处理后的信号在缺陷回波外信号起伏较大,不利于波形起振位置的识别,而本文方法得到的信号无论从波形完整性还是波形平滑性上来看均优于其他算法,去噪后信号仅含有少量的干扰信号。其中EMD-WTD法处理后的缺陷回波发生了一定程度的畸变,波峰突变比较严重,波形的平滑度也有所下降；EMD-DT法在抑制干扰信号的同时也削弱了缺陷回波幅值,但可以清楚识别回波的起振位置,有利于更精准地定位缺陷；经WTD法和本文方法处理后的回波信号畸变程度较小,回波幅值和波形完整性也得到较好保留,相较于WTD法,本文方法处理后的信号能更好抑制缺陷回波以外的干扰信号,更有利于缺陷回波的精确检测。

5 结 论

(1)针对传统基于EMD的方法在对缺陷检测信号去噪时未充分考虑超声检测信号在时域以及时频域有稀疏特性的缺点,引入低秩稀疏分解,提出了EEMD和低秩稀疏分解相结合的去噪方法。

(2)分析了含噪信号中噪声和缺陷信号的特性,并以STFT为基础构建了EEMD重构信号的低秩稀疏特性矩阵,使用低秩稀疏分解算法成功的提取出了矩阵中的稀疏成份,改善了EEMD重构信号被杂波干扰的状况。

(3)利用仿真信号和铝块缺陷检测信号对本文所提方法和WTD,EMD-DT以及EMD-DWT方法进行了试验,并应用均方差和相关系数量化评价了本文所提方法的信号去噪性能。充分证明了本文方法可以显著改善信号质量,明显抑制信噪混叠现象。