基于自适应低秩稀疏分解的超声缺陷回波检测方法*

2022-01-17周航锐孙徐红伟缪存坚

传感技术学报 2021年11期

周航锐孙坚*徐红伟缪存坚

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

超声技术在无损检测中得到了广泛的应用，但由于电子噪声和相干噪声的存在(材料的晶粒散射回波和噪声影响存在)，在实际工程使用超声对金属材料中的微小缺陷检测时仍具有挑战性。国内外学者提出了许多从含噪超声信号中检测缺陷回波的方法，如小波变换[1-2]、经验模态分解[3]、稀疏表示[4]等。

近年来，低秩稀疏分解(Low-Rank Sparse Decomposition，LRSD)作为一种新的矩阵分解理论被提出并广泛用于视频前背景分离[5]，图像去噪[6]、语音增强[7]和故障诊断[8]等领域。

LRSD算法在充分考虑数据矩阵特征的基础上，可以鲁棒地将矩阵分解为低秩部分和稀疏部分，又称为鲁棒主成分分析(Robust Principal Component Analysis，RPCA)[9]。传统的RPCA模型是非凸的，且是一个NP(non-deterministic polynomial)难问题，一般不易直接求解。在RPCA的基础上，一系列改进的算法被众多学者提出，Zhou等人考虑到在干扰较多的环境中RPCA模型分解得到的结果不理想，提出去分解(Go Decomposition，GoDec)算法认为观测矩阵是由前景、背景和噪声组成，同时使用了双边随机投影(bilateral random projections，BRP)算法加速矩阵前背景分离[10]。该算法虽然在一定程度上加快了模型求解的速度，但是存在难于选择合适参数对秩、稀疏度范围约束的缺陷。由此Zhou等人优化了GoDec算法，提出半软去分解(Semi-Soft GoDec，SSGD)算法[10]，稀疏度约束值可以通过软阈值处理的方式自适应确定，但是该算法依旧需要对背景矩阵的秩进行估计。在背景矩阵的秩估计方面，Li Li等人通过计算待分解矩阵的信息熵线性拟合估计背景矩阵的秩，但信息熵与背景矩阵的秩不完全满足线性关系会影响秩估计结果[11]。Jain P等人提出基于奇异值投影(Singular Value Projection，SVP)的秩估计方法在背景矩阵秩较低时会出现迭代不连贯[12]。

本文在李安冬等人研究低秩度参数和重构精度关系的基础上[13]，将基于误差重建的背景矩阵秩估计方法引入到半软去分解算法中，提出了一种自适应低秩稀疏分解的超声缺陷回波检测方法，对金属材料超声缺陷检测信号进行数据处理，可以有效减少无用信息并提高缺陷识别率。该方法首先对含噪信号进行短时傅里叶变换(short time Fourier transform，STFT)并提取幅度谱和相位谱；然后采用误差重建的方法估计幅度谱背景矩阵的秩，基于估计的低秩度参数执行SSGD算法将带噪信号幅度谱分解为稀疏、低秩和噪声3部分并舍弃噪声部分；最后通过时频掩蔽分离出缺陷信号幅度谱并结合相位谱逆STFT变换得到缺陷回波信号。

1 自适应低秩矩阵分解

1.1 SSGD算法

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

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

对于式(2)可以交替求解以下两个参数的子问题:

式中:‖·‖F表示Frobenius范数。针对GoDec算法的参数r选择过大会导致一部分噪声泄露到稀疏部分S中的缺陷，Zhou等人提出了半软去分解算法[10]:

式中:λ是软阈值参数，‖·‖1表示l1范数。不同于GoDec算法对低秩部分L和稀疏部分S都进行硬阈值处理，SSGD对式(1)中S项采用软阈值处理，通过软阈值λ自适应确定约束“card(S)<＝r”中的稀疏度r的值。对于式(4)可以通过式(3)表述的交替优化方法求解。此外，为了避免对大规模矩阵进行奇异值分解，SSGD算法引入双边随机投影(Bilateral Random Projections，BRP)来求解式(3)，从而进一步提高了计算效率。

1.2 低秩矩阵构建

低秩稀疏分解算法处理的对象是含有潜在低秩和稀疏结构特征的矩阵，考虑到超声信号的时频谱图能够突出信号在不同时间下频率能量分布特征，而缺陷检测中的非声学噪声信号多以自相似的背景混响为主，在时频谱图上具有相似或重复的结构特点，超声检测中缺陷反射信号则为不平稳的信号，其在时频域的稀疏性是一个广泛接受的观点[4]。因此，本文采用短时傅里叶变换构建低秩矩阵。

通过有限长的窗口移动加权对非平稳的原始信号进行交叠分段得到短时平稳的数据帧，然后对预处理后的每一帧信号进行快速傅里叶变换(fast Fourier transformation，FFT)，其计算公式为

式中:N为信号长度，W为帧长，h(n)(n＝0，…，W-1)为归一化窗函数，R＝t i-t i-1为相邻两帧交叠的点数，计算得到的̇Y即为描述原始信号的时频谱。

1.3 自动秩估计

SSGD算法在求解时频谱的低秩和稀疏成份前，需要预估低秩成份的秩作为迭代求解的先验信息[13]，准确的秩估计值可以避免分解出的稀疏成份受到噪声成份的干扰。但是受采样率和杂波的影响，背景实际的秩往往难以准确估计。当预估背景矩阵的秩明显小于矩阵实际的秩时，分解得到的低秩成份L会有大量干扰信号泄露到稀疏成份S和噪声成份N中，则缺陷回波信号可能会被虚假回波干扰导致微小信号漏检；当设定的秩值逐渐递增但依旧小于实际的秩之前，低秩成份L中包含的信号量会逐渐增加，记重建误差为‖N‖F＝‖Y-(S+L)‖，其变化量为ΔN，此时重建误差的变化量较大；当背景矩阵的秩设定接近或者大于实际的秩时，分解得到的低秩成份L包含了大量具有重复结构特征的干扰信号，此时重建误差较之前有明显的下降后，同时其变化量也趋于稳定，但是随着背景矩阵的秩设定越来越大时，稀疏部分S中的缺陷回波信号又会混入低秩成份中导致漏检。因此根据信号重建误差的变化趋势自适应确定SSGD算法中的参数k，可以保证幅度谱中的缺陷回波和干扰信号有效区分。逐步递增设定的秩值并分解幅度谱，得到重构信号和重构误差。当重建误差平稳时，有‖ΔN‖F≤η。选择重建误差平稳后对应的第一个秩值作为k值输入到SSGD算法。算法具体流程如图1所示。

图1 低秩度参数估计流程图

1.4 增强信号波形重构

对时频谱取模即可得到幅度谱Y，在使用SSGD算法估计出Y中分别代表稀疏特征和重复特征的幅度谱S和L后，需要采用二值时频掩蔽法[14]处理以得到最终需要的缺陷回波幅度谱。二值时频掩蔽法假设S和L的混合幅度谱中每个时频点只可能包含一种信源，故可以通过保留只包含稀疏成份的时频区域而去除包含低秩成份的区域，最大限度的消除稀疏成份中的干扰信号。本文的二值时频掩蔽计算公式如下:

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

结合原始信号STFT时记录下的相位信息∠·，对估计出幅度谱进行逆STFT就可以得到最终的回波信号。

1.5 基于自动秩估计的SSGD矩阵分解算法

基于自适应低秩矩阵分解的超声缺陷回波检测步骤如下:

步骤1 根据式(5)使用汉明窗(Hamming)对观测信号进行STFT得到时频谱，取模得到幅度谱，记为Y。

步骤2 在各低秩度参数下计算重建误差‖N‖F，记录满足误差变化量小于η的第一个低秩度参数k。

步骤3 根据式(3)和式(4)在秩为k的条件下对幅度谱Y进行分解，得到噪声部分N、低秩部分L和稀疏部分S，记去除噪声部分N后的信号幅度谱为Y L S

步骤4 根据式(6)计算掩蔽值M，再根据式(7)对幅度谱Y LS进行二值掩蔽后可得到缺陷信号的幅度谱。

步骤5 最后对做逆STFT得到降噪后的时域回波信号。

算法流程图如图2所示。

图2 方法流程图

2 实验

2.1 仿真缺陷回波检测

超声波是一种非平稳时变信号，缺陷检测接收的信号包括缺陷回波、材料晶粒散射形成的结构噪声和非声学噪声等。在脉冲反射法中，单个缺陷回波s(t)可以表示为。

式中:β为幅值；τ为到达时间；φ为初相位；a(t)为单位峰值的包络函数；f c为中心频率。

考虑实际的缺陷回波包络可能不对称(通常为上升沿陡峭而下降沿平缓)，引入调幅回波模型(Amplitude-modulated echo model，AMEM)[15]，取

式中:α为带宽因子；m为双曲正切函数阶次；r为不对称因子-1

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

式中:b为大于零的常数；γ为材料衰减因子；K为为超声探测声束范围内的晶粒数；σk为第k个晶粒散射回波的强度系数；ωk为第k个晶粒散射回波的频率漂移；τk为第k个晶粒散射回波延时；ωc为超声缺陷回波信号的中心频率。

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

为简化表述缺陷回波信号s i(t)的各个参数，记λi＝(βi，αi，r i，mi，τi，fci，φi)为s i(t)的参数向量。

利用以上模型生成超声无损检测信号，参数设置如下:缺陷回波个数N＝2，λ1＝(1，35，0.5，10，0.8，5，0)，λ2＝(0.8，35，0，10，3.2，5，0)构造回波信号s(t)。在信号长度范围内材料的晶粒数量K设为2 500，衰减系数γ为10-28，σk服从瑞利分布[4]，ωk和τk服从均匀分布，构造u(t)。再添加能量大小为-3 dBW的随机噪声n(t)得到信噪比为-10.5 dB的含噪信号x(t)。仿真信号的采样频率取作f s＝500 MHz，采样长度为2 000，仿真结果如图3所示。

图3 仿真信号

从图3(b)可见，仿真信号x(t)信噪比很低，缺陷回波被淹没。图4为提取幅度谱后经自适应秩估计得到的不同低秩度参数下低秩稀疏分解算法重建矩阵产生的误差变化量曲线，其中误差变化量经归一化处理，η根据大量仿真实验在本文中取＝0.1。图5展示了部分低秩度参数下自适应矩阵分解得到的幅度谱低秩部分L、稀疏部分S和噪声部分N经逆STFT的时域波形。从图中可以看出，低秩度参数设定低于实际值时，有大量干扰信号泄露到稀疏成份中；低秩度参数越高，分解得到的低秩矩阵包含的干扰信号越多，但稀疏矩阵中的缺陷信号幅值逐渐减弱甚至出现了丢失第二个幅值较小的缺陷回波的情况。所以准确估计待分解矩阵背景的秩，可以保证检测率的同时可以有效降低泄露到缺陷回波的噪声。对于本文的仿真信号，当k设定由2增加为3时，相应的误差变化量首次低于η值，所以选择低秩度参数为3进行低秩稀疏分解。对得到的低秩和稀疏部分进行二值时频掩蔽处理，最终提取的缺陷回波信号如图6所示，对比原始回波信号可以看到噪声基本被消除，回波显著，畸变程度较小。降噪后的信号幅值虽有所降低，但对于强噪声下的缺陷检测和定位并无影响。

图4 重建误差变化曲线

图5 不同秩参数约束条件下分解得到的L、S和N对应的时域波形

图6 去噪信号(t)和回波信号s(t)对比图

为了评估本文所提方法在不同强度的噪声下对缺陷检测信号处理性能，改变仿真信号中噪声的强度，分别得到信噪比为-12.17 dB、-8.56 dB、-5.55 dB、-3.09 dB四种不同的含噪信号x(t)，同时引入小波阈值(Wavelet thresholding denoise，WTD)，经验模态分解(Empirical Mode Decomposition，EMD)，EMD离散小波变换(EMD discrete wavelet transform，EMDDWT)与本文方法进行比较，处理结果如图7所示。由时域波形图可以看出，随着信噪比的提升，四种方法的性能均逐步提升，但是前三种方法处理后的信号均产生了较大程度的畸变，尤其是在信噪比较低时，波形中仍含有较多噪声，起振位置难以识别。其中EMD-DWT在信噪比较高时可以达到和本文所提方法相同的性能，但是其波形前后仍含有少量噪声。本文所提方法处理后的信号整体效果最好，波形畸变程度最小，在缺陷回波之外几乎不含有干扰杂波，信号失真更小，两个缺陷回波的起振位置和振幅以及其他细节信息更易于识别，当信噪比处在较低水平时也依旧可以保持稳定的性能。

图7 不同噪声强度下四种方法的处理结果对比

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

式中:和xorg分别是处理后的信号和原始信号，〈·〉表示内积运算。COR用来评价重构缺陷信号与原始回波信号之间的相似程度。

RMSE的计算公式为:

式中:N为用于分析的信号采样长度。RMSE用于衡量重构缺陷信号(t)与原始回波信号s(t)之间的误差。其中COR越大，RSME越小，去噪效果越好。去噪后信号的各项指标如表1所示。

分析表1可以得到和时域波形对比类似的结论，即相对于原含噪信号x(t)，在经低秩矩阵分解处理后，信噪比和互相关系数的提升量最大，均方差的下降量最大，说明低秩矩阵分解可以在很大程度上去除噪声，同时还能保持回波信号失真较小。

表1 不同噪声强度四种方法处理结果评价指标

2.2 实际缺陷回波检测

应用本文所提方法处理实测超声缺陷回波信号。采用基于脉冲反射法的水浸超声纵波对缺陷进行检测，检测示意图如图8(a)所示，收发一体式超声探头产生具有一定频率间隔和持续时间的脉冲进入工件，在声阻抗变化处发生不同程度的反射，其中按照入射路径返回的部分脉冲被探头接收。实验中检测对象为如图8(b)所示的6061铝合金，在侧面依次钻有两个直径2 mm的孔作为缺陷，使用其中距离上平面d＝60mm的孔F1作为代表进行检测，材料纵波传播声速为c＝6 300 m/s。实验装置如图9所示，超声换能器中心频率5 MHz，示波器采样频率f s为36 MHz。

图8 缺陷检测实验系统

图9 缺陷检测实验系统

实验采集的超声信号如图10(a)所示，可以观察到回波信号受到了较多高频非声学噪声干扰，分析其频谱(图10(b))发现有大量接近或低于探头中心频率的频率成分，主要由结构噪声和缺陷回波引起[17]。

图10 实验采集的超声信号及其频谱

图10 (a)中信号幅值首次和末次近似达到1 V对应的采样点分别是86和120，取均值作为始波位置，记为n＝103。通过传播声速c和缺陷实际位置d计算得到缺陷回波的大概时间值t f＝n/f s+2d/c＝2.86μs+19.05μs＝21.91μs，底波可以通过同样的方式计算得到，为了简化分析，本文仅截取始波和底波间800个采样点的数据进行分析，如图11所示。

图11 用于分析的信号段

最终处理结果如图12所示，缺陷回波在图中的采样点对应的时间值在22.08μs附近，与实际回波时间点偏差很小。同时可以观察到，四种方法都能较大程度抑制干扰信号，但前三种方法处理后的信号在缺陷回波外信号起伏较大，不利于波形起振位置的识别，而本文方法得到的信号无论从波形完整性还是波形平滑性上来看均优于其他算法，去噪后信号仅含有少量的干扰信号，回波幅值和波形完整性得到较好的保留，更有利于缺陷回波的精确检测。