引 言

近年来，信息科学发展迅速，提取真实世界的信息是它的重要任务之一，比如在控制工程、通信、军用电子和其他范畴的应用中，总是需要根据某种规律来确定微弱周期信号^[1]（即有用信号）是否存在。通常所提的微弱周期信号不仅指幅度很小的信号，也指存在于噪声背景下的信号。在实际工程应用中，有用信号的幅度和噪声信号相比要小得多，因此对有用信号进行检测是一个难题，为此人们进行了长期的研究。以往人们所研究的微弱周期信号检测方法是通过消除噪声^[2]或应用滤波，例如小波分析^[3,4]、相关积分等，但这些方法在信噪比较低时容易失效并且对有用信号有一定的损害。目前使用比较频繁的方法是Duffing振子混沌检测系统，不同于以往方法的是它对有用信号进行提取时并未对噪声进行抑制^[5,6,7]。当应用Duffing振子混沌检测系统对微弱周期信号进行检测时，因为Duffing振子对微弱周期信号特别敏感，使其易于检测有用信号；同时Duffing振子对噪声信号免疫，所以可以使用Duffing振子混沌检测系统来完成对混杂在噪声中的微弱有用信号的检测。

应用Duffing振子混沌系统检测有用信号的前提条件是确定系统相变的阈值，这直接决定该检测系统的有效性和准确性。文献[8]通过直接观察相轨迹图来确定相变阈值（即临界阈值），缺乏量化判断依据。文献[9]所给出的定量检测手段，是在忽略噪声影响的前提下实现的，而实际上噪声会对检测系统有一定的影响。文献[10]所提出的微弱周期信号检测方法的可检信噪比为-33.2 dB，但在更低的信噪比情况下不一定能检测出微弱周期信号或者对检测结果形成误判。文献[11]所构造的检测系统仿真模型，虽然该仿真模型能够检测出被覆盖在强噪声中的微弱周期信号，获得了相对满意的仿真效果，但这种仿真模型是针对频率 ω=1 rad/s进行检测的，没有讨论待检信号角频率 ω≠1 rad/s的情况。上述文献均是基于正向相变检测（从混沌状态到大尺度周期状态），但正向相变在由混沌状态跃迁为大尺度周期状态时极其容易受到过渡带影响。

根据以上所提到对微弱周期信号检测的不足，提出Duffing振子逆相变的检测方法，同时给出了改进后的检测系统仿真模型。改进后的仿真模型可以检测出任意已知频率的微弱周期信号，并利用最大Lyapunov特性指数得到检测系统产生相变的临界阈值，然后将改进后的Duffing振子混沌检测系统与高阶累积量相结合，通过求待检信号的高阶累积量对角切片能够降低噪声对系统的影响。通过理论分析和仿真实验，证明该方法可以检测出任意已知频率的微弱周期信号，其检测信噪比约-50.97 dB。

1 基于Lyapunov指数的微弱信号检测原理

1.1　Duffing振子检测原理

为了检测微弱周期信号，构造一个对微弱周期信号极其敏感的检测系统是非常有必要的，这是信号检测的前提条件。近几年常用的是著名的由Duffing方程构成的Duffing振子混沌检测系统。其形式如式(1)所示。

式中：k为检测系统的阻尼系数，αx(t)+βx³(t)为检测系统的非线性项，γcos(ωt) 为检测系统的内置周期策动力，γ 和 ω 分别为其幅度和频率。

为了检测微弱周期信号，应用的是对式(1)改进后的Duffing系统方程^[12]，表达式为

为了说明式(2)所对应Duffing系统的工作原理，在k取定值0.5，ω=1 rad/s的情况下，当内置周期策动力幅值 γ 产生变化时检测系统状态输出相图也会随之改变，依次为同宿轨道状态、分叉状态、混沌状态、大尺度周期状态。当 γ=0时，检测系统会出现两个焦点(+1,0)，(-1,0)和一个鞍点(0,0)，系统的相点最终会停在这两个焦点中的一个；随着周期策动力 γ 的继续加大，且当γ大于混沌阈值 γ_c 时，系统会出现混沌状态，系统的输出相图会在之后的一段时间内仍然保持这一状态，此时再增加周期策动力幅值 γ 使其大于另一阈值 γ_d 时，检测系统的状态就会发生改变即变成大尺度周期状态。系统的输出相图在这两个状态的差异非常明显，所以一般选取检测系统的输出相图从混沌状态向大尺度周期状态的转变来检测微弱信号。

Duffing振子混沌系统传统的检测方法是：调整内置周期策动力的幅值 γ=γ_d，这一时刻检测系统输出相图为临界状态；之后把淹没在噪声中的微弱周期信号作为内置周期策动力的补充加入到检测系统中，如果微弱周期信号是与内置周期策动力同频的且叠加后的幅值超过临界阈值 γ_d，那么检测系统的输出相图就会呈现大尺度周期状态，故能够通过相图是否发生变化来检测微弱周期信号的存在。

由于正向相变检测（混沌状态到大尺度周期状态）容易受到过渡带的影响，因此提出逆相变检测。先调整内置周期策动力幅值 γ 为临界阈值 γ_d，使检测系统的输出轨迹相图呈现临界状态。之后把待检信号传进到检测系统中，若待检信号中存在和内置周期策动力频率一样相位相反的微弱信号，那么相加后的幅值就会比 γ_d 小，检测系统的输出轨迹相图就会变为混沌状态，即产生了逆相变，从而可以检测出微弱有用信号。检测系统由临界状态转变为混沌状态的相图如图1所示。

图1 Duffing振子检测相图

Fig.1 Duffing oscillator detection phase diagram

1.2　基于Lyapunov指数计算临界阈值

Duffing振子混沌系统检测微弱周期信号的前提是确定临界阈值 γ_d，目前确定临界阈值的方法有直接观察相轨迹图法、Melnikov方法^[13]和Lyapunov特性指数方法。采用直接观察相图法的缺点是如果仿真时间过短会造成误判，很有可能会把仿真时刚开始呈现的形态错认为周期状态，如果再增加仿真时间才会呈现混沌状态，但如果加长仿真时间则会降低工作效率；基于Melnikov方法来计算临界阈值的准确性比较低并且误差也比较大；所以采用基于Lyapunov特性指数来准确计算临界阈值 γ_d，这是一种定量的判断标准，便于实际应用。利用Lyapunov指数判定检测系统从混沌状态转变到大尺度周期状态的原理是：当Lyapunov_max>0时，Duffing振子检测系统输出混沌状态的相轨迹图；当检测系统的Lyapunov_max>0跳变为Lyapunov_max<0时，检测系统的输出相轨迹图由混沌态转换为大尺度周期态，这一时刻所对应的内置周期策动力幅值即为所求的临界阈值。

Duffing方程形成的是一个二维非自治系统，采用变量代换令z=t将其改成相应的三维自治系统，故将式(2)改写为式(3)。其线性变分方程为

式中Jacobi矩阵J=0103x2-5x4-k-γsin(z)000。

然后再对F进行QR分解，即

式中：Q为一个正交矩阵，R为一个上三角矩阵。将式(5)代入式(4)得到

在式(6)的等式两端左乘Q^T，右乘R^-1可得

因为Duffing振子系统是一个连续系统，Geist等人结合Lyapunov特性指数定义^[14]证明出

由以上的推导过程，基于Lyapunov特性指数计算临界阈值 γ_d 的具体步骤如下：

（1） 给定时间步长；

（2） 选取内置周期策动力幅值 γ 的范围，修改其步长Δγ；

（3） 计算每次迭代中的J，Q；

（4） 对每次出现的F进行QR分解；

（5） 计算R·/Rii=(QTJQ)ii；

（6） 计算每次的LCEi=LCEi+(ln(Rii))'；

（7） 找到最大的LCE_max；

（8） 如果LCE_max>0，就返回步骤(2)，反之跳出循环，就可以得到临界阈值 γ_d。

在初始值设置完成并计算出临界阈值 γ_d 之后，将待检信号作为内置周期策动力的补充输入到Duffing振子检测系统得到输出时间序列，基于FFT算法来获取该模型输出时间序列的平均周期；再将模型的输出时间序列，通过C-C方法来计算延迟时间和嵌入维数并进行相空间重构；最后利用重构的相空间，再次计算检测系统的最大Lyapunov特性指数，因为运用的是逆相变过程，所以是通过判断检测系统的最大Lyapunov指数是否小于零。如果其值小于零，就调整周期策动力的幅值，并且重新计算调整后的Lyapunov特性指数，直到某一时刻所对应的最大Lyapunov特性指数刚好由小于零变成大于零，这样就计算出输入待检信号后新的周期策动力幅值 γ_d′。

2 改进方法

2.1　Duffing方程改进

上述过程讨论了微弱周期信号频率 ω=1 rad/s的情况，但所要检测的微弱周期信号频率不都为1 rad/s。所以有必要对方程进行改进，使得改进后的方程能够检测频率不是1 rad/s的微弱周期信号。对式(2)进行改进，暂不考虑待检信号。如

写成状态方程

根据式(10)所构造的Simulink仿真模型如图2所示。

图2 Simulink仿真模型

Fig.2 Simulink simulation model

在图2的仿真模型中，Fcn是函数发生器，Sine wave可以生成内置周期策动力，Add 1，Add 2是加法器；Integrator 1,Integrator 2是积分器；Gain 1,Gain 2是增益器；XY Graph是相轨迹显示器。下面讨论测量不同的待检频率 ω 值，临界阈值是否受到影响。

讨论Duffing方程x¨(t)+kx˙(t)-x3(t)+x5(t)=γcos(t)，现令t=ωτ，则

然后将式(11)代入Duffing方程中，可以得到系统方程为

可见此时式(12)是以时间尺度 τ 为自由变量的一个方程，由以上的推导过程可以看出Duffing振子混沌检测系统中的参数 ω 对临界阈值没有影响，若待检信号频率不等于1 rad/s，只要通过调整模型中的 ω 即可。比如当想测量 ω=10 rad/s的微弱周期信号时，只需要在模型中将放大器Gain 1的放大倍数设置为100，放大器Gain 2中放大倍数设置为5即可，表明给出的模型可以检测任意已知频率的微弱信号。

2.2　噪声对检测系统的影响

在检测过程中是基于混沌系统对微弱周期信号的敏感并且对噪声免疫的特性，因此在检测微弱信号时经常会忽视噪声对系统的影响，但经过大量的实验发现噪声会使检测系统的相轨迹图变得粗糙，如图3所示。并且当噪声功率比较大时即使待检信号里含有微弱周期信号也不能被检测出来。

图3 噪声对检测系统影响

Fig.3 Influence of noise on detection system

为了提高检测系统的性能，可以先对信号进行预处理，即通过求待检信号的高阶累积量^[15]降低噪声功率之后，再作为周期策动力摄动输入到检测系统。那么需要保证待检信号经过抑制高斯白噪声之后，信号的特征并未发生变化或者按照已知的可逆的规律改变，才能从所测得的信号参量中恢复出抑制高斯白噪声之前的原始有用的微弱周期信号（谐波信号），从而达到信号检测的目的。

待检信号y(n)=x(n)+s(n)，表达式y(n)中的s(n)为高斯白噪声，x(n)为要检测的微弱周期信号。具体情况分析如下：

（1） 高斯白噪声的二阶累积量不为零，但是其三阶和四阶累积量等于零；

（2） 微弱周期信号的二阶和四阶累积量不为零，然而其三阶累积量等于零。

通过上面的分析可知不能利用二阶累积量和三阶累积量。只可以利用四阶累积量，四阶累积量对噪声的抑制非常明显，但仅仅对噪声有抑制作用是远远不够的，还必须保证对信号无根本上的影响。下面探讨四阶累积量对微弱周期信号和高斯白噪声的影响。

取x=[x₁, x₂, …, x_k]^T为一个随机矢量，并且s=[s₁, s₂, …, s_k]^T，因此随机矢量x的第一特征函数^[16]（即矩生成函数）为

由式(13)可得到x的累积量生成函数表达式为

x的(v₁, v₂, … , v_k)阶矩和累积量分别可由第一特征函数和累积量生成函数在s1v1,s2v2,…,skvk处进行泰勒级数处展开得到，即

式中满足v= v₁+v₂+…+v_k。

在v₁=v₂=…=v_k=1的条件下，随机矢量x的矩和累积量分别等于

随机过程的高阶矩和高阶累积量便是基于式(17)进行定义的。

定义：若{x(n)}为一个k阶随机过程，则其k阶矩和k阶累积量分别为

根据式(18)可知，随机过程的k阶矩和k阶累积量即为x₁=x(n), x₂=x(n+τ₁), …, x_k=x(n+τ_k-1)之后的随机矢量[x(n), x(n+τ₁), …, x(n+τ_k-1)]^T在v₁=v₂=…=v_k=1条件下的矩和累积量。

依据该定义，可以得到随机过程{x(n)}的k阶矩和k阶累积量，如式(19—22)所示。

设x(n)=∑m=1pamcos(ωmn+ϕm)为一谐波过程，根据四阶累积量的计算公式可得其计算结果为

当不存在频率耦合时，只有前3项，如式(24)所示。

但是由式(24)可以看出此时四阶累积量的计算仍然非常复杂。因此为了方便计算同时不能损失有用信息，经常通过计算其对角切片。四阶累积量对角切片^[17]实际上是让c_4x(τ_1,τ_2,τ₃)表达式中τ₁=τ₂=τ₃=τ 所得的四阶累积量，简化表示为c_4x(τ)，将τ₁=τ₂=τ₃=τ 代入式(24)得到

s(n)为一高斯过程，根据表达式(19—22)并结合高斯过程的特点，能够证明其矩和累积量为

从式(27)可以看出，高斯过程的四阶累积量为零。

由高阶累积量的性质可知

由以上分析可以看出，四阶累积量对角切片对微弱周期信号的影响只是幅度的影响，并且知道幅度变化的规律，但经过计算四阶累积量可以使高斯白噪声得到抑制，所以可以利用求待检信号的四阶累积量对角切片来消除高斯白噪声对检测系统的影响，从而可以提高检测信噪比。然后再通过谐波信号幅度变化式(28)来计算实际的谐波信号幅度。

基于以上分析本文检测系统的设计方案如图4所示。

图4 检测系统总体方案

Fig.4 Overall plan of detection system

3 仿真实验

根据上述讨论，首先取k=0.5，不加入待检信号，先计算临界阈值。然后按照求临界阈值的方法，设置时间步长h=0.01，初值设置成（0,0,0），初始周期策动力幅值取0.720，迭代步长Δγ=0.001，以上参数全部设置完成后开始运行程序，结果见图5。

图5 周期策动力幅值γ和最大Lyapunov指数的关系图

Fig.5 Relationship between the amplitude of periodic power γ and the maximum Lyapunov index

运行程序求得 γ_d =0.732 0。将基于Lyapunov指数求得的临界阈值所用时间和另两种方法进行比较，结果如表1所示。

由表1可以看出，基于Lyapunov指数求临界阈值不仅精确度比较高，所用时间也最短，再一次说明了Lyapunov的高效性和准确性。

在求出临界阈值 γ_d 之后将待检信号y(t)=x(t)+s(t)作为周期策动力的补充加入到检测系统。其中，x(t)有3种情况，分别为x₁(t)=0.01cos(2πt)，x₂(t)=0.02cos(2πt)，x₃(t)=0.05cos(2πt)。s(t)是均值为0、方差为2.5的高斯白噪声（Noise-92^[18]），3组实验均基于MATLAB平台进行仿真。具体过程如下：

（1）先对待检信号进行预处理即通过计算待检信号的四阶累积量对角切片，此处理方法可以很好地滤除高斯白噪声。

（2）在仿真模型中，先不输入经过处理的待检信号，只加入内置周期策动力，设置其幅值 γ=γ_d =0.732 0，此时检测系统的输出相图是临界状态。

（3）将经过处理后的待检信号作为周期策动力的补充加入到逆相变检测系统中，此时检测系统输出相图就会变成混沌状态。经过对检测系统的输出时间序列进行相空间重构并再次计算检测系统的最大Lyapunov指数，求出当最大Lyapunov指数由小于零到大于零符号转变的那一刻所对应的策动力幅值 γ_d′，微弱信号为x₁(t)时所对应的 γ_1d′=0.732 000 003 72，通过计算可得输入信号幅值A₁为3.72×10^-9；微弱信号为x₂(t)时所对应的 γ_2d′=0.732 000 053 117，通过计算可得输入信号幅值A₂为5.311 7×10^-8；微弱信号为x₃(t)时所对应的γ_3d′=0.732 002 517 1，通过计算可得输入信号幅值A₃为2.517 1×10^-6。

（4）由输入信号幅值为A=0.375a⁴，可以分别求出待检信号幅值a₁，a₂，a₃。和用基于正相变未进行高阶累积量预处理的直接观察法得到的幅度估计值进行对比，结果如表2所示。由表2可看出本文方法和直接观察法相比对待检信号幅值的估计结果相对误差更小。

（5）计算最低信噪比SNR=10log0.012/22.52=-50.97dB。

要检测信号的原始波形和含有噪声的信号波形以及抑制高斯白噪声之后的信号波形如图6所示。

图6 信号波形

Fig.6 Waveform of the signal

由图6可以看出，对待检信号进行求四阶累积量一维切片可以消除高斯白噪声的影响，使噪声的功率降低，从而提高检测的信噪比。这里有一个问题，在实际工程应用中所检测的微弱周期信号幅值根本达不到10^-2，量级一般为10^-7或10^-8，此时的检测系统不可以直接用来检测这样的信号，由前面的推导过程，10^-2的信号经过计算四阶累积量之后，就变为10^-8。但是可以先将处理后的待检信号经过一个信号放大器，这样幅度就得到了放大，之后再加入到检测系统进行检测。

然后将没有经过处理的待检信号也加入到逆相变的检测系统中，发现检测系统的状态没有发生变化仍然是临界状态。将经过处理后的待检信号输入到逆相变系统，此时检测系统的输出相图由临界状态变成混沌状态，结果如图7所示。由图7可以看出，经过处理后的待检信号可以被有效地检测。在这之后又进行多次实验再次证明了本文所提方法的有效性。将待检信号没有被处理和经过处理的检测结果进行比较，结果如表3所示。由表3能够看出和没有经过计算高阶累积量的传统Duffing振子检测系统相比，本文所提出的方法可以检测更低信噪比条件下的微弱信号。

图7 采用高阶累积量对信号处理前后检测系统输出相图

Fig.7 Output phase diagram before and after signal processing by high-order cumulants

4 结束语

本文是将常规的微弱信号检测方法（高阶累积量）和Duffing振子方法相结合并且基于逆相变对微弱信号进行检测。首先选取的是改进后的仿真模型，该仿真模型可以检测任意已知待检信号频率，在推导出仿真模型中的相关参数和待检信号频率之间的关系后设置该待检信号所对应的参数值。然后基于最大Lyapunov指数得到检测系统的临界阈值，该方法提高了求解临界阈值的速度和精确性；确定临界阈值后将待检信号输入到检测系统中进行检测，此时的待检信号是经过处理的已经降低了噪声对检测系统的影响。最后的仿真结果证明所提方法可以有效检测出被强噪声淹没的微弱信号，为检测任意未知的微弱信号提供出新的思路。

参考文献