内置频率对Duffing振子微弱二进制相移键控信号盲检测影响*

吴彦华，马庆力

(国防科技大学 电子对抗学院， 安徽 合肥 230037)

二进制相移键控(Binary Phase Shift Keying， BPSK)具有较高的频带利用率、较强的抗噪声干扰能力以及具备展宽信号的带宽的优势，难以被探测。

对于第三方无线电监测来说，对微弱BPSK信号的盲检测是一个重要课题。当前，针对PSK信号的盲检测已有很多方法：文献[1]采用高阶累积量和循环谱检测方法，文献[2]采用联合频谱估计与循环矩检测方法，文献[3]采用最大自然准则估计方法，文献[4]采用基于星座图检测方法，文献[5]采用随机共振PSK信号参数估计方法，以及早期采用各种时频分析检测方法等。这些方法中大部分需要在信噪比大于0 dB时才能取得良好效果，对于信噪比在-10 dB以下的微弱PSK信号检测没有太多帮助。

Duffing振子对与系统策动力同频率的小信号具有敏感性，而对高斯噪声信号具有极强免疫力。利用Duffing振子在混沌态到周期态的分岔行为，可以判断强噪声中是否存在微弱周期信号。

当前，采用Duffing振子检测微弱信号理论和方法的研究成果主要集中在微弱周期信号的检测方面，在Duffing振子检测PSK(包括BPSK)信号的研究方面，文献资料比较少。文献[6]建立了识别PSK信号识别模型，但它们只针对PSK信号频率与Duffing系统内置信号频率同频的情况，不适用于PSK信号盲检测。文献[7]建立了基于Duffing振子的PSK信号调制识别算法，但同样需要知道载频频率的先验知识，并且信噪比要求较高(-5 dB以上)。

本文推导了Duffing系统内置频率对检测灵敏度和过渡带时长影响的关系表达式，建立了采用Duffing振子和S变换对任意频率微弱BPSK信号进行盲检测的模型。

1 问题描述

1.1 Duffing振子微弱周期信号检测

考虑Holmes型Duffing方程[8]：

(1)

式中：k为阻尼率；-x+x3为非线性恢复力；Fcos(ωt)为系统内置信号即系统策动力，F表示内置信号幅值，ω表示内置信号角频率即振子的固有频率。

为了对任意频率信号进行处理，式(1)也可写为：

(2)

固定k值，F≠0条件下，F较小时，Duffing系统相轨迹逐渐收敛到两个焦点(±1,0)中的一个。随着F的增大，相轨迹由周期运动逐渐演化为混沌运动，并在较大范围内保持混沌运动。当F大于某一个阈值Fr(Fr称为临界值)时，系统相轨迹由混沌态跃迁为大尺度周期态。

对弱信号检测时，系统内置信号幅值设定为临界值Fr，输入信号si(t)为待测信号s(t)和高斯白噪声n(t)的叠加。弱信号检测模型为：

=Frcos(ωt)+s(t)+n(t)

(3)

若取

s(t)=Acos[(ω+Δω)t+φ]

此时，Duffing系统周期驱动力r(t)可以写为：

r(t) =Frcos(ωt)+Acos[(ω+Δω)t+φ]

=B(t)cos[ωt+θ(t)]

只考虑微弱信号情况，即A≪Fr，可以得到式(4)[9]。

(4)

在没有频差，即Δω=0情况时，有：

(5)

式中，k为整数。

由式(5)可知：当-π/2+2kπ<φ<π/2+2kπ时，输入信号幅度B(t)大于系统临界值Fr，系统处于周期态；当π/2+2kπ≤φ≤3π/2+2kπ时，输入信号幅度B(t)小于等于系统临界值Fr，系统处于混沌态。

仿真实验如下。在式(3)中取参数k=0.5，初始值x(0)=0.1，y(0)=x′(0)=0，采样速率fs=40 MHz，内置信号频率f=ω/(2π)=1 MHz。待检信号取s(t)=0.01cos(2πfct)，fc=1 MHz，噪声为高斯白噪声，SNR=-30 dB。采用四阶Runge-Kutta算法对式(2)进行求解，步长为1/fs，去除暂态后得到Duffing系统在Fr=0.825 V处只有噪声以及有待检信号和噪声输入时的相轨迹，如图1所示。

(a) 只有噪声时，混沌态(a) Chaotic state with noise only

(b) 有待检信号和噪声时，周期态(b) Periodic state with signal and noise图1 微弱周期信号的Duffing系统相图Fig.1 Phase figures of Duffing system with the weak input signal

1.2 Duffing振子微弱BPSK信号时序图法检测

双极性BPSK信号的时域表达式可以表示为：

(6)

其中，ωc为信号载波频率。

若BPSK信号“1”码和“-1”码等概率出现，可以得到BPSK信号的功率谱密度：

(7)

式中，Td为基带信号码元周期。

从式(7)可以看出，BPSK信号功率谱中无载波分量，实际是抑制载波的调制信号，这给BPSK信号的检测带来难度。

由式(5)可知，当Duffing系统输入信号为微弱BPSK信号时，由于BPSK信号相位的不断变化，在Δω=0时，Duffing系统相图随着BPSK信号相位的变化处于间歇性混沌状态，周期态和混沌态交替出现。此时，只有和时间密切关联，才能准确判断间歇性混沌状态的存在，进而判断微弱信号的存在。

对混沌特性的判定，目前通常采用的方法主要包括相平面法、Lyaponov指数法和时序图法[6]。前两种方法反映的都是一段时间内Duffing系统输出的整体特征，而时序图法反映了在不同时间Duffing系统输出的状态，它和时间关联紧密。本文采用时序图法对Duffing系统输出状态进行描述。

Duffing系统输出的混沌态和周期态的包络有比较明显的差别：在混沌态下，Duffing系统输出包络起伏较大；而在周期态下，包络恒定。为了更清晰地反映时序输出包络变化，本文采用S变换提取Duffing系统输出时序包络。

S变换是在连续小波变换和短时傅里叶变换基础上发展起来的一种时频分析方法，同其他时频方法相比，具有多分辨特性、无交叉项、计算速度快等特点[10-11]。

信号h(t)的离散S变换形式为：

式中：a为常数，其决定了S变换对应的窗函数；H(·)为信号h(t)的离散傅里叶变换；N为采样点数；k=0,1,…,N-1；n=1,…,N-1。为了进一步提高运算速度，采用作者文献[11]提出的改进离散S变换算法提取Duffing系统输出时序包络。

仿真实验如下。取BPSK信号仿真采样点数N=10 000，A=0.01，fc=1 MHz，码元速率fd=1/Td=40 kHz，采样速率fs=40 MHz。Duffing系统参数的设定同图1，SNR=-30 dB。S变换中，为了更好提取f=1 MHz的Duffing系统时序输出包络和反映包络变化，S变换窗函数在f=1 MHz处的时窗宽度取为4.68μs。此时，BPSK信号与Duffing系统内置频率同频，得到Duffing系统时序输出和在f=1 MHz处S变换包络提取结果，如图2所示。

(a) 调制码元波形(a) Time wave of modulation symbol

(b) Duffing系统时序输出(b) Time sequence output of Duffing system

(c) f=1 MHz 处S变换结果(包络)(c) Extracted envelope at f=1 MHz by using S transform图2 BPSK信号Duffing振子输出时序图及包络Fig.2 Time sequence and envelope of BPSK signal’s Duffing system output

从图2可以看出，Duffing系统输出包络基本反映了BPSK信号的状态:输入信号调制码元为“1”时，Duffing系统呈周期态，输出包络恒定；输入信号调制码元为“-1”时，呈混沌态，输出包络起伏较大。

仔细观察图2(c)，可以看出：Duffing系统输出状态的变化时刻并不能与调制码元变化时刻完全对应，总是存在一个时间延迟。即在混沌态和周期态之间转换时存在过渡带。

这种延迟或过渡带的存在，对BPSK信号的判断是有害的。当BPSK调制码元速率较高时，如果此时过渡带时间过长，Duffing系统输出将很难被分辨处于周期态还是混沌态，从而严重影响BPSK信号的检测。

2 内置频率对Duffing系统性能影响分析

对于具有非线性恢复力和一个周期外力的Duffing方程的求解，文献[12-13]分析了依赖阻尼、内置频率、振幅和相差的分支结构和动力学行为，应用二次平均方法给出了周期解和三次以下谐波解，应用Melnikov方法分析了m(m>3)阶次谐波解和混沌的存在条件，得到了一些有益的结论。本节在这些结论的基础上，依据方程(1)，在阻尼系数和非线性恢复力确定的情况下，重点分析内置频率对Duffing系统性能的影响。

2.1 内置频率对过渡带的影响

在周期态和混沌态转换之间存在过渡带，对这方面研究的文献资料比较少。文献[10]和文献[14]提出了从周期态向混沌态的逆向相变检测方法，减少相变过渡带的影响，但对过渡带与什么因素有关没有进行分析。

为了求解方程(1)，可以将其看成是未扰动系统和扰动系统的组合。其未扰动系统为：

(8)

存在三个不动点：鞍点(0,0)，中心(±1,0)。方程的解为周期解。初始条件为x(0)=D，y(0)=x′(0)=0。采用文献[15]的结论，可以推导出无阻尼Duffing方程(8)的解为：

(9)

式中，

(10)

式中，

(a) 混沌态(a) Chaotic state

(b) 周期态(b) Periodic state图3 混沌态和周期态下Melnikov曲线示意图Fig.3 Melnikov curves in chaotic and periodic state

外部输入信号为BPSK信号时。在A≪Fr和Δω=0情况下，并且不考虑噪声，由式(4)可以得到式(3)的Duffing系统周期驱动力。

其中，

Duffing系统将随着F的变化在周期态和混沌态之间变换：F=F1时，系统进入周期态；F=F2时，系统进入混沌态。此时可以采用在混沌态时的零解和周期态时的最大解对应的时间差的最小值来衡量混沌态和周期态转换的过渡时间。

在混沌态时的零解对应时刻tc满足：

解得：

(11)

式中，m为整数。

在周期态时的最大解对应时刻tp满足：

解得：

(12)

式中，n为整数。

由式(11)～(12)可知，用来衡量混沌态和周期态转换的过渡带时间的tc和tp之间时间差的最小值tg为：

(13)

由式(13)可知，过渡带时间仅与Duffing系统内置频率有关。

直接分析式(13)比较复杂，为此仿真实验画出时间差tg与系统内置频率f之间的曲线。取内置频率从1～10 MHz，取样间隔0.1 MHz，其他实验条件同图1。计算式(13)可以得到tg与Duffing系统内置频率f的关系仿真结果，如图4所示。由图4可以看出，tg与ω基本成反比关系，Duffing系统内置频率越大，过渡带时间越短。

为了进一步验证上述结论，在不同频率下，采用Duffing系统对BPSK信号进行处理，得到Duffing系统输出时序包络和调制码元对应关系，如图5所示。实验中，频率分别取1 MHz、1.75 MHz、3.25 MHz和5 MHz，SNR=-15 dB，其他实验条件同图2。从图5可以明显看出，随着Duffing系统内置频率的增加，周期态和混沌态之间的过渡时间明显缩短。

图4 过渡带时间与内置频率关系Fig.4 Curve of transition zone time length and forcing frequency

(a) f=1 MHz

(b) f=1.75 MHz

(c) f=3.25 MHz

(d) f=5 MHz图5 不同频率下BPSK信号Duffing系统输出时序包络与码元变化对应关系Fig.5 Duffing system output envelopes of BPSK signal at different frequencies

2.2 内置频率对灵敏度的影响

在Duffing系统幅频特性研究方面，文献[18]给出了周期态下系统内置频率作为因变量的频率-幅度表达式，并没有给出输出谐波幅度作为因变量的表达式，无法得到系统幅频响应曲线。但其首次分析了Duffing系统最佳内置策动力频率，估算出最佳策动力频率约为0.5 rad/s。

当Duffing系统有同频外部激励信号并处于周期态时，方程(1)的谐波解主要包含一次谐波，其周期解如式(14)所示[12-13]。

(14)

式中，(x0,0)为未扰动系统式(8)的中心。

应用谐波平衡法对处于临界态的方程(1)进行分析，可以得到Duffing系统输出幅度a与频率ω的关系[19]。

a2[ω2+1-0.75a2]2+k2ω2a2=F2

(15)

令z=a2，同时将ω2看作系数，则式(15)可以看作一个标准的一元三次方程：

Ez3+Bz2+Cz+D=0

(16)

式中，E=0.752，B=-1.5(ω2+1)，C=ω4+(2+k2)ω2+1，D=-F2。

由盛金公式求解式(16)，令

在式(15)中取k=0.5，F=Fr=0.825 V。可以得到，当ω>1 rad/s时，一元三次方程根的判别式Δ=q2/4+p3/27>0，式(15)的实数解为：

(17)

当ω≤1 rad/s时，根判别式Δ存在小于0的情况。Δ<0时，式(15)的解为：

(18)

综合式(15)、式(17)和式(18)，可以看出，不同频率下Duffing系统输出幅度不同，即Duffing系统对弱信号的灵敏度是内置频率ω的函数。而当ω>1 rad/s时，Duffing系统输出幅度a是内置频率ω的递减函数。

仿真实验如下。ω取值范围为0.1～10 rad/s，间隔0.1 rad/s，F=Fr=0.825 V，k=0.5。对式(15)求解，得到Duffing系统输出幅度a与内置频率ω的关系曲线，如图6所示。

图6 周期态时Duffing系统输出幅度与内置频率关系曲线Fig.6 Curve of the output amplitude and the forcing frequency of Duffing system in periodic state

从图6可以看出，Duffing系统的幅频动态响应特性在低频段具有良好的特性，并且在ω≈0.6 rad/s处取得最大值，获得最佳幅频响应。而后随着内置信号频率的增大，幅频动态响应特性将变差。Duffing系统只有小频率参数条件下才能有较好的动态特性和检测效果。

3 对任意频率微弱BPSK信号盲检测

由前文可知，对微弱BPSK信号检测，一方面，BPSK信号本身相位变化，引起Duffing系统输出在周期态和混沌态之间切换，此时要求过渡带时间尽量短，需要提高Duffing系统内置频率；另一方面，由Duffing系统本身幅频响应特性决定，只有在低频段才能获得良好的动态特性。这两方面因素构成了一对矛盾。

为了将Duffing振子应用于任意频率的微弱信号检测，文献[20-21]等提出了变尺度的微弱信号检测方法，对待测信号在时间尺度上进行压缩或放大。但这种方法需要知道待测信号频率，不能应用于未知频率信号的盲检测。文献[18,22]提出了变步长型搜索检测未知信号的方法，其实质是采用变尺度方法，在|Δω/ω|<0.03频段范围内实现对未知频率信号的检测。

变尺度方法的基本思路是：对于频率为ω，采样速率为fs的待测信号，取微分方程数值计算的步长为1/fs。若将数值计算步长人为增大R倍，相当于相应信号的频率被压缩至原值的1/R，从而获取Duffing系统良好的低频特性。

对于未知频率信号，假设与内置频率存在频差Δω≠0，由式(4)可以得到：

φ一定，-π/2+2kπ<Δωt+φ<π/2+2kπ时，输入信号幅度B(t)大于系统临界值Fr，输出周期态；π/2+2kπ≤Δωt+φ≤3π/2+2kπ时，输入信号幅度B(t)小于等于系统临界值Fr，输出混沌态。令Δω=2πΔf，B(t)将以周期2π/Δω(即1/Δf)大于和小于临界值Fr，使得Duffing系统也将以周期1/Δf在周期态和混沌态之间转换，且周期态和混沌态持续时间都为π/Δω(即1/(2Δf))，系统进入周期性混沌状态。

对BPSK信号，无频差情况下，信号码元的变化将导致Duffing系统输出处于间歇性混沌状态。如果再考虑频差的影响，则Duffing系统输出肯定处于间歇性混沌状态。

从BPSK信号相位变化影响来说，即码元速率影响方面，频差引起的周期态或混沌态持续时间tspan=1/(2|Δf|)与调制信号码元长度td=1/fd满足2tspan>td，即|Δf|

从频差影响方面来说，输入信号与Duffing系统内置信号频差越小，对间歇性混沌状态的贡献越小；频差越大，贡献越大。频差太大，将使混沌态和周期态持续时间很短，无法找到持续一段时间的恒定包络，不能判断系统状态。文献[23]表明，当|Δf|/f≤0.03时，有规则的周期混沌现象才可以被辨认出来。

综合以上两点，Duffing系统识别BPSK信号的工作带宽B=2|Δf|应该满足：

B

(19)

通过变尺度算法，内置频率为f的Duffing振子可以实现对f±(B/2)范围内的任意频率BPSK信号的盲检测，并且获得良好的幅频动态特性。此时，对Duffing系统内置频率的选择依据主要考虑过渡带时长。转换时间过长，将会影响判断系统输出状态，因此，内置频率f不宜选取过低。

仿真实验如下。取BPSK信号频率fc=3.505 MHz，Duffing系统内置频率f=3.5 MHz，其他实验条件同图2。得到Duffing系统时序输出和在f=3.5 MHz处S变换包络提取结果，如图7所示。可以看出，受频差和相位变化双重影响，BPSK信号Duffing系统输出状态复杂，其周期态、混沌态之间切换时刻与码元变化时刻不能一一对应，甚至有时状态截然相反，但仍然可以清晰分辨出系统输出为间歇性混沌状态，从而判断BPSK信号的存在。

图7 有频差BPSK信号Duffing振子输出时序包络Fig.7 Duffing system output envelope of BPSK signal having frequency difference

在整个频段范围内判断是否有微弱BPSK信号存在，可以将处理频段分为n个通道，每个通道控守的带宽为B，采用通道分段控守，实现对整个频段的全覆盖。建立Duffing振子微弱BPSK信号盲检测模型，如图8所示。若任一通道内有间歇性混沌状态输出，则判断在对应通道内有微弱BPSK信号存在。

图8 基于Duffing振子和S变换微弱BPSK信号盲检测模型Fig.8 Blind detection model of weak BPSK signal by using Duffing oscillators and S transform

为验证变尺度方法的有效性，分别以原始方法和变尺度方法观测Duffing系统的幅频特性。仿真接收机中频输出的BPSK信号载频以0.5 MHz步进从0.5～5 MHz变化，共计10批实验。信道输入高斯白噪声，在每个频率点上(即每批实验中)信噪比以1 dB步进从-10～-45 dB变化，每批共计36组实验。取BPSK信号码元速率fd=1/Td=40 kHz，设定Duffing系统过渡时间阈值为Th=Td/2，当信噪比从高到低变化时，求得满足过渡带时长tg

图9 变尺度和原始方法Duffing系统检测性能对比Fig.9 Comparison of detection performance by using original and scale-transformation Duffing oscillator

4 结论

针对Duffing系统输出混沌态和周期态之间存在过渡带问题，分析发现过渡带时间和Duffing系统内置频率有关，给出了过渡带时间与内置频率之间关系的表达式，指出内置频率越大，过渡带时间越短。

针对周期态下Duffing系统输出与内置频率的幅频响应问题，给出了输出幅度作为因变量、内置频率作为自变量的关系表达式，指出ω≈0.6 rad/s处为最佳幅频响应点，此后随着内置频率增大，系统幅频特性下降。

将变尺度方法与检测阵列结合，分析了内置频率和检测阵列控守频段范围的选取原则，建立了基于Duffing振子和S变化的任意频率微弱BPSK信号盲检测模型。