Duffing振子变尺度未知多频微弱信号提取方法研究

黄继尧，陈长兴，凌云飞，蔺向阳

(空军工程大学基础部，陕西西安 710043)

0 引言

Duffing振子是一种用于信号处理的非线性混沌振子模型，它在混沌临界状态下对噪声具有免疫性，但对初始条件极为敏感[1]。可以利用Duffing振子的这些特性来检测强噪声背景下的微弱信号。自20世纪后期引入混沌振子系统以检测微弱信号[2-4]以来，诸多学者在此领域进行了深入的研究。谢涛等证明了混沌振子检测微弱信号的可靠性，为混沌检测微弱信号提供了理论论据[5]。

李月等分析了不同类型噪声对混沌系统的影响，证明了混沌系统对噪声的免疫性[6]；秦卫阳等分析了Duffing方程的派生系统在同步过程中能够实现检测2种信号中的微小差别[7]；刘海波等提出一种利用单一振子实现以Duffing振子大周期态为基础的测量频率的方法[8]。常见的混沌状态判定方法有Lyapunov指数法、庞加莱截面法、分维数计算法、Kolmogorov熵法、Melnikov法等[9]。与传统的线性信号检测方法相比，基于Duffing振子的微弱信号检测方法可以实现非常低的信噪比(SNR)，信噪比门限最低可达到-100 dB，并可得到传统检测方法无法得到的微弱待测信号，这为低信噪比的微弱信号检测提供了一种新的方法[10]。

然而目前大多数基于Duffing振子的微弱信号检测研究都是针对单频谐波信号的检测识别[11-14]，而很少有研究关注微弱多频信号的检测。因此对基于Duffing振子的微弱多频信号检测进行深入研究是十分有必要的。本文提出了一种结合信号频率截取预处理和Duffing振子检测微弱多频信号的方法。

1 Duffing振子微弱信号检测原理

Holmes型Duffing方程[15]一般形式为

(1)

式中：k为阻尼比；F为周期性策动力幅值；ω为角频率；Fcosωt为周期性策动力。

当k保持不变时，随着F的增大，存在一个临界策动力振幅Fd。当F小于Fd时，系统处于混沌状态；当F大于Fd时，系统处于大周期态。在Duffing振子方程式(1)中，临界值Fd≈0.825。混沌和大周期态的变化及策动力临界值如图1所示。

(a)混沌态

(b)周期态

(c)临界点图1 系统历经的状态变化及分岔临界图

当Duffing系统处于从混沌状态转变为大周期态的临界状态时，就统计特性而言，它不受高斯白噪声的影响。由于对参数 的敏感性和对噪声的免疫性，式(1)可以进一步表示为

(2)

式中：s(t)为待测信号；n(t)为噪声；s(t)+n(t)为带有噪声的待测信号。

由式(2)可知，当式(1)中的F略小于临界幅值Fd时，s(t)+n(t)输入Duffing振子系统，n(t)作为扰动项，则可以使用式(2)检测待测信号。sn(t)=s(t)+n(t)是含噪待测信号。当s(t)包含角频率为ω的信号分量时，系统将从混沌状态转换为大周期态，这可以从系统响应的相平面轨迹获得，然后便可以获得待测信号的频率。

傅里叶变换是一种分析信号的方法，可以把具有周期特征的非正弦信号分解为多个叠加的正弦信号，然后通过对时间连续信号的等间隔采样，得到采样频率fs，并把采样值转换为数字序列。连续时间序列x(t)与其对应的连续非周期频谱函数X(f)之间的转换关系为：

(3)

当设定一个待测信号为hcos(ωt+φ)，其采样频率为fs。ω为确定的角频率，φ是初始相位。假设采样信号向坐标轴左侧移动dt=1/fs。根据公式T/dt=2π/dθ，就可以得到dθ=dt·2π/T=ω/fs，T为信号的周期，dθ为对应相位的增量。因此，将采样信号向左挪动时间j等同于将待测信号由hcos(ωt+φ)变成hcos(ωt+φ+jdθ)。

在不考虑噪声，即n(t)=0的情况下，设定待测信号尺度变换后(尺度变换系数R=ω)与策动力有相同的角频率，则检测待测信号的Duffing方程为

(4)

改变式(6)右侧项的形式

(5)

从式(4)中，可以推出只有当(φ+jdθ)=2mπ(m是整数)时，与待测信号和策动力相加的待测信号的振幅是最大的。

文献[16]提出了一种基于Duffing振子的尺度变换微弱信号检测模型，有效地解决了这些问题。检测模型如式(6)所示：

(6)

这4个方程的输入信号初始相位的检测范围能够覆盖整段频率(-π，π)，避免了单一方程只能检测某一段初始相位的限制。

式(6)解决了基于Duffing振子微弱信号检测的一些关键问题，然而目前仅对含噪声的单频谐波信号进行分析，而不涉及实际工程应用中常见的微弱多频信号的检测。因此本文将研究基于Duffing振子的微弱多频信号检测。

2 基于频率截获Duffing振子的微弱多频信号检测

2.1 信号频率截取方法

为了将式(6)中的检测模型应用于微弱多频信号，提出了一种信号频率截取的预处理方法。其原理如下:根据Shannon采样定理，试验信号sn(t)的采样频率为fs，则其截止频率为fs/2，为了便于后续分析，首先在频域中截取待测信号。截取频率的下限为fmin，上限为fmax，其满足0

首先，利用快速傅里叶变换(FFT)将试验信号sn(t)变换为频域信号sn(f);其次，信号sn′(f)是由sn(f)通过频率截取得到的，即离散谱数据的频率范围为(fmin,fmax);最后，将信号sn′(f)通过快速傅里叶逆变换(IFFT)转变为时域信号sn′(t)。

上述频率截获过程可以表示为:

(7)

(8)

(9)

频率截取过程后，sn′(t)只包含满足fmin

2.2 待测信号提取模型

将信号频率截获方法与式(6)检测模型相结合，针对微弱多频信号，提出了一种新的待测信号提取模型。其过程具体如下:

(1)对原始待测信号进行加噪处理，得到被噪声污染的待测信号sn(t)，采样频率为fs， 由强度为D的背景噪声和微弱待测信号s1(t),s2(t),…,sn(t)构成，微弱待测信号的角频率和振幅分别为ω1,ω2,…,ωn和h1,h2，…,hn。

(2)确定原始待测信号的频率测试范围[ωmin,ωmax]。然后设定ω=ωmin，接着设置频率增量dω，令ωi+1=ωi+dω。

(3)然后对待测信号sn(t)进行傅里叶变换后得到信号sn(f)再进行频段为[ωmin,ωmax]的频率截取预处理。新的信号在预处理后为sn′(f) ;

(4)设定尺度系数Ri=ωi(即ω的值每增长一次，R的值都会相应的增长并与之相等)，然后对信号进行频率/时间尺度变换，之后得到在时间尺度t′下新的待测信号sn′(t)。变换后的采样频率为fsr=fs/Ri；

(5)频率截取预处理和尺度变换后，将新获得的信号sn′(t)输入到式(6)中。然后利用四阶Runge-Kutta算法求解各方程，绘制各输出响应的相平面轨迹图。只要任何一个相平面轨迹处于大周期态，则意味着待测信号sn(t)中的频率分量ωk被检测到，然后重复迭代过程，直到频率满足ω≥ωmax为止，过程如图2所示。

图2 信号检测框图

从而检测到该频率范围内的所有信号频率分量，实现对未知多频率信号的检测。如果dω太小，检测未知信号的计算量将大幅增加。选择角频率为ω=20 rad/s的待测信号进行测试，设定Δω为0.2 rad/s。选择不同的Δω的倍数与计算时间的关系如图3所示，选择dω为2Δω有较快的运算速度，因此选择dω为2Δω。

图3 不同的Δω倍数与计算时间的关系图

所以当检测过程完成后，可以覆盖完整的频段[ωmin,ωmax]。迭代过程如图4所示：

图4 未知多频信号检测的迭代过程

3 仿真与分析

设定待测信号的方程为

(10)

式中：ω1,ω2,ω3和h1,h2,h3分别为sn(t)中包含的3个频率分量的角频率和振幅；g(t)为均值为0，方差为1的高斯白噪声。

设定待测信号的初始相位为0，待测信号的采样频率fs=100 Hz，h1=h2=h3=0.1，ω1、ω2、ω3分别为20，30，50 rad/s，噪声强度D=1，Δω为0.2 rad/s。设定[ωmin,ωmax]为(0,60)。绘制待测信号的时域图和频域图如图5所示，波形显示3 000个采样点。

由于在实际工程问题中噪声一般是不可忽略的，因此在时域图和频域图中，可以轻易发现由于噪声的污染，很难找到待测信号的频率的存在。而相图的检测也同样容易受到噪声的干扰，如图6所示。

将ω=20 rad/s的待测信号送入式(6)中的第一个方程，可以发现图6(a)所显示的相图在噪声的影响下，较难分辨是出于混沌状态还是大周期态。而图6(b)所显示的则是在经过尺度变换后消除了噪声的影响，相图是处于大周期态的。在下文中，都将以经过尺度变换后的相图阐述。

采用本文提出的优化模型进行检测，在ω=10 rad/s,20 rad/s,30 rad/s,40 rad/s,50 rad/s,60 rad/s，6个点处进行取值，如图7所示。

(a)时域图

(a)未经过尺度变换的相平面轨迹

(b)经过尺度变换的相平面轨迹 图6 ω=20 rad/s的待测信号的相平面轨迹

通过实验可以得知，图7(b)、(c)、(e)所显示的是分别为ω=20 rad/s、ω=30 rad/s和ω=50 rad/s的频率分量所对应的相图，相图是明显处于大周期态的。而图7(a)、(d)、(f)所显示的分别是在ω=10 rad/s、ω=40 rad/s和ω=60 rad/s的频点所对应的相图，相图是处于混沌状态的，符合实验设定的预期条件。

实验发现：

(1)未使用尺度变换，就很难将噪声污染的相图分辨出来，存在对相图的误判可能，尺度变换可以有效地消除噪声对相图的干扰，尺度系数R的存在有利于待测信号的提取。

(2)使用本方法成功检测到了ω=20 rad/s,ω=30 rad/s，ω=50 rad/s的频率分量。而在ω=60 rad/s这两个频点则相图未发生状态转移，证明不存在该频率的分量，具有普遍性和实用性。

(a)ω=10 rad/s的相平面轨迹

(b)ω=20 rad/s的相平面轨迹

(c)ω=30 rad/s的相平面轨迹

(d)ω=40 rad/s的相平面轨迹

(e)ω=50 rad/s的相平面轨迹

(f)ω=60 rad/s的相平面轨迹

(3)本方法实现的微弱信号检测最低信噪比门限可以达到-24.2 dB，低于文献[16]中提出的原检测模型最低信噪比门限的-13.01 dB。

(4)Δω的选取如果过大，则会出现漏检测现象；Δω的选取如果过小，则检测效率就较低，因此需要找一个平衡点。一般采用比待测信号频率低3个数量级的Δω可以实现限定频率范围内全部频率的检测。本文所提出的方法在实际实验中采用的Δω是0.2 rad/s，检测到了3个待测信号的频率。

(5)本文所提出的方法与Duffing振子并行处理方式的区别在于，本文所提出的方法是对频域加带通滤波器的基础上进行而且仅用一个Duffing振子便可检测到限定频域内全部的待测信号，而Duffing振子并行处理方式一般采用阵列的方式进行检测，即需要较多振子。本文提出的方法节省了振子的数量，对运算条件需求较低，在实际实验处理中也可以达到较为理想的效果。

4 结束语

本文研究了一种未知多频微弱信号检测方法，该方法针对原检测模型只能同时检测某一频率的待测信号，克服了噪声诱导对检测结果的不利影响，提出了一种可同时检测多个未知频率微弱信号的Duffing振子优化模型，可以实现提取未知多频微弱信号各信号分量的频率参数，具有实用性和普遍性。仿真结果表明：该方法为微弱多频信号处理提供了一种解决方案，拓宽了Duffing振子在微弱信号检测领域的应用。