CN106556647A - 一种冲击回波数据处理方法 - Google Patents

Abstract

本发明公开了一种冲击回波数据处理方法。本发明的技术要点是：利用冲击回波仪在待测混凝土外侧激发出应力波，并接收冲击回波数据；通过计算相对均方根误差极大值来确定所添加噪声的上限截止频率f_p；根据所确定的噪声上限截止频率，对回波数据进行互补集总经验模态分解(CEEMD)，得到一组本征模态函数IMF；将所有IMF进行希尔伯特变换，求得瞬时频率，得出希尔伯特谱；根据上一步骤的结果在希尔伯特谱中提取出与混凝土内部质量相关的特征。采用本发明的方法使集总次数减少，在需要处理大量回波数据时能大大缩短计算时间；求得的频率峰值更精确，更能反映缺陷的位置信息；图形更直观、清晰。

Description

一种冲击回波数据处理方法

技术领域

本发明属于声波信号处理技术领域，具体涉及一种冲击回波数据处理方法。

背景技术

冲击回波法是一种无损评价混凝土的方法，在诸多无损检测方法中，冲击回波法曾被列为混凝土检测“最有发展前途的现场检测方法之一”，该方法具有只需单面检测，测试过程简便、快捷等优点，可用于精确地、非破损的、经ASTM认可的测量混凝土板块和路面的厚度；也可用于探测素混凝土、钢筋混凝土和后张钢筋混凝土中缺陷的位置和大小，且不受钢筋存在的影响。其原理是：利用冲击产生应力波，该应力波穿入混凝土结构内部，在其内部的缺陷、构件的界面或底端界面处发生多重的反射，来回的反射引起结构的瞬态共振，表现出不同的模态形式，当反射波到达冲击表面，引起表面位移，通过邻近冲击位置的换能器记录下位移-时间信号，传感器的信号输出表现为与位移成比例的电压-时间信号，通常，在计算机里利用数据处理方法—傅里叶变换将电压-时间信号转换为频域波形，不同的模态在频谱上表现为不同的峰值频率，这些峰值频率即与应力波在结构中的多重反射有关。

冲击回波法的基本公式为D＝βV_P/2f，其中，f为某处回波峰值频率，V_p为被测介质的纵波波速，β为截面形状系数，D为回波对应的深度(缺陷到冲击表面的距离或密实结构的厚度)。通常，通过傅里叶变换将回波信号转化为频率信号，进而测得峰值频率来求得缺陷的位置信息。

然而，用以上方法仍然存在如下问题：其一，冲击回波法产生的回波信号为非稳态信号，而傅里叶变换只能针对平稳信号才有意义，并不适用于回波信号的频谱分析。其二，在实际检测中，人们往往不仅需要知道缺陷的位置信息，更需要知晓缺陷的相对尺寸信息，这也是传统傅里叶方法无法直观显示的。因而需要有效的时频分析方法。

较典型的时频分析方法有短时Fourier变换、Wigner分布和小波变换，它们对非平稳信号的处理做出了较大的贡献，在工程实际中也获得了较广泛的应用，但它们是以傅立叶变换为其最终的理论依据。傅立叶变换理论中表征信号交变的基本量是与时间无关的频率，基本时域信号是平稳的简谐波信号。这些概念是全局性的，因而用它们分析非平稳信号容易产生虚假信号和假频等矛盾现象。

对非平稳信号比较直观的分析方法是使用具有局域性的基本量和基本函数。瞬时频率是容易想到的具有局域性的基本量，也是很早就已提出的概念。瞬时频率比较直观的定义是解析信号相位的导数，但以往这一定义会产生一些虚假的结果，导致基于瞬时频率的时频分析方法和理论始终未真正建立和发展起来。而基于EMD(Empirical ModeDecomposition)的时频分析是一种新型分析方法，结合EMD分解和Hilbert谱技术可以将复杂信号分解为有限个本征模态函数(IMF)，从而赋予瞬时频率合理的物理意义，得到信号准确的、定量的时频表示。由于EMD方法存在模态混叠现象，提出了一种噪声辅助数据分析方法——集总经验模态分解(EEMD)，但EEMD方法存在噪声残留和计算量大的问题，互补集总经验模态分解(CEEMD)解决了这些问题。

发明内容

本发明的目的在于提供一种基于互补集总经验模态分解(CEEMD)方法，能更加减少噪声残留、减少计算量的冲击回波数据处理方法。

本发明的冲击回波数据处理方法，包括如下步骤：

(1)利用冲击回波仪在待测混凝土外侧激发出应力波，并接收冲击回波数据；

(2)通过计算相对均方根误差极大值来确定所添加噪声的上限截止频率f_p；

(3)根据步骤(2)所确定的噪声上限截止频率，对回波数据进行互补集总经验模态分解(CEEMD)，得到一组本征模态函数IMF；

(4)将所有IMF进行希尔伯特变换，求得瞬时频率，得出希尔伯特谱；

(5)根据步骤(4)的结果在希尔伯特谱中提取出与混凝土内部质量相关的特征。

具体的，步骤(2)所述噪声的上限截止频率f_p的确定包括如下步骤：

(1)对于回波时域信号x(t)，设定一高斯白噪声ω(t)，其上限截止频率为f_p′＝m·f_s，m取10至20，f_s为x(t)的采样率，以频率f_p′为重采样率，将x(t)进行三次样条插值计算，得到新信号x′(t)；

(2)对信号x′(t)做加噪声的集总经验模态分解，集总次数为2次，所添加的噪声的上限频率为f_p'，噪声标准差为原信号标准差的0.01倍，在不同的噪声上限频率下计算RRMSE的值，公式如下：

其中，c_max(t)为与原始信号相关性最高的IMF成分；

当RRMSE达到最大值时所对应的噪声频率即设为本次CEEMD的最佳噪声上限截止频率f_p。

具体的，步骤(3)所述对回波数据进行互补集总经验模态分解(CEEMD)包括如下步骤：

(1)对回波时域信号x(t)分别施加一对相位相反、幅值相同的高斯白噪声ω(t)，设定噪声标准差为原信号标准差的0.01倍，噪声的上限截止频率为f_p，得到：

x⁺(t)＝x(t)+ω(t)

x^-(t)＝x(t)-ω(t)；

(2)分别对x+(t)和x-(t)进行经验模态分解运算，得到第i阶IMF分量c_i(t)，直到得到两组模态分解的结果，记为：c_i ⁺和c_i ^-，求其均值c_i＝(c_i ⁺+c_i ^-)/2。

具体的，步骤(4)所述得出希尔伯特谱包括如下步骤：

(1)对每个IMF分量c_i(t)作希尔伯特变换，可得数据序列c_i％(t)：

(2)将c_i(t)作为实部、作为虚部构成一个复序列解析信号z_i(t)：

其中，

(3)从而得到瞬时频率：

若忽略残余分量，则原信号x(t)即为z_i(t)的实部，表示为：

(4)以时间、瞬时频率为自变量，幅值在时频平面内以等高线表现，则可将三者的关系表示成一个三维图形，定义为HHT谱，记做H(ω,t)：

本发明提供的冲击回波数据处理方法，与现有技术相比，具有以下优点：

(1)采用本发明的方法使集总次数减少，在需要处理大量回波数据时能大大缩短计算时间。

(2)使用本发明方法求得的频率峰值更精确，更能反映缺陷的位置信息。

(3)能对缺陷的尺寸进行定量分析。

(4)图形更直观、清晰。

附图说明

图1为本发明方法的算法流程框图。

图2为本发明实施例中混凝土梁预应力注浆管道模型的纵向截面示意图。

图3为本发明实施例中混凝土梁预应力注浆管道模型的横向截面示意图。

图4为本发明实施例中注浆密实度为0％时的冲击回波HHT谱。

图5为本发明实施例中注浆密实度为30％时的冲击回波HHT谱。

图6为本发明实施例中注浆密实度为50％时的冲击回波HHT谱。

图7为本发明实施例中注浆密实度为70％时的冲击回波HHT谱。

图8为本发明实施例中注浆密实度为100％时的冲击回波HHT谱。

具体实施方式

下面结合附图和实施例对本发明作进一步详细的描述。

以检测混凝土梁预应力管道注浆密实度为例，利用数据模拟的方法来说明本发明的实用性。采用ANSYS有限元分析软件建立混凝土梁预应力注浆管道模型，如图2、图3所示，尺寸为2m*0.3m，图2中，1为混凝土，2为水泥，3为塑料波纹管；图3中，D为波纹管3的直径，h为波纹管3的注浆厚度；波纹管3位于正中间部分内部的注浆体预设五种不同比例的脱空缺陷，其注浆密实度分别是0％、30％、50％、70％、100％，利用冲击回波仪在五种脱空缺陷外侧激发主频为10KHz的半正弦应力波，并在激励点附近接收回波，利用本发明的方法对数据进行处理。

参见图1，本发明方法的具体步骤如下：

(1)利用冲击回波仪在待测混凝土1的外侧激发出应力波，并接收冲击回波数据；

(2)通过计算相对均方根误差极大值来确定所添加噪声的上限截止频率f_p；

具体过程如下：

(a)对于回波时域信号x(t)，设定一高斯白噪声ω(t)，其上限截止频率为f_p′＝m·f_s，m取10至20，f_s为x(t)的采样率，以频率f_p′为重采样率，将x(t)进行三次样条插值计算，得到新信号x′(t)；

(b)对信号x′(t)做加噪声的集总经验模态分解，集总次数为2次，所添加的噪声的上限频率为f_p'，噪声标准差为原信号标准差的0.01倍，在不同的噪声上限频率下计算RRMSE的值，公式如下：

其中，c_max(t)为与原始信号相关性最高的IMF成分；

当RRMSE达到最大值时所对应的噪声频率即设为本次CEEMD的最佳噪声上限频率f_p。

(3)根据步骤(2)所确定的噪声上限截止频率，对回波数据进行互补集总经验模态分解(CEEMD)，得到一组本征模态函数IMF；

具体过程如下：

(a)对回波时域信号x(t)分别施加一对相位相反、幅值相同的高斯白噪声ω(t)，设定噪声标准差为原信号标准差的0.01倍，噪声的上限截止频率为f_p，得到：

x⁺(t)＝x(t)+ω(t)

x^-(t)＝x(t)-ω(t)；

(b)分别对x+(t)和x-(t)进行经验模态分解运算，得到第i阶IMF分量c_i(t)，直到得到两组模态分解的结果，记为：c_i ⁺和c_i ^-，求其均值c_i＝(c_i ⁺+c_i ^-)/2。

(4)将所有IMF进行希尔伯特变换，求得瞬时频率，得出希尔伯特谱；

具体过程如下：

(a)对每个IMF分量c_i(t)作希尔伯特变换，可得数据序列

(b)将c_i(t)作为实部、作为虚部构成一个复序列解析信号z_i(t)：

其中，

(c)从而得到瞬时频率：

若忽略残余分量，则原信号x(t)即为z_i(t)的实部，表示为：

(d)以时间、瞬时频率为自变量，幅值在时频平面内以等高线表现，则可将三者的关系表示成一个三维图形，定义为HHT谱，记做H(ω,t)：

(5)根据步骤(4)的结果在希尔伯特谱中提取出与混凝土内部质量相关的特征。

参见图4至图8，为本实施例不同注浆密实度模型下，根据本发明的上述方法得出的HHT谱图，对比各图可发现，各图都出现了4k-6k的模态，这是对应底端界面的回波，随着注浆密实度越低，波形散射、反射越严重，对应HHT图可见，10KHz附近的频率越分散。因而本方法可以有效地检测出注浆密实度信息。

最后有必要在此说明的是：以上实施例只用于对本发明的技术方案作进一步详细的说明，不能理解为对本发明保护范围的限制，本领域的技术人员根据本发明的上述内容作出的一些非本质的改进和调整均属于本发明的保护范围。

Claims (4)

1.一种冲击回波数据处理方法，其特征在于包括如下步骤：

(1)利用冲击回波仪在待测混凝土外侧激发出应力波，并接收冲击回波数据；

(2)通过计算相对均方根误差极大值来确定所添加噪声的上限截止频率f_p；

(3)根据步骤(2)所确定的噪声上限截止频率，对回波数据进行互补集总经验模态分解(CEEMD)，得到一组本征模态函数IMF；

(4)将所有IMF进行希尔伯特变换，求得瞬时频率，得出希尔伯特谱；

(5)根据步骤(4)的结果在希尔伯特谱中提取出与混凝土内部质量相关的特征。

2.根据权利要求1所述的冲击回波数据处理方法，其特征在于：步骤(2)所述噪声的上限截止频率f_p的确定包括如下具体步骤：

(a)对于回波时域信号x(t)，设定一高斯白噪声ω(t)，其上限截止频率为f_p′＝m·f_s，m取10至20，f_s为x(t)的采样率，以频率f_p′为重采样率，将x(t)进行三次样条插值计算，得到新信号x′(t)；

(b)对信号x′(t)做加噪声的集总经验模态分解，集总次数为2次，所添加的噪声的上限频率为f_p'，噪声标准差为原信号标准差的0.01倍，在不同的噪声上限频率下计算RRMSE的值，公式如下：

$RRMSE=\sqrt{\underset{t=1}{\overset{N}{\Σ}}{\left(x\right(t)-c_{max}(t\left)\right)}^2/\underset{t=1}{\overset{N}{\Σ}}x{\left(t\right)}^2};$

其中，c_max(t)为与原始信号相关性最高的IMF成分；

当RRMSE达到最大值时所对应的噪声频率即设为本次CEEMD的最佳噪声上限截止频率f_p。

3.根据权利要求1所述的冲击回波数据处理方法，其特征在于：步骤(3)所述对回波数据进行互补集总经验模态分解(CEEMD)包括如下具体步骤：

(a)对回波时域信号x(t)分别施加一对相位相反、幅值相同的高斯白噪声ω(t)，设定噪声标准差为原信号标准差的0.01倍，噪声的上限截止频率为f_p，得到：

x⁺(t)＝x(t)+ω(t)

x^-(t)＝x(t)-ω(t)；

(b)分别对x⁺(t)和x^-(t)进行经验模态分解运算，得到第i阶IMF分量c_i(t)，直到得到两组模态分解的结果，记为：c_i ⁺和c_i ^-，求其均值c_i＝(c_i ⁺+c_i ^-)/2。

4.根据权利要求1所述的冲击回波数据处理方法，其特征在于：步骤(4)所述得出希尔伯特谱包括如下具体步骤：

(a)对每个IMF分量c_i(t)作希尔伯特变换，可得数据序列

(b)将c_i(t)作为实部、作为虚部构成一个复序列解析信号z_i(t)：

其中，

(c)从而得到瞬时频率：

若忽略残余分量，则原信号x(t)即为z_i(t)的实部，表示为：

$x\left(t\right)=\mathrm{Re}\[\underset{i=1}{\overset{n}{\Σ}}{a}_i\left(t\right)\exp (j\∫{\mathrm{\ω}}_i(t\left)dt\right)\];$

(d)以时间、瞬时频率为自变量，幅值在时频平面内以等高线表现，则可将三者的关系表示成一个三维图形，定义为HHT谱，记做H(ω,t)：

$H(\mathrm{\ω},t)=\underset{i=1}{\overset{n}{\Σ}}{a}_i\left(t\right)\left|\exp \left({\mathrm{j\ω}}_i\right(t\left)dt\right)\right|.$