基于VMD方法的混凝土缺陷超声成像噪声处理研究

庄 政，任宏伟，田博文，郑元成，路 乾

(济南大学 自动化与电气工程学院，济南 250024)

0 引言

混凝土材料在建筑工程领域依靠其自身的高强度的稳定性、可塑性、经济实用性等优势被广泛的运用。在国家工程建设中，如桥梁，楼房、水坝等建设领域，混凝土结构已成为建筑材料中不可缺少的一部分，和我们的生活息息相关。正因如此，工程建设安全性的相关研究显得尤为重要。由于施工技术的差异和内置钢材的腐蚀，会使混凝土内部产生蜂窝状缺陷或者是裂缝。因此，需要在不破混凝土结构的情况下，对其内部结构的损伤进行评估。采用超声波对混凝土进行无损检测[1]，根据缺陷位置和其他部位的介质区别引起的反射、散射等信号特性，在不破坏被检结构的正常工作基础上，获得内部结构的损伤信号。对采集的缺陷信号进行处理分析，判断出缺陷位置并对缺陷进行评估。对缺陷信号进行判断分析需要专业的检测人员进行，带有一定的主观性，利用超声信号进行成像处理，能够直观的判断出缺陷位置和损伤程度。

TR-MUSIC算法[2]是基于时间反演(Time Reversal)和MUSIC(Multiple Signal Classification)算法结合的一种超声成像技术，它可以在不需要任何先验性信息的情况下，实现缺陷目标的高分辨率成像。通过压电传感器向混凝土结构发射超声波，声波在被测物体内部经反射、折射后，传感器接收到含缺陷信息的时域信号。通过反演算法将信号逆向传播到缺陷所在位置，即将信号进行时间反转并发送出去，声波将会聚焦源位置。接下来采用MUSIC算法将含缺陷的超声信号进行缺陷目标处理，确定目标缺陷在物体内部的位置和形状，通过计算信号的特征向量和谱分解，来提取物体内的信号。对其信号的阵列响应矩阵进行奇异值分解，通过奇异值划分噪声子空间和信号子空间[3]，分别对应为背景噪声和确定物体的信号分量。由于存在噪声和背景信号的影响，噪声子空间和信号子空间在实际计算中并不是完全正交的。因此，噪声子空间中的向量与信号子空间中的格林函数向量内积的结果并不总是为零，计算时，区域的向量内积的结果取倒数，并将这些结果用于成像。此算法成像分辨率高，计算速度快，但由于计算方式等原因，对噪声的敏感度较高。

在工业无损检测领域，噪声是超声成像检测技术应用的主要限制因素之一。其中，除了电磁干扰、震动干扰和温度干扰等因素产生的电噪声外，对于混凝土超声检测来说，混凝土密度较低，能量衰减大，内部包含各种尺寸的沙石骨料和杂质，结构复杂，且声阻抗较低，与空气和探头之间的匹配度较差，会导致声波反射和衍射较强，进一步增强了噪声的干扰。

通过相控阵超声设备采集在混凝土介质缺陷处的回波信号，因混凝土能量衰减严重，采用高频激励信号探测深度不够，而激励信号的频率与超声成像的分辨率为正相关，为保证探测距离，且能够提高成像分辨率，一般采用100～500 kHZ范围内超声信号进行激励，并进行TR-MUSIC算法超声成像，在激励信号不改变保证探测距离的情况下，提高混凝土缺陷成像分辨率。但是，此算法成像极易收到噪声干扰，混凝土结构因材质的原因，噪声干扰又尤为严重。

针对TR-MUSIC算法抑制噪声问题，目前大多采用多频时间反转多信号分类法(Multi-Frequency TR-MUSIC)进行处理，利用给定频率带宽的频率组进行成像[3]，但由于混凝土内部结构复杂，噪声干扰强度大，此种方式效果并不明显。

为解决混凝土缺陷检测进行TR-MUSIC超声成像噪声敏感度高的问题，提出了一种基于VMD的噪声处理方式，此外将成像处理过程进行改进，从信号处理和成像算法方面互补，抑制噪声对混凝土缺陷成像的干扰问题。

在无损检测的信号处理领域，降噪处理的常规方法有小波降噪、经验模态分解EMD降噪、变模态分解VMD降噪[4]等方式。其中小波降噪方式根据通过分解和设定阈值去除噪声部分进行重构实现降噪处理，能够很好的处理数据相关性，得到广泛应用，但如何设置分解层数和设定阈值对降噪有很大影响。经验模态分解EMD降噪方式其核心为希尔伯特-黄变换，此类相关改进方法广泛应用在信号处理当中。2014年，VMD变模态分解由Dragomiretskiy等人提出，此方法结合Hilbert变换、Wiener滤波器和频率处理方式[5]。进行变分问题的求解，得到各模态分量的带宽和中心频率，此方法是进行迭代搜寻变分模型最优解来实现的[6-7]，对待分解信号的频域和各模态分量进行自适应分解。相对于其他两种降噪方式，VMD方法能够有效的避免端点效应和模态混叠等问题，更有利于对缺陷信号的处理。

1 基本原理

1.1 VMD分解原理

变模态分解方式进行降噪处理，是一种基于Wiener滤波器、Hilbert变换等为基础的信号降噪方法，其本质是将待处理信号分解为多个局部频率模态，并通过优化一个带约束的能量函数来获取这些模态，得到各个模态分量的中心频率和带宽，对模态分量进行分类和降噪处理[8-11]。

VMD的约束变分模型如式(1)所示：

(1)

式中，f为待分解信号，k为分解VIMF分量的个数，ωk、uk分别为对应模态分量的中心频率和振幅。

为求得约束变分问题得最优解，引入拉格朗日乘子，得到拉格朗日函数，其表达式为：

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

(2)

式中，α惩罚参数，λ为拉格朗日乘子。

运用交替方向乘子法(ADMM)更新处理(2)式寻找ωk、uk，计算约束变分模型最优解。

(3)

1.2 TR-MUSIC算法

时间反演多信号分类法将传统的时间反转TR和MUSIC结合应用在超声成像领域内，对缺陷目标检测成像提高成像分辨率效果显著，实现超分辨率的混凝土内部缺陷成像。

将采集到的含缺陷的N*N(对应超声传感器序号)个时域信号进行傅里叶变换得到阵列响应矩阵Ki，j(ω)，其表达式为：

(4)

式中，F(ω)为发射传感器和接收传感器传输的传递函数，τl为散射子的散射强度，φ(ω)是频率相关的传感器的相位响应，v(ω)表示噪声项。

通过阵列响应矩阵K和发射矩阵E、接收矩阵R之间的关系，进行时间反转，经2n次迭代后得到新的输入矩阵，以此作为发射信号在缺陷点聚焦实现时间反演过程。

接下来对阵列响应矩阵进行奇异值分解可得到

(5)

因为阵列响应矩阵K为复对称矩阵，奇异值对应各奇异向量相互正交，理想情况下，响应矩阵进行奇异值分解(SVD)后，非零特征值对应的特征向量和散射子相对应，信号子空间和噪声子空间互相正交，可以得到噪声子空间和对应格林函数向量的内积为0，但在实际计算时候，因存在背景噪声等因素影响，导致内积结果不为0，但值一定接近于0，取其结果的倒数，散射子位置的数值会远大于其它区域。其区域内点的声场强度计算方式如(6)式所示：

(6)

式中，r表示区域点的声压幅值；ωc为传感器激励信号的中心频率取值。

1.3 Savitzky-Golay滤波

Savitzky-Golay滤波[6]是一种数字信号处理技术，通过应用滑动窗口并对该窗口内的数据拟合多项式函数来平滑有噪声的数据。多项式的系数用于估计窗口内数据的平滑值。其表达式为：

s.t.Yi=c0+c1i+c2i2+…+cpip

(7)

通过最小二乘法求解，在每个窗口上，我们将其内的数据点看作是一个多项式函数上的采样点，通过最小二乘法拟合一个多项式函数，得到该窗口内的平滑结果。在整个数据序列上，通过将窗口移动并对每个窗口进行拟合，得到数据序列的平滑结果，可以减少噪声的影响，提高信号的精度。

1.4 能量熵和重心熵

能量熵和重心熵都是用于衡量信号在频域内特性的指标。通过能量熵和重心熵的加权和的方式进行VMD分解后的噪声分量和非噪声分量的分类。

重心熵Hc是用于在频域内能量分布的平衡程度，用于描述信号能量在不同频率成分之间的均衡性。能量熵He是一种用于衡量信号在频域内能量集中的指标。He、Hc可以通过以下公式计算：

(8)

其中：p(x)表示随机变量x在不同取值的概率分布。一般来说，根据采集的含缺陷时域信号能量主要集中在小部分频率的特点，能量熵和重心熵加权和越低，其分量信息越接近混凝土含缺陷的无噪声信号。

2 信号处理过程

2.1 信号采集

采用COMSOL仿真软件搭建含缺陷的混凝土仿真模型。配置填充沙石材质骨料[11]和混凝土，构建尺寸400 mm×400 mm混凝土结构体，随机投放一定数量的不同尺寸沙石骨料，搭建随机骨料混凝土结构模型；在结构体内的(186 mm，200 mm)和(216 mm，200 mm)两个点设置半径10 mm的近距离孔洞缺陷；设置各类材料属性，其中混凝土和沙石骨料属性参数如表1所示。

表1 仿真模型材料参数

设置12个压电超声传感器，压电材料选择PZT-5H型号，置于结构体表面进行超声波无损缺陷检测，传感器宽度为0.63 mm，阵元间隔取波长的一半，构建1-D线性相控阵形式，图1为搭建的混凝土结构超声检测模型。将相控阵内传感器进行1～12标号，1号传感器作为激励阵元，所有传感器(包含1号传感器)作为接收阵元进行信号采集；采集完成后，2号传感器作为激励阵元，所有传感器进行接收。依次类推，采集得到12×12个含缺陷信息的时域信号。其中传感器激励信号为：

图1 随机骨料混凝土缺陷检测模型

S(t)=1e-3*gp1(t)*sin(2*pi*f0*t)

(9)

式中，gp1(t)为高斯脉冲函数，f0为中心频率，取值为150 kHz。

如图2所示，取第8号传感器为激励源，6号传感器为接收源的缺陷采集信号，其中位置1为Lamb波，能量迅速衰减，位置2为缺陷反射信号，位置3、4为混凝土边界的反射信号。为更好的模拟实际混凝土采集的缺陷信号，将所有信号加入热噪声、机械噪声和多路径噪声。机械噪声可以模拟超声传感器或采集设备的机械振动和噪声。这种噪声可能是由于机器的震动或传感器接触到杂散物体造成的。通过添加低频噪声来模拟这种噪声；多路径噪声可以模拟信号传输过程中的多路径传播和散射。通过添加多个不同的信号路径和相位差异来模拟这种噪声。

图2 采集信号和加噪声信号图

图3 变模态分量及对应频谱图

2.2 信号处理

首先将加噪后的含缺陷信号进行VMD分解，选择合适的分解层数K=4，惩罚因子α=2 000。图 3为分解后的变模态分量和对应的频谱图；随后计算各模态分量的能量熵和重心熵，混凝土缺陷超声成像主要是根据中心频率所对应的能量进行成像，根据图2缺陷采集信号的特点，无噪声情况下，其信号能量主要集中在某一小部分固定频率范围内。因此可通过能量熵和重心熵加权和的方法评定噪声分量和非噪声分量；对噪声分量使用小波阈值降噪方法[15]处理，以提高信号质量和提取信号的特征。这里选择了db4小波作为小波基函数，将信号分解为5层，使用自适应阈值选取方法(基于信号局部能量的阈值选取方法)对噪声分量进行处理。将处理后的噪声分量和非噪声分量进行Savitzky-Golay平滑处理并重构得到降噪后的信号及频谱图。如图4所示，对比原信号、含噪声信号和降噪后信号，可以看出原始信号和降噪后的信号的频谱图中心频率和对应幅值接近，降噪后的信号在高频区域的能量明显降低，降噪效果良好。

图4 信号及对应频谱图对比图

通过信噪比(SNR)、均方误差(MSE)评价指标进行降噪效果评估如表2所示；其中SNR表示信号和噪声之间的比率，信号和噪声比例越大，信号质量越高；MSE指用来表示估计值和真实值之间差的平方的平均值，MSE越小，估计值和真实值间的误差越小。经对比含噪声信号、VMD降噪方法、改进后的VMD算法各评价指标。相对于VMD降噪处理方式，改进后算法降噪效果提升明显。

表2 降噪效果评估结果

2.3 时间反演成像

将12×12个含缺陷时域信号进行降噪处理后，构建反演模型，配置信号经时间反演后发射位置和结构空间，采用TR-MUSIC算法进行成像。

首先，对信号进行傅里叶变换，通过降噪信号的频谱图进行中心频率的选取，指定降噪后中心频率=550 kHz为工作频率，提取对应的频谱值。将信号构建为12×12的阵列响应矩阵，由于阵列的第i个发射j个接收的信号和第j个发射i个接收的信号相同，可以得到为复对称矩阵。对阵列响应矩阵进行奇异值分解，并将奇异值进行归一化处理，得到的12个奇异值，奇异值对应混凝土内部的缺陷信息，将非零的奇异值对应的奇异值向量作为信号子空间，其它奇异值向量作为噪声空间。在无噪声的条件下，信号子空间对应的奇异值个数一般对应混凝土内缺陷个数[13]。降噪后的信号虽然无法完全将噪声对应的奇异值变为0，但明显能够和信号子空间区分开来，可通过设置固定阈值进行处理，进行空间划分。经过上述处理，采用式(6)的函数进行TR_MUSIC成像，通过噪声子空间和对应格林函数向量的内积接近0的特点取其倒数，根据区域内声场强度不同进行成像。

3 实验结果与分析

通过COMOSL软件建立的仿真模型进行超声采集实验，混凝土模型缺陷为间隙20 mm的孔洞；经上文提到压电传感器激励信号和阵元采集方式对混凝土进行检测，得到144个含缺陷信息的时域超声回波信号；将采集到的时域信号通过TR-MUSIC算法进行成像处理，设置了对比试验，分别是对采集到的144个信号不添加噪声进行成像；添加噪声干扰进行成像；添加噪声并使用改进的变模态分解方法降噪后再进行成像。此外进行TR-MUSIC算法成像过程中设置划分子空间的统一阈值。

经3种处理方式的信号进行时间反演成像过程中，各类信号归一化后所提取的奇异值点图，如图5所示，可明显看出信号子空间和噪声子空间对应奇异值的区分度提高。降噪后的噪声子空间对应的奇异值更接近于0，设置区分空间的阈值，设定阈值15%进行子空间的划分。

图5 各类信号奇异值对应点图

成像结果如图6所示，含噪声信号进行TR-MUSIC成像，无法定位缺陷点位置，受噪声干扰，产生伪像问题。经信号降噪和设定阈值划分子空间的方法抑制了噪声对成像的干扰，去除了大部分伪像，能够精准定位两个缺陷位置。

图6 各类信号TR-MUSIC成像图

4 结束语

针对混凝土进行TR-MUSIC算法超声成像易受噪声干扰的问题，从信号噪声处理和设定划分空间阈值的方式来抑制噪声对成像的干扰。本文基于VMD对噪声的敏感性，提出一种改进VMD的降噪处理算法，通过重心熵和能量熵加权和的方式进行噪声模态分量和非噪声模态分量的划分，对噪声分量进行小波阈值降噪，并对所有分量进行Savitzky-Golay滤波。对于此类采用TR-MUSIC算法进行超声成像的时域信号，此种方能够根据原信号的能量集中在比较少的频率成分上特点，更好的进行降噪处理。从信噪比、均方误差和结构相似度等评估指标进行对比，改进的VMD降噪算法效果更优。经降噪后的信号进行TR-MUSIC超声成像，在通过归一化的奇异值划分噪声子空间和信号子空间的过程中，改善了因噪声引起奇异值区分度不高的问题。提高了噪声子空间和信号子空间对应奇异值大小的区分度，并通过设定阈值进行空间划分。经降噪后的图像和原图像对比，此算法改善了混凝土缺陷超声成像的噪声干扰问题，能够更精准的识别缺陷目标，抑制了伪像的产生。