一种基于变分模态分解的冲击回波信号分析方法与流程

本发明涉及一种基于变分模态分解的冲击回波信号分析方法，属于土木结构工程检测技术领域。

背景技术：

目前，采用冲击回波法对混凝土构件厚度与缺陷检测有了一定的研究，一般采用fourier变换求解出振动信号的频谱特征，再根据主频和表观速度计算出混凝土构件的厚度。但实际的冲击回波信号并非平稳信号，却是非平稳的信号，在实际工程中冲击回波不但从底板反射而且从缺陷处发生反射，使得fourier变换频谱常出现多重峰值，因而采用传统的fourier变换很难正确识别。hilber-huang变换将非平稳信号分解为有限个本征模态函数，然后对本征模态函数进行hilbert变换，再构造出信号的hilbert振幅谱，最后通过对hilbert谱积分得出边际谱，边际谱相对于fourier频谱更加准确反映信号的实际频谱，但是hilber-huang变换经验模态分解法常出现模态混叠和端点效应等问题。

技术实现要素：

本发明所要解决的技术问题是克服现有技术的缺陷，提供一种基于变分模态分解的冲击回波信号分析方法，根据变分模态分解和hibert变换方法原理，提出通过冲击回波信号的边际谱识别混凝土构件的缺陷的检测方法，该方法在强噪声影响下能很好的压制噪声干扰，相对于传统fourier变换具有更高的分辨率。

为解决上述技术问题，本发明提供一种基于变分模态分解的冲击回波信号分析方法，包括以下步骤：

1)根据现场测量环境要求，选择测量参数和设备工作参数，采集冲击回波信号；

2)根据冲击回波信号，设定变分模态分解中的参数；

3)利用变分模态分解方法将采集的冲击回波信号分解为n个不同尺度的本征模量函数；

4)利用hilbert变换，得出本征模量函数的hilbert时频谱；

5)将不同频率的hilbert时频谱在时间域进行积分得到最终的边际谱。

前述的步骤1)中，测量参数和设备工作参数包括采样频率，采样点及震源。

前述的步骤2)中，变分模态分解中的参数包括拉格朗日乘子及中心频率。

前述的步骤3)中，变分模态分解包括小波变换、短时傅立叶变换和经验模态分解。

前述的步骤3)的分解过程具体如下：

3-1)假设将采集的冲击回波信号分解为n个本征模量函数：

其中，f(t)为冲击回波信号，uk(t)表示第k个本征模量函数，t表示t时刻；

3-2)对每个本征模量函数进行hilbert变换得到解析信号：

其中，δ(t)是狄利克来函数，j是虚数符号；

3-3)通过移频方式，将各解析信号的频谱变换到基带上：

其中，ωk表示中心频率；

3-4)计算解析信号的欧式距离，计算各本征模量函数的带宽，变分问题如下：

3-5)利用二次罚函数项和lagrange乘子将变分问题转化为无约束问题：

其中，α为惩罚因子，λ(t)为lagrange乘子，表示范数，<·>表示内积运算；

3-6)利用乘子交替方向算法求取方程(5)的无约束变分问题，最终得到多个不同尺度的本征模量函数。

前述的步骤4)中，本征模量函数的hilbert时频谱表示为：

其中，p为柯西主值，τ为时间变量。

前述的步骤5)，求解边际谱，具体如下：

5-1)由本征模量函数的hilbert时频谱和本征模量函数构造出解析信号，形式为：

a[uk(t)]＝uk(t)+ih[uk(t)](7)

其中，a[uk(t)]表示解析信号；

5-2)求解解析信号的振幅和相位：

其中，ak(t)表示振幅，θk(t)表示相位；

5-3)对相位进行求导，计算得出解析信号的瞬时频率：

其中，ωk(t)表示瞬时频率。

5-4)将振幅ak(t)表示在联合的时频平面上，得出本征模量函数uk(t)的hilbert时频谱：

其中，hk(t,ω)表示hilbert时频谱，ω为频率；

5-5)对冲击回波信号f(t)进行整体hilbert谱分析，hilbert时频谱表示为：

5-6)根据式(12)的hilbert时频谱进一步定义出边际谱：

其中，h(ω)表示边际谱。

本发明所达到的有益效果：

1)本发明针对非稳定频率的冲击回波信号，采用边际谱识别混凝土构件的缺陷，在强噪声影响下能很好的压制噪声干扰，相对于传统fourier变换具有更高的分辨率。

2)本发明首次提出用变分模态分解法对冲击回波信号进行分析，该方法相对现有的多尺度分析方法有更完善的数学基础及实际效果，该方法采用温纳滤波理论不仅能分解出不同尺度多模态，而且能压制噪声。

3)本发明并不仅限于对冲击回波信号的特征分析，也可以用于其它振动检测信号分析，例如桩基低应变检测信号的分析。

附图说明

图1为本发明方法的流程图；

图2为冲击回波检测混凝土底板几何模型；

图3为底板模拟信号，(a)为原数值模拟信号，(b)为加入白噪声的强干扰信号；

图4为对强干扰信号处理后的频谱图，(a)为fourier变换频谱，(b)为vmd边际频谱。

具体实施方式

下面对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明的基于变分模态分解的冲击回波信号分析方法，如图1所示，包括以下步骤：

1、根据现场测量环境要求，选择合适的测量参数和设备工作参数，如设置采样频率，采样点及震源，采集冲击回波信号。

2、根据冲击回波信号，设定变分模态分解中的参数；参数包括拉格朗日乘子及中心频率。

3、利用变分模态分解方法将采集的冲击回波信号分解为若干个不同尺度的本征模量函数，求解变分问题，具体如下：

3-1)假设将冲击回波信号分解为n个本征模量函数：

其中，f(t)为冲击回波信号，uk(t)表示第k个本征模量函数，t表示t时刻；

变分模态分解，可以是传统的信号分析方法，包括小波变换、短时傅立叶变换和经验模态分解等方法。

3-2)对每个本征模量函数进行hilbert变换得到其解析信号：

其中，δ(t)是狄利克来函数，j是虚数符号；

3-3)通过移频方式，将各解析信号的频谱变换到基带上：

其中，ωk表示中心频率；

3-4)计算上述解析信号的欧式距离，估计各模态带宽，其变分问题如下：

3-5)利用二次罚函数项和lagrange乘子将其转化为无约束问题：

式中，α为惩罚因子，λ(t)为lagrange乘子，表示范数，<·>表示内积运算；

3-6)再利用乘子交替方向算法求取方程(5)无约束变分问题，最终得到多个不同尺度的本征模量函数。

4、利用hilbert变换，得出本征模量函数的hilbert时频谱；变换形式表示为：

其中，p为柯西主值，τ为时间变量。

5、将不同频率的hilbert时频谱在时间域进行积分得到最终的边际谱，具体过程如下：

5-1)由本征模量函数的hilbert时频谱式(6)和本征模量函数构造出解析信号，其形式为：

a[uk(t)]＝uk(t)+ih[uk(t)](7)

a[uk(t)]表示解析信号。

5-2)对公式(7)的解析信号，求出其振幅和相位：

其中，ak(t)表示振幅，θk(t)表示相位。

5-3)对相位进行求导，计算得出解析信号的瞬时频率：

其中，ωk(t)表示瞬时频率。

5-4)将振幅ak(t)表示在联合的时频平面上，得出本征模量函数uk(t)的hilbert时频谱：

其中，hk(t,ω)表示hilbert时频谱，ω为频率。

5-5)对f(t)进行整体hilbert谱分析，hilbert时频谱可表示为：

5-6)根据式(12)的hilbert时频谱进一步定义出边际谱：

其中，h(ω)表示边际谱。

实施例

冲击回波是基于瞬态p波进行的结构无损检测，构件的厚度按下列公式计算：

其中，t为混凝土结构构件的厚度，vp为p波的波速，ωc为振幅谱主频。

p波波速vp与结构的杨氏模量e，泊松比ν及结构密度ρ存在下列关系：

以一个混凝土底板作为本次试验对象，混凝土底板尺寸1.50m×1.50m，厚度为0.20m。杨氏模量e，泊松比ν及结构密度ρ分别为37.8gpa，2500kg/m3，0.2。根据公式(15)可得出纵波波速4099m/s。冲击回波厚度根据公式(14)计算，fourier主频理论值应为10247hz。

对于冲击回波表面响应加速度值，根据模型参数通过有限元软件进行仿真。试验中将混凝土底板中心作为冲击球捶击点，冲击点与加速度接收点如图2所示。图中，带箭头为冲击点，另一点为加速度接收点。冲击球的冲击点与加速度检波器之间距离为0.06m，加速度采样频率取1mhz，采样长度取2ms。冲击球荷载的由下列函数确定：

p(t)＝(sin(πt/40))³0≤t≤40μs(16)

通过上述参数设置，模拟出底板表面接收点的加速度随时间的变化情况，如图3(a)所示，取初始信号的信噪比snr＝-8，将原数值模拟信号图3(a)加入白噪声，得到加入强干扰后的结果(图3(b))。

将含有强干扰的信号分别利用fourier和vmd变换，求解出fourier频谱和边际谱，变换结果如图3所示，图4(a)为fourier变换频谱，图4(b)为vmd边际频谱。

根据图4(a)分析结果，在强白噪声背景下采用fourier变换得到的频谱图中显示多个振动主峰，呈锯齿状，不包含单一的主频率峰值，主频率特征不明显。而在图4(b)采用vmd方法得出的边际谱，背景噪声对主频峰值的识别影响较小，即使在强白噪声背景下主频峰值也很明显。因此，vmd方法得出的边际谱能很好的压制干扰的影响，相对传统的fourier频谱具有显著优势。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。