基于Duffing-Holmes系统的材料非线性的量化检测方法

柏 林， 唐 滔， 刘小峰， 韦代平

(重庆大学 机械传动国家重点实验室，重庆 400044)

针对材料非线性检测，有研究表明，非线性Lamb波的传播会由于材料自身的非线性以及其结构微裂纹的局部接触而激发非线性Lamb波的非线性变化[1-2]，因此可利用高次谐波信号对材料性能进行实时程度评价。大量试验证明，通过测量Lamb波非线性效应可以有效反映材料结构的微观变化，以此实现早期损伤或性能退化的检测[3-4]。但目前的研究成果中，针对非线性Lamb波的信号的特征增强以及抗噪能力的研究很少。并且在试验中，由于Lamb波自身的扩散角以及频散特性等，致使接收信号极易受噪声污染从而失真。为了提高超声波检测技术的精度，采用非线性信号传统的降噪方式，如Kalman滤波[5]，LMS自适应滤波[6]，以及小波阈值降噪[7]等，都已经难以满足对微弱信号的检测要求。鉴于非线性Lamb波对于材料特性检测的优势性[8]但易受外界干扰，实现非线性Lamb波抗噪量化分析对对板构件的早期性能退化识别具有重要意义。

由于杜芬混沌振子对外加扰动具有初值敏感性及噪声免疫性，在微弱信号的检测中得到了广泛的运用。但目前研究都主要停留在混沌系统的对微弱信号的定性分析上，而针对非线性Lamb波的抗噪定量分析并未涉及。运用杜芬振子对微弱非线性Lamb的定量研究也尚属于发展阶段。论文利用在不同幅值激励下杜芬系统混沌程度的差异性，并通过Lyapunov指数对混沌特性进行量化，建立激励幅值与系统输出的线性模型用于实现对噪声干扰下的微弱信号量化分析。

1 非线性Lamb波传播原理

当Lamb波经过发生换能器输入材料中，材料在超声波作用下内部空间受到扰动，在内部微小裂纹以及材料的非完全弹性的共同作用下，激发出高阶谐波信号。对于纵向传递的Lamb波，即u0=A1sin(ωt)，联立一维波动方程

(1)

式中：x为纵向传播距离;c为纵向传播波速;β′即为Lamb波非线性系数。当只考虑波在传递过程中存在的频散现象时，式(1)有唯一解为

u=A1sin(kx-ωt)+0.125β′(kA1)2xsin 2(kx-ωt)

(2)

式中:k为基波波数;ω为基波圆频率;A1为基波幅值。可以看出对于二次谐波幅值A2，在波的传递过程中材料的非线性具有累积效应，并与非线性系数β′存在线性关系

β′=(8A2)/(kA1)2x

(3)

在检测距离x恒定的情况下，重构非线性系数β

β=A2/(A1)2

(4)

在实际非线性Lamb波检测过程中，二次谐波幅值与基波相比非常小，常常受到噪声信号的干扰。因此，β指数抗噪能力不佳。

2 Duffing-Holmes系统与Lyapunov指数谱

为了实现杜芬系统在不改变系统参数的情况下，对任意频率谐波检测的普适性。本文采用改良后的Duffing-Holmes系统，利用待测频率ω带系统的时间尺度进行放缩处理，其三维自治状态方程为

(5)

式中:k为阻尼比;x3-x5为系统非线性恢复力;Fcos(ωt)为系统内周期策动力;ω即为待测频率;F为内策动力幅值。(s(t)=Acos(ωt+φ)+n(t))外部系统扰动，其中:Acos(ωt)为待检信号;n(t)为待检信号中的噪声成分。

Duffing-Holmes系统相轨迹具有如下特性：当系统阻尼k固定，内策动力幅值F从0逐渐增加临界点F0时，相轨迹状态由周期内轨运动过渡到倍周期分叉最后进入混沌状态。一旦F超过F0，系统动力学特性发生突变，整体进入大周期状态，即发生正相变。针对式(5)，待测信号s(t)中存在与内策动力同频谐波信号，虽然同频谐波的合成并不影响频率成分，但改变了系统的策动力幅值与相位，触发临界系统正相变，即系统由混沌状态向大周期状态跳变。另一方面，当输入信号为噪声信号或不含待测成分时，将不会对系统的周期策动力F造成稳态影响，Duffing系统基本免疫。

在混沌理论中，Lyapunov指数可以在一定范围内定量描述非线性系统的动力学行为特征，并能在一定时间范围内对不便于观察的微弱信号进行幅值估计。其中，Lyapunov指数沿某一方向的正负与大小，分别表征系统在该方向上的相邻轨迹的平均发散与收敛程度。当系统的最大Lyapunov指数为正时，Duffing系统在相空间中相邻轨迹成指数发散，系统呈现混沌状态。

对于式(5)的三维Duffing自治系统方程的Lyapunov指数，如式(2)，此时系统的Jacobian矩阵为

(6)

式中:s′(t)为系统外部扰动关于时间的一阶导数;ω为待测频率。在系统初始阶段，令初始值x0=y0=t0=0，并将式(5)表示为3维连续动力学系统

(7)

此时，令Φ(t)为式(5)的基础解空间，则该系统的线性变分方程可表示为

(8)

式中：J(t)为Duffing系统t时刻的3维Jacobian矩阵;I3为3×3的单位阵。利用4阶Runge-Kutta法可分别计算变分方程在不同时刻的Jacobian矩阵。并对基础解空间Y(t)进行QR分解，即Φ(t)=Q(t)R(t)。在文献[9]中指出，系统的变分方程式(7)对应的Lyapunov指数满足

(9)

式中：Rii(t)为上三角矩阵R(t)的正对角项。

3 基于谐波逼近的信号幅值估计

临界Duffing-Holmes系统具有参数敏感性以及噪声免疫性[10]。对于混沌振子系统，系统扰动s(t)可分为由系统频率决定的慢变周期扰动，以及由高频信号决定的快变周期扰动，已有研究表明慢变周期扰动能使不稳定的相轨迹进入稳定周期[11]，而快变信号仅能驱动系统的局部振荡。对于线性放缩后的杜芬系统，Acos(ωt+φ)为慢变周期成分，噪声成分n(t)相对于系统频率ω仍属于快变周期成分。鉴于Duffing系统动力学特征——Lyapunov指数——与待测慢变周期成分幅值A具有一定的相关性，利用此间关联对噪声下的信号成分进行幅值定量分析，具体步骤如下：

步骤1 根据Lamb波的采样频率fs及待测频率ω，设置杜芬振子系统的缩放系数ω及分析步长dt=1/fs；

步骤2 根据系统分岔图，确定混沌临界阈值F0；

步骤3 由于接收到的Lamb波属于瞬态信号，而杜芬系统一般用于稳定谐波检测，因此需要对采集到的Lamb波进行周期延拓处理；

步骤4 为降低一次谐波对二次谐波检测的影响，采用高通滤波方法滤除一次谐波，通过幅频特性确定滤波器对二次谐波幅值影响近似忽略；

步骤5 将滤波后信号进行γ=0.4-1.6的等比放缩，并作为系统扰动信号输入混沌系统。利用最小二乘拟合建立放缩系数γ与Lyapunov指数之间的线性方程模型

LE=a×γ+b

(10)

式中：γ即为参考信号幅值A0与待测幅值Ax之比。

步骤6 为减小噪声带来的局部振动的影响，通过小波软阈值法提取噪声信号，并与相同的采样频率下的参考谐波信号融合，其中，仿真谐波的幅值已知，即为式(10)中的参考幅值A0；

步骤7 测试参考信号在经过相同演化时间后Lyapunov指数，并结合式(10)，求解待测幅值Ax。

4 试验测试

本章通过试验验证前文提出的方法的可行性。在谐波幅值定量分析试验中，通过收集非线性Lamb波在非线性材料中传播一定距离后的信号，测量其基波幅值与二次谐波幅值，以此获得近似无噪的初始信号。

图1 试验设备Fig.1 Experimental equipments

试验使用的非线性超声测试系统为RAM-500 SNAP (RITEC Inc., Warwick, RI), 如图1所示，具有检测超声衰减和波速的能力。相较于传统线性检测法，RAM-500具有更好的灵敏度。图1中为本次试验对象1.5×620×2 500 mm的无损铝板，铝板平行放置在海绵块上，无固定约束。在基准线两端分别布置压电式发射换能器与PZT接收换能器。本次试验距离700 mm，激发信号200 kHz，试验周期数30，采样频率为10 MHz。由于换能器与铝板之间存在空隙，试验通过甘油作为超声耦合剂排除间隙空气，以便Lamb波的传递。

调整Duffing-Holmes系统内策动力幅值F进入临界状态，并设置式(5)中的系统参数k=0.5[12]，此时系统的临界阈值F0=0.725[13]。分析试验采集信号，提取出采集到的Lamb波包作为初始信号，并对原始信号做零相位带通滤波后得到如图2，分析二次谐波与基波的时域关系。

图2 谐波时域对比图Fig.2 Comparisons of original wavepacket and extracted harmonic

根据图2，在有限时域采样长度N内，由于频散特性以及波在传递过程中的扩散，二次谐波在时域空间中严重发散。根据信号的单边幅值谱定义

(11)

对于瞬态信号，其对应频率的幅值由信号采样长度N以及瞬态信号长度N1共同影响。但相比与基波长度N，二次谐波长度N1难以确定。为防止由于截断时长带来的幅值谱差异，对待测信号统一补零到相同时窗长度后计算幅值谱。取接收到的Lamb波包信号进行分析，如图3(a)所示；为模拟试验干扰，在接收信号中加入标准差为1×10-3的白噪声，即图3(b)。

(a)原始信号

(b)加噪信号图3 原始信号与加噪信号对比图Fig.3 Comparison of the received signal and noised signal

从图3(a)可知，接收到的Lamb波包中二次谐波非常微弱，与一次谐波幅值相差2个数量级。在试验中，即使存在微弱干扰，也会使得谐波幅值剧烈波动。在图3(b)中二次谐波幅值依然清晰可见，但根据表1可知，此时二次谐波幅值明显失真，传统非线性指数β在干扰情况下已失效。

采用本文方法首先对Lamb波包进行整周期延拓滤波后，并乘以缩放系数γ=1后作为系统扰动参数s(t)输入系统。通过4～5阶变步长Runge-Kutta法求解，得出系统的相轨迹图与Lyapunov指数谱如图4(a)所示。

表1 不同算法条件下的β指数对比Tab.1 Comparison of calculated β with different algorithms

在图4(b)中，当系统无扰动输入时，混沌系统存在奇怪吸引子。对于三维Duffing-Holmes自治系统，LEx，LEy和LEt分别表示系统沿x轴，y轴，时间轴t方向的Lyapunov指数。

(a)γ=1时杜芬系统输出

(b) 无输入时杜芬系统输出图4 系统输出响应对比Fig.4 Comparison of chaotic system outputs

此时沿某一方向相轨迹必发散，即图中沿x轴方向LEx>0。在图4(a)中，由于系统扰动项中存在稳定待测频率成分，混沌系统出现普通吸引子——极限环，在Lyapunov指数谱中，有且仅存沿时间轴方向动力学指数为0，其余指数皆小于0。经过试验分析，随着待测成分幅值Ax增大，杜芬系统整体趋向稳定速度加快，在相同演化时间内与沿x轴方向Lyapunov指数呈现负相关趋势。本文利用LEx分析系统动力学特性。

调整扰动信号的缩放系数γ∈[0.4,1.6]，在经历相同的演化时间后，记录对应的沿x轴方向的Lyapunov指数LEx，如表2所示。

表2 放缩系数γ与Lyapunov指数对应关系Tab.2 The relation between the coefficient of contraction γ and the Lyapunov exponent

图5中，利用最小二乘法对得到的离散数据进行一次线性拟合，得到放缩系数γ与LE之间的线性模型(拟合确定度R-square: 0.968 1)

LE=-0.004 52γ-0.078 1

(12)

图5 γ与LE的一阶线性拟合Fig.5 Linear fitting between γ and LE

在整个时域空间内，混沌系统的整体趋势主要受到待测频率成分幅值影响，而噪声成分仅造成局部的波动。因此，对于扰动信号的放缩，可以等价视为不同幅值的无噪谐波信号直接作用于混沌系统，此时系统状态将表现出于待测信号缩放系数之间的直接相关性，即式(12)，系统状态差异通过Lyapunov指数进行量化。由式(10)可知，当系统线性模型确定，为定量分析噪声下的准确二次谐波幅值Ax，先引入参考信号(波数30，频率0.4 MHz，采样频率10 MHz的正弦谐波信号)，此时仿真信号并没有考虑波的扩散，将参考信号通过补零的方式延拓至与原始信号相同的采样长度，分析此时幅值谱，如图6所示。

图6 同频率基准信号Fig.6 Reference signal of the same frequency

图7中，在相同的采样长度下，参考信号幅值为5.508×10-5V。考虑噪声带来的局部振动，对原始信号进行小波软阈值信噪分离提取噪声成分，如图7所示。

图7 提取噪声信号Fig.7 Extracted noise

将参考信号与噪声成分结合延拓后作为系统扰动输入杜芬系统，经过相同的演化时间后，沿X轴方向动态指数LEx=-0.083 8，如图8所示。联合线性方程式(10)与式(12)，代入参考幅值A0=5.508×10-5V。代入求解参考信号与真实信号之间的放缩系数γ=A0/Ax=1.26，即预测二次谐波幅值为Ax=4.371 4×10-5V。

图8 参考信号的Lyapunov指数Fig.8 Lyapunov exponent of the reference signal

通过表1，当加入微量的噪声，此时基波幅值基本维持不变，但二次谐波参数剧烈波动。通过小波阈值法可以看出，当噪声信号与待测信号出现在同一频率上时，小波并不能实现很好的分离。Kalman滤波器在信噪分离阶段，分离出的信号并不是最佳估计，信号失真严重。LMS自适应滤波收敛速度较慢，信噪分离效果不佳。而此时混沌预测达到了较好的估计效果。

5 结 论

本文利用Duffing-Holmes振子与Lyapunov指数谱在噪声干扰下，构建待测信号幅值与系统动力学特性之间的线性模型并对材料非线性进行量化分析。从试验结果可以看出，该方法利用杜芬系统对噪声免疫的特性，实现在白噪声干扰下的微弱信号幅值估计。从而避免了对带噪信号抑噪而引入的失真，体现了该方法检测微弱信号幅值的优越性。