基于频谱分解和相谱增强方法的兰姆波相速度频散评估

Phase Velocity Dispersion Evaluation of Lamb Wave Based on Spectrum Decomposition and Phase Spectrum Enhancement Methods

摘要

兰姆波因其传播距离长、 衰减小、 对材料缺陷和特性灵敏度高等优势， 被广泛应用于无损检测和结构健康监测中。相速度是兰姆波的基本特性之一， 可用于缺陷检测、 材料性能估计和频散评估等。现有的相速度频散重建方法大多需要布置多个传感器阵列， 成本高且布置复杂。鉴于此， 提出了一种基于频谱分解和相谱增强的混合方法， 仅需两个相邻传感器即可重建兰姆波的相速度频散曲线。在铝板和碳纤维复合材料板中进行了仿真验证， 并在碳纤维复合材料板中进行了实验验证， 最后与半解析有限元法的理论结果进行了比较。研究结果表明， 兰姆波A0模态和S0模态的相速度频散曲线均被成功重建， 且重建结果与理论值吻合度达到90%以上。该方法不仅简化了测量过程， 还提高了频散曲线重建的精确度和稳定性， 为无损检测和结构健康监测提供了新的技术手段。

Abstract

Lamb waves are widely used in non-destructive testing and structural health monitoring due to their advantages of long propagation distance， low attenuation， and high sensitivity to material defects and properties. Phase velocity， as one of the fundamental characteristics of Lamb waves， can be utilized for defect detection， material property estimation， and dispersion evaluation. Most existing phase velocity dispersion reconstruction methods require the deployment of multiple sensor arrays， which are costly and complex to arrange. This paper proposes a hybrid method based on spectrum decomposition and phase spectrum enhancement， which reconstructs the phase velocity dispersion curves of Lamb waves using only two adjacent sensors. The proposed method was validated through simulations on both aluminum and CFRP plates， and further verified experimentally on CFRP plates. The results were compared with theoretical values obtained from the semi-analytical finite element method. The findings demonstrate that the phase velocity dispersion curves of both the A0 and S0 modes of Lamb waves were successfully reconstructed， with a consistency rate of over 90% compared to the theoretical values. This method not only simplifies the measurement process but also enhances the accuracy and stability of dispersion curve reconstruction， providing a new technical approach for non-destructive testing and structural health monitoring.

Graphical abstract

关键词

兰姆波 / 相速度 / 频散重建 / 频谱分解 / 相谱增强

Key words

Lamb wave / phase velocity / dispersion reconstruction / spectrum decomposition / phase spectrum enhancement

0 引　言

兰姆波作为一种在固体板中传播的弹性波， 因其长传播距离、 低衰减以及对材料缺陷和特性具有高灵敏度的显著优势， 在无损检测和结构健康监测领域中受到广泛关注^［1-3］。兰姆波能够穿透结构材料检测材料内部的缺陷， 如裂纹、 腐蚀和脱粘等， 不会对结构造成任何损害。此外， 兰姆波的传播特性使其能够提供关于材料状态和结构完整性的重要信息， 这对于预测潜在的结构失效和制定维护策略至关重要。

然而， 兰姆波的应用也面临着一些挑战。兰姆波在传播过程中表现出明显的频散性和多模态性， 这意味着不同频率的波将以不同的速度传播， 且在结构中可能存在多种模式的波。这些特性使得兰姆波的接收信号复杂， 难以解释^［4‐5］。频散现象是兰姆波的一个主要特征， 它不仅影响波信号的分析和处理， 还增加了信号扩展的持续时间， 并根据距离降低信号幅度^［6‐7］。因此， 对于应用中使用的兰姆波， 应该使用能够评估频散现象并将频率与相速度和群速度联系起来的信号处理技术。基于信号幅度或功率谱测量的传统信号处理方法不能解决这些任务， 需要新的信号处理技术来估计相速度值， 将它们与频率联系起来， 并重建频散曲线。

近年来， 研究人员提出了一系列基于时频分析的信号处理方法来确定频散曲线， 这些方法包括短时傅立叶变换、 小波变换、 Wigner-Ville分布和Choi-Williams分布等。然而， 这些方法在分析兰姆波信号时遇到了一些挑战， 包括时频分辨率的极限和内部干扰交叉项等问题^［8‐9］。二维快速傅立叶变换（Two-Dimensional Fast Fourier Transform， 2D-FFT）是相速度频散估计中最广泛使用的方法之一， 它在时间和空间上执行双重FFT^{［10‐12］}。尽管这种方法在某些情况下是有效的， 但也有一些已知的限制。例如， 它需要沿着波传播方向进行多点信号采集， 以重建频散曲线的片段， 这在实际操作中是繁琐且复杂的^［13］。因此， 有必要寻求新的方法和途径以更简单、 更准确地分析和评估非平稳兰姆波信号。最近提出的重建相速度频散曲线的解决方案之一是仅使用两个测量信号。Zeng等^［13］提出了使用时变沃尔德卡尔曼滤波器的短时间线性调频傅立叶变换方法。然而， 这种方法在兰姆波的生成中使用了啁啾信号， 这意味着需要同时在时域和频域中实现高分辨率， 算法复杂， 并且对于相速度频散曲线的重建不够精确和稳定。

针对上述挑战， 本文提出了一种基于频谱分解和相谱增强的混合方法。与现有方法相比， 本文方法的技术创新点主要体现在以下方面： 该方法将频谱分解技术与相谱增强技术相结合， 通过频谱分解提高小幅值分量的灵敏度， 并通过相谱增强精确估计相位差， 从而实现了高精度的相速度重建； 传统方法通常需要多个传感器阵列， 而本文方法仅需两个相邻传感器， 显著降低了实验成本和复杂性； 通过带通高斯滤波器的设计， 本文方法能够在宽频范围内重建相速度频散曲线， 特别是在低频区域的相速度估计上表现出较高的精度； 本文方法不仅在各向同性铝板中得到了验证， 还在各向异性碳纤维复合材料（Carbon Fiber Reinforced Polymers， CFRP）板中展现了良好的适用性， 与半解析有限元法（Semi-Analytical Finite Element Method， SAFE）的结果进行了比较， 吻合度在90%以上。通过这些研究， 本文旨在为兰姆波相速度频散评估提供一个新颖、 有效的解决方案， 以促进其在无损检测和结构健康监测领域的应用。

1 兰姆波传播模型

1.1 复合材料的弹性模量

通常， 复合材料的弹性常数共有36个独立分量， 可由一个\mathrm{6}\times \mathrm{6}的矩阵\mathit{C}表示。本文所研究的CFRP具备正交对称性， 根据对称等价条件， 其有效弹性模量可以简化为9个。依据胡克定律， 应力、 应变与弹性模量之间的关系可以表述为

\begin{array}{c}\left[\begin{array}{c}\begin{array}{l}{\sigma }_{\mathrm{11}}\\ {\sigma }_{\mathrm{22}}\\ {\sigma }_{\mathrm{33}}\end{array}\\ \begin{array}{l}{\sigma }_{\mathrm{23}}\\ {\sigma }_{\mathrm{13}}\\ {\sigma }_{\mathrm{12}}\end{array}\end{array}\right]=\left[\begin{array}{cccccc}{C}_{\mathrm{11}}& C_{\mathrm{12}}& C_{\mathrm{13}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ C_{\mathrm{12}}& C_{\mathrm{22}}& C_{\mathrm{23}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ C_{\mathrm{13}}& C_{\mathrm{23}}& C_{\mathrm{33}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& C_{\mathrm{44}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& C_{\mathrm{55}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& C_{\mathrm{66}}\end{array}\right]\left[\begin{array}{c}\begin{array}{l}{\epsilon }_{\mathrm{11}}\\ {\epsilon }_{\mathrm{22}}\\ {\epsilon }_{\mathrm{33}}\end{array}\\ \begin{array}{l}\mathrm{2}{\epsilon }_{\mathrm{23}}\\ \mathrm{2}{\epsilon }_{\mathrm{13}}\\ \mathrm{2}{\epsilon }_{\mathrm{12}}\end{array}\end{array}\right],\end{array}（1）

式中： {\sigma }_{ij}为应力； {\epsilon }_{ij}为应变； C_{ij}为弹性矩阵\mathit{C}的元素， 通常为复数形式， 实部为弹性矩阵， 虚部为粘滞矩阵， 可以解释兰姆波在复合材料中的衰减。为了简化模型， 本文中的C_{ij}只取实部， 忽略了复合材料的粘弹性。

需要注意的是， 式（1）描述的是材料坐标系下的弹性常数矩阵， 而材料坐标系与空间全局坐标系并不重合。因此， 在计算时， 需要基于复合材料层的纤维坐向角度旋转坐标系， 得到统一旋转坐标系的正确表示。假设复合材料层合板相对于全局坐标系的单层纤维坐向角度为\theta， 那么正确全局表示可以通过以下变换关系旋转每层局部坐标系得到。

{\mathit{C}}_{\theta }={\mathit{R}}_{\mathrm{1}}\mathit{C}{\mathit{R}}_{\mathrm{2}}^{-\mathrm{1}},（2）

\begin{array}{c}{\mathit{R}}_{\mathrm{1}}=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& c^{\mathrm{2}}& s^{\mathrm{2}}& -\mathrm{2}cs& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& s^{\mathrm{2}}& c^{\mathrm{2}}& \mathrm{2}cs& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& cs& -cs& c^{\mathrm{2}}-s^{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& c& s\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -s& c\end{array}\right],\end{array}（3）

\begin{array}{c}{\mathit{R}}_{\mathrm{2}}=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& c^{\mathrm{2}}& s^{\mathrm{2}}& -cs& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& s^{\mathrm{2}}& c^{\mathrm{2}}& cs& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}cs& -\mathrm{2}cs& c^{\mathrm{2}}-s^{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& c& s\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -s& c\end{array}\right],\end{array}（4）

式中： s=\mathrm{s}\mathrm{i}\mathrm{n}\theta ，c=\mathrm{c}\mathrm{o}\mathrm{s}\theta。

经坐标变换后的弹性矩阵{\mathit{C}}_{\theta }可以用于复合材料板内兰姆波传播模型的建模， 需要注意的是， CFRP板内每一层的铺层角度\theta可能是不同的， 需要通过独立的矩阵变换到空间坐标系下计算。

1.2 兰姆波半解析有限元模型

半解析有限元法在兰姆波的传播方向上使用波动方程求解， 在板的厚度方向上使用一维谱单元代替横截面的二维近似， 一个谱单元通常包含3~5个节点， 且可以非均匀分布， 每个节点有3个自由度^［14］。

SAFE方法建模的过程可以概括为： 首先， 使用一系列谱单元离散化垂直于兰姆波传播方向的板截面， 并结合单元最小位能原理， 推导单元刚度矩阵和质量矩阵； 然后， 求解总刚度矩阵和质量矩阵， 得到其特征方程； 最后， 通过求解特征方程的特征值和特征向量， 得到固定条件下兰姆波频散特性分布曲线的数值解。

兰姆波在板内传播的SAFE模型如图 1 所示。

对于任意的传播方向角\beta， 对应的波矢量\mathit{k}可以表示为

\mathit{k}=k\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \right)\widehat{\mathit{z}}-k\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta \right)\widehat{\mathit{y}},（5）

式中： \widehat{\mathit{z}}和\widehat{\mathit{y}}为坐标轴方向的单位矢量； k为对应的波数。

通常， 一个谱单元e的位移场\mathit{u}是时间t和空间(x,y,z)的函数，

\mathit{u}\left(x,y,z,t\right)=\mathit{N}{\mathit{q}}^{\left(e\right)}{\mathrm{e}}^{\mathrm{j}\left[\omega t+k\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta \right)y-k\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \right)z\right]},（6）

\mathit{N}=

\left[\begin{array}{cccccccccc}{\phi }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& {\phi }_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \dots & {\phi }_n& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\phi }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& {\phi }_{\mathrm{2}}& \mathrm{0}& \dots & \mathrm{0}& {\phi }_n& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\phi }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& {\phi }_{\mathrm{2}}& \dots & \mathrm{0}& \mathrm{0}& {\phi }_n\end{array}\right],（7）

式中： \mathit{N}为由形函数{\phi }_i组成的插值矩阵， 此处{\phi }_i采用Gauss-Lobatto-Legendre节点的拉格朗日插值多项式^［15］； n为谱单元内节点数量； \mathit{q}为包含节点位移矢量的列向量。

谱单元e的应变场\mathit{\epsilon }与其位移场\mathit{u}有关，

\mathit{\epsilon }={\mathit{L}}_{\mathit{x}}\frac{\partial \mathit{u}}{\partial x}+{\mathit{L}}_{\mathit{y}}\frac{\partial \mathit{u}}{\partial y}+{\mathit{L}}_{\mathit{z}}\frac{\partial \mathit{u}}{\partial z},（8）

{\mathit{L}}_{\mathit{x}}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right],{\mathit{L}}_{\mathit{y}}=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right],

{\mathit{L}}_{\mathit{z}}=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\mathrm{。}（9）

将式（8）所示的应变-位移关系代入式（6）进一步推导可得

\mathit{\epsilon }=\left[{\mathit{B}}_{\mathrm{1}}-\mathrm{j}{k}_y{\mathit{B}}_{\mathrm{2}}-\mathrm{j}{k}_z{\mathit{B}}_{\mathrm{3}}\right]{\mathit{q}}^{\left(e\right)}{\mathrm{e}}^{\mathrm{j}\left[\omega t+k\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta \right)y-k\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \right)z\right]},（10）

式中： {\mathit{B}}_{\mathrm{1}}={\mathit{L}}_{\mathit{x}}{\mathit{N}}_{\mathit{x}}，{\mathit{B}}_{\mathrm{2}}={\mathit{L}}_{\mathit{y}}\mathit{N}，{\mathit{B}}_{\mathrm{3}}={\mathit{L}}_{\mathit{z}}\mathit{N}；{\mathit{N}}_{\mathit{x}}为\mathit{N}对x的导数。

基于材料的弹性矩阵{\mathit{C}}_{\theta }可以得到SAFE模型的质量矩阵\mathit{M}和刚度矩阵{\mathit{K}}_{mn}。

{\mathit{m}}^{\left(e\right)}=\underset{\left(e\right)}{\int }{N}^{\mathrm{T}}\rho N\mathrm{d}x,\mathit{M}=\stackrel{n_e}{\underset{e=\mathrm{1}}{\cup }}{\mathit{m}}^{\left(e\right)},（11）

{\mathit{k}}_{mn}^{\left(e\right)}=\underset{\left(e\right)}{\int }{\mathit{B}}_m^T{\mathit{C}}_{\theta }{\mathit{B}}_{n}dx,{\mathit{K}}_{mn}=\stackrel{n_e}{\underset{e=\mathrm{1}}{\cup }}{\mathit{k}}_{mn}^{\left(e\right)},（12）

式中： \rho为材料密度； n_e为谱单元数量。

基于质量矩阵\mathit{M}和刚度矩阵{\mathit{K}}_{mn}， 式（10）可以进一步改写成特征方程的形式。

\left[\mathit{A}-\gamma \mathit{B}\right]\mathit{V}=\mathrm{0},（13）

\begin{array}{c}\mathit{A}=\left[\begin{array}{cc}\mathrm{0}& {\mathit{K}}_{\mathrm{11}}-{\omega }^{\mathrm{2}}\mathit{M}\\ {\mathit{K}}_{\mathrm{11}}-{\omega }^{\mathrm{2}}\mathit{M}& \begin{array}{l}-cc{\mathit{K}}_{\mathrm{13}}-ss{\mathit{K}}_{\mathrm{21}}+\\ ss{\mathit{K}}_{\mathrm{12}}+cc{\mathit{K}}_{\mathrm{31}}\end{array}\end{array}\right],\end{array}（14）

\begin{array}{c}\mathit{B}=\left[\begin{array}{cc}{\mathit{K}}_{\mathrm{11}}-{\omega }^{\mathrm{2}}\mathit{M}& \mathrm{0}\\ \mathrm{0}& \begin{array}{l}s{s}^{\mathrm{2}}{\mathit{K}}_{\mathrm{22}}+c{c}^{\mathrm{2}}{\mathit{K}}_{\mathrm{33}}-\\ ccss{\mathit{K}}_{\mathrm{23}}-ccss{\mathit{K}}_{\mathrm{32}}\end{array}\end{array}\right],\end{array}（15）

\begin{array}{c}\mathit{V}={\left[\begin{array}{cc}\mathit{q}& \gamma \mathit{q}\end{array}\right]}^{\mathrm{T}},\end{array}（16）

式中： \gamma =\mathrm{i}k，ss=\mathrm{s}\mathrm{i}\mathrm{n}\beta， cc=\mathrm{c}\mathrm{o}\mathrm{s}\beta。

给定材料相关参数后， 通过求解全局控制方程， 获得导波全部模态的相应频率范围内的波数变化情况， 进而得到理论相速度， 即可获得兰姆波传播的理论相速度频散曲线， 具体公式为

c_p=\frac{\omega }{k}\mathrm{。}（17）

2 兰姆波相速度频散重建方法

所提相速度重建方法采用经典的相移方法来估计与接收信号幅度谱峰值相对应的特定频率处的速度值， 同时结合频谱分解方法， 提高小幅值分量的灵敏度， 从而实现在宽带上重建速度值。所提算法的一般流程图及其主要步骤如图 2 所示， 可以概括为以下步骤：

1） 测量获得在特定距离x_{\mathrm{1}}和x_{\mathrm{2}}处的两个相邻兰姆波信号u_{x\mathrm{1}}(t)和u_{x\mathrm{2}}(t)。

2） 测量信号u_{x\mathrm{1}}(t)和u_{x\mathrm{2}}(t)的频谱U_{x\mathrm{1}}(f)和U_{x\mathrm{2}}(f)分别为

U_{x\mathrm{1}}\left(f\right)=FT\left[u_{x\mathrm{1}}\left(t\right)\right],U_{x\mathrm{2}}\left(f\right)=FT\left[u_{x\mathrm{2}}\left(t\right)\right],（18）

式中： FT表示傅里叶变换； U_{x\mathrm{1}}\left(f\right)\mathrm{、}{U}_{x\mathrm{2}}\left(f\right)分别为u_{x\mathrm{1}}(t)\mathrm{、}{u}_{x\mathrm{2}}(t)的频谱。

3） 频谱由具有预定义参数的k个带通高斯滤波器进行滤波。

\begin{array}{c}{U}_{x\mathrm{1},k}\left(f\right)=U_{x\mathrm{1}}\left(f\right)\cdot B_k\left(f\right),\\ U_{x\mathrm{2},k}\left(f\right)=U_{x\mathrm{2}}\left(f\right)\cdot B_k\left(f\right),\end{array}（19）

式中： U_{x\mathrm{1},k}\left(f\right)和U_{x\mathrm{2},k}\left(f\right)分别为分解后的子频谱； B_k\left(f\right)={\mathrm{e}}^{\mathrm{4}\mathrm{l}\mathrm{n}\left(\mathrm{0.5}\right){\left(\frac{f-f_L-\left(k-\mathrm{1}\right)df}{∆B}\right)}^{\mathrm{2}}}为第k个带通滤波器的频率响应； k=\mathrm{1,2},\dots ,K， K为滤波器总数； f_L和f_H定义了滤波器的中心频率滑动区间； ∆B为滤波器带宽； df=\left(f_H-f_L\right)/(K-\mathrm{1})为两个相邻滤波器中心频率之间的频域移动步长^［16］。

4） 使用傅里叶逆变换重建滤波后的信号，

\begin{array}{c}{u}_{x\mathrm{1},k}\left(t\right)=IFT\left(U_{x\mathrm{1},k}\left(f\right)\right),\\ u_{x\mathrm{2},k}\left(t\right)=IFT\left(U_{x\mathrm{2},k}\left(f\right)\right),\end{array}（20）

式中： IFT表示傅里叶逆变换。

5） 为了避免相位展开过程中的不确定性， 每个波形u_{x\mathrm{1},k}\left(t\right)和u_{x\mathrm{2},k}\left(t\right)在时域中偏移t_{m\mathrm{1}}和t_{m\mathrm{2}}。在中等频散情况下， 可以根据希尔伯特包络的最大值来偏移波形。

e_{x\mathrm{1},k}\left(t\right)=u_{x\mathrm{1},k}\left(t+t_{m\mathrm{1}}\right),e_{x\mathrm{2},k}\left(t\right)=u_{x\mathrm{2},k}\left(t+t_{m\mathrm{2}}\right),（21）

\begin{array}{c}{t}_{m\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\left(\underset{t}{\mathrm{m}\mathrm{a}\mathrm{x}} \left[HT\left|u_{x\mathrm{1},k}\left(t\right)\right|\right]\right),\\ t_{m\mathrm{2}}=\mathrm{a}\mathrm{r}\mathrm{g}\left(\underset{t}{\mathrm{m}\mathrm{a}\mathrm{x}} \left[HT\left|u_{x\mathrm{2},k}\left(t\right)\right|\right]\right),\end{array}（22）

式中： HT表示希尔伯特变换； t_{m\mathrm{1}}和t_{m\mathrm{2}}为偏移时间， 其对应于希尔伯特包络的最大值， 以补偿由于相速度引起的信号延迟。

6） 时移后信号的复频谱计算为

\begin{array}{c}{E}_{x\mathrm{1},k}\left(f\right)=FT\left[e_{x\mathrm{1},k}\left(t\right)\right],\\ E_{x\mathrm{2},k}\left(f\right)=FT\left[e_{x\mathrm{2},k}\left(t\right)\right]\mathrm{。}\end{array}（23）

7） 对于给定的频带f， 估计偏移信号e_{x\mathrm{1},k}\left(t\right)和e_{x\mathrm{2},k}\left(t\right)之间的相位差\Delta {\phi }_k(f)为

\Delta {\phi }_k\left(f\right)=\left({\alpha }_{x\mathrm{1},k}\left(f\right)-{\alpha }_{x\mathrm{2},k}\left(f\right)\right)，（24）

\begin{array}{c}{\alpha }_{x\mathrm{1},k}\left(f\right)=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left[\frac{\mathrm{I}\mathrm{m}\left[E_{x\mathrm{1},k}\left(f\right)\right]}{\mathrm{R}\mathrm{e}\left[E_{x\mathrm{1},k}\left(f\right)\right]}\right],\\ {\alpha }_{x\mathrm{2},k}\left(f\right)=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left[\frac{\mathrm{I}\mathrm{m}\left[E_{x\mathrm{2},k}\left(f\right)\right]}{\mathrm{R}\mathrm{e}\left[E_{x\mathrm{2},k}\left(f\right)\right]}\right],\end{array}（25）

式中： \mathrm{I}\mathrm{m}和\mathrm{R}\mathrm{e}分别为复频谱的实部和虚部。注意到， 相位{\alpha }_{x\mathrm{1},k}\left(f\right)和{\alpha }_{x\mathrm{2},k}\left(f\right)是在\left[-\mathrm{\pi },\mathrm{\pi }\right]范围内计算的， 这意味着， 如果相位超过\pm \mathrm{\pi }弧度的极限， 将出现一些不连续性。因此， 相位{\alpha }_{x\mathrm{1},k}\left(f\right)和{\alpha }_{x\mathrm{2},k}\left(f\right)必须被展开。

8） 相速度和频率通过式（26）选取和计算。

c_p\left(k\right)=\frac{\mathrm{2}\mathrm{\pi }{f}_{k}d}{\Delta {\phi }_k\left(f_k\right)-\mathrm{2}\mathrm{\pi }{f}_k\left(t_{m\mathrm{1}}-t_{m\mathrm{2}}\right)},（26）

式中： f_k为对应于幅度谱|E_{x\mathrm{1},k}\left(f\right)|的峰峰值频率； k=\mathrm{1,2},\dots ,K， K为滤波器总数量； d为两传感器之间的距离， d=x_{\mathrm{2}}-x_{\mathrm{1}}。

最终， 重构的相速度频散曲线段被定义为计算所得频率点和相速度点\left(f_k,c_p\left(k\right)\right)的集合。

3 仿真验证

为了验证所提出方法的有效性， 在COMSOL中分别建立了各向同性铝板和各向异性CFRP板的仿真模型， 其结构和空间坐标系如图 3 所示。致动器和传感器均使用压电材料PZT-5H制成， 建模为长5 mm、 厚0.5 mm的矩形， 使用轻薄的线弹性层与板粘结在一起； 致动器和传感器下表面均设置为接地， 致动器上表面施加激励信号， 激励信号为一高斯窗调制的余弦脉冲， 表达式为

g_{\mathrm{0}}\left(t\right)=A{\mathrm{e}}^{\frac{{\left(t-t_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\mathrm{\pi }{f}_c\left(t-t_{\mathrm{0}}\right)\right),（27）

式中： A=1； f_c=200 kHz； t_{\mathrm{0}}=\mathrm{12.5}\mathrm{\mu }\mathrm{s}； \sigma=\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{6}}。

激励信号在时域和频域上的表示如图 4 所示。

3.1 各向同性铝板验证

在COMSOL中建立各向同性铝板模型， 其厚度为4 mm， 材料的力学性能如表 1 所示， 两个压电元件传感器分别以距离致动器350 mm和360 mm连接到铝板上表面。COMSOL仿真记录的导波时域信号如图 5（a） 所示， A0模态和S0模态相速度频散曲线重建的结果如图 5（b） 所示。从图中可以看出， A0模态和S0模态的相速度频散曲线均被成功重建， 且与SAFE法求得的理论结果吻合良好。

3.2 各向异性CFRP板验证

建立铺层为［0/90］_4s的CFRP层合板模型， 其厚度为2.3 mm， 每一层碳纤维都是单独建模后形成联合体。材料的力学性能如表 2 所示， 两个压电元件传感器分别以距离致动器227 mm和234 mm连接到复合材料层合板上表面。COMSOL仿真记录的导波时域信号如图 6（a） 所示， A0模态和S0模态相速度频散曲线重建结果如图 6（b） 所示。从图中可以看出， 即使在复杂的各向异性复合材料中， 本文方法仍能有效重建A0模态和S0模态的相速度频散曲线， 表明本文方法对于各向异性材料仍然有效且具有较高精度。

3.3 与传统方法的对比分析

为了进一步体现所提方法的优越性， 本文将其与传统2D-FFT方法进行了详细的对比分析。本文方法在多个方面均优于传统2D-FFT方法： 在传感器数量方面， 2D-FFT方法需要多个传感器阵列， 而本文方法仅需2个传感器， 显著降低了实验成本和复杂性； 在重建频率范围方面， 本文方法可实现对相速度频散曲线的宽带重建， 而2D-FFT方法一般为窄带重建； 在计算精度方面， 本文方法计算精度高， 而2D-FFT方法计算精度较低； 在计算效率方面， 本文方法的计算效率显著高于2D-FFT方法。

4 实验验证

为了证明所提方法在实际应用中的有效性， 采用的实验样件是在航空航天中普遍应用的 ［0/90］_4s铺层碳/环氧树脂层合板。其各项力学性能与仿真验证所用复合材料板一致， 尺寸为1 000 mm×1 000 mm， 厚度为（2.3±0.1） mm， 密度为1 630.4 kg/m³。

实验所用仪器如下：

函数/任意波形发生器： 型号为Tek AFG3052C， 用于产生目标激励电压信号。该发生器具有高精度和宽频带特性， 能够生成复杂的波形信号。

宽频功率放大器： 型号为Aigtek ATA-2022H， 用于放大波形发生器产生的激励信号， 以确保能激发压电致动器。

数字示波器： 型号为RIGOLL MSO8104四通道数字示波器， 用于记录传感器接收到的超声导波信号。其高采样率和带宽确保了信号的高保真度。

压电致动器和传感器： 采用PZT-5H材料制成， 直径为5 mm， 厚度为1 mm。致动器和传感器通过一层薄的环氧树脂胶与复合材料板粘结， 确保良好的声学耦合。

电脑配套装置： 用于数据采集和分析。

图 7 为搭建的复合材料层合板超声导波检测平台。波形发生器用于产生目标激励电压信号， 随后， 电压信号被放大并传输到压电致动器PZT； 通过压电效应， PZT将产生超声导波并随后被压电传感器PZT检测到； 最后， 使用数字示波器记录检测到的超声导波。

传播方向选择为0°方向， 激励信号与仿真验证时所用激励信号一致， 为高斯窗调制的余弦脉冲。经过功率放大器将电压放大， 再以示波器同时记录距离220 mm和230 mm处传感器PZT接收到的导波信号波形， 采样频率设置为10 MHz。

记录导波信号后， 执行兰姆波频散曲线重建算法， 结果如图 8 所示。分析实验结果可以看出， 重建的散点与理论频散曲线基本一致， A0模态和S0模态的相速度频散曲线被成功重建。

分析实验的误差来源， 可能有传感器粘贴位置的微小偏差导致的相位差计算误差， 及信号中存在的噪声影响相位差的精确计算。

为了验证实验的可重复性， 本文在相同实验条件下进行了多次重复实验， 每次实验均重新粘贴压电致动器和传感器， 并重新设置仪器参数。实验结果表明， 重建的相速度频散曲线在不同次实验中的一致性较高， 重建结果与理论频散曲线的吻合度均在90%以上， 表明所提方法具有较好的稳定性和可重复性。

此外， 本文还对不同位置的传感器进行了实验， 以验证方法的空间适用性。结果表明， 在不同位置重建的相速度频散曲线也具有较高的一致性， 进一步证明了该方法的鲁棒性。

5 结　论

本文提出了一种基于频谱分解和相谱增强方法的兰姆波相速度频散曲线重建技术， 并通过仿真和实验验证了其有效性。主要结论如下： 1） 所提方法能够成功重建兰姆波A0模态和S0模态的相速度频散曲线， 且重建结果与半解析有限元法的理论值吻合良好， 吻合度达到90%以上； 2） 与传统方法相比， 本文方法仅需两个相邻传感器即可完成频散曲线的重建， 显著简化了测量过程， 降低了实验成本和复杂性； 3） 通过频谱分解和相谱增强技术的结合， 本文方法提高了小幅值分量的灵敏度， 实现了宽带范围内的相速度重建， 具有较高的精度和稳定性； 4） 该方法在各向同性铝板和各向异性CFRP板中均得到了验证， 表明其适用于不同材料的兰姆波频散曲线重建； 5） 通过多次重复实验和不同位置的测量， 验证了该方法具有良好的可重复性和空间适用性。

尽管本文方法在铝板和CFRP板上展现了良好的有效性， 但其在更广泛实际工程应用场景中的普适性还有待进一步研究验证。未来的工作将集中在以下方面： 1） 传播方向的影响： 研究传播方向对兰姆波相速度频散曲线重建的影响， 验证本文方法在不同传播方向下的有效性； 2） 复杂工况的验证： 在实际工程应用中， 材料可能会受到温度、 湿度、 应力等复杂工况的影响， 未来的研究将探究这些复杂工况对本文方法有效性的影响， 并进一步优化算法以提高其鲁棒性。

参考文献

[1] XU C B， YANG Z B， CHEN X F， et al. A guided wave dispersion compensation method based on compressed sensing［J］. Mechanical Systems and Signal Processing， 2018， 103： 89-104.

[2] VOß M， ILSE D， HILLGER W， et al. Numerical simulation of the propagation of Lamb waves and their interaction with defects in CFRP laminates for non-destructive testing［J］. Advanced Composite Materials， 2020， 29（5）： 423-441.

[3] DU Y， ZHOU S， JING X， et al. Damage detection techniques for wind turbine blades： a review［J］. Mechanical Systems and Signal Processing， 2020， 141： 106445.

[4] ZUO H， YANG Z， XU C， et al. Damage identification for plate-like structures using ultrasonic guided wave based on improved MUSIC method［J］. Composite Structures， 2018， 203： 164-171.

[5] KUZNETSOV S V. Lamb waves in stratified and functionally graded plates： discrepancy， similarity， and convergence［J］. Waves in Random and Complex Media， 2021， 31（6）： 1540-1549.

[6] KUZNETSOV S V. Abnormal dispersion of Lamb waves in stratified media［J］. Zeitschrift Für An-gewandte Mathematik und Physik， 2019， 70（6）： 175.

[7] JIA H， ZHANG Z， LIU H， et al. An approach based on expectation-maximization algorithm for parameter estimation of Lamb wave signals［J］. Mechanical Systems and Signal Processing， 2019， 120： 341-355.

[8] BOASHASH B. Time frequency signal analysis and processing： a comprehensive reference［M］. Second edition. Amsterdam： Academic Press， 2016.

[9] KHAN N A， SANDSTEN M. Time-frequency image enhancement based on interference suppression in Wigner-Ville distribution［J］. Signal Processing， 2016， 127： 80-85.

[10] ZHANG P， TANG Z， LÜ F， et al. Numerical and experimental investigation of guided wave propagation in a multi-wire cable［J］. Applied Sciences， 2019， 9（5）： 1028.

[11] KAUFFMANN P， PLOIX M A， CHAIX J F， et al. Multi-modal leaky Lamb waves in two parallel and immersed plates： theoretical considerations， simulations， and measurements［J］. The Journal of the Acoustical Society of America， 2019， 145（2）： 1018-1030.

[12] TIWARI K A， RAISUTIS R， TUMSYS O， et al. Defect estimation in non-destructive testing of composites by ultrasonic guided waves and image processing［J］. Electronics， 2019， 8（3）： 315.

[13] ZENG L， HUANG L， CAO X， et al. Determination of Lamb wave phase velocity dispersion using time-frequency analysis［J］. Smart Materials and Structures， 2019， 28（11）： 115029.

[14] KUDELA P， RADZIENSKI M， FIBOREK P， et al. Elastic constants identification of woven fabric reinforced composites by using guided wave dispersion curves and genetic algorithm［J］. Composite Structures， 2020， 249： 112569.

[15] KUDELA P， KRAWCZUK M， OSTACHOWICZ W. Wave propagation modelling in 1D structures using spectral finite elements［J］. Journal of Sound and Vibration， 2007， 300（1/2）： 88-100.

[16] DRAUDVILIENE L， AIDER H AIT， TUMSYS O， et al. The Lamb waves phase velocity dispersion evaluation using an hybrid measurement technique［J］. Composite Structures， 2018， 184： 1156-1164.

基金资助

国家重点研发计划资助项目(2023YFC3010303)

国家自然科学基金资助项目(52272433)

江苏省市场监督管理局技术项目(KJ2024024)