利用阵发混沌现象测定未知信号参数

翟笃庆 刘崇新 刘 尧 许 喆

(西安交通大学电气工程学院,西安 710049)

利用阵发混沌现象测定未知信号参数

翟笃庆†刘崇新 刘 尧 许 喆

(西安交通大学电气工程学院,西安 710049)

(2009年5月12日收到;2009年6月17日收到修改稿)

利用Duffing方程对频率的极端敏感性产生阵发混沌现象,研究了一种利用该现象定量检测未知的微弱周期信号的各项参数的新方法,通过理论分析和实例仿真证明了该方法的可行性,并针对仿真结果提出了改进措施,提高了检测精度.

Duffing方程,阵发混沌,检测信号参数

PACC:0545

1.引言

信号处理领域的一个重要研究课题是对未知信号的各项参数的测定,包括频率、幅值和相位,这在工程应用中有着重要的实际意义.近年来,混沌同步的思想被科学工作者们应用于未知参数的辨识中.文献[1]利用状态观测器实现了不确定的Lü系统的混沌控制和参数辨识.文献[2]对一类不确定参数的复杂网络提出了自适应同步方法.文献[3—6]分别以不确定参数的Chen系统、超混沌Chen系统、Lorenz-like系统和Rössler系统为例,在实现自适应同步的同时,辨识了未知参数.文献[7]则以Gray-Scott系统产生时空混沌,实现混沌同步和参数辨识.文献[8]利用自适应法,实现了Duffing等系统的混沌同步和参数辨识.文献[9]研究了不确定Liu系统的混沌同步和参数辨识.文献[10]则把自适应同步和参数辨识的原理和思想应用到混沌保密通信中,获得加密系统的参数.这些成果不仅实现了不确定混沌系统之间的同结构和异结构同步,完善了混沌同步的理论,也为不确定系统的未知参数辨识找到了方法.但目前该方法的研究主要集中于自治系统的未知系数(包括非自治系统的周期信号的幅值)的辨识,对于非自治系统的频率和相位参数,尚没有行之有效的确定方法.本文针对Duffing方程,基于另一种思想,寻求找到确定Duffing系统中时间变量的频率和相位参数的方法.另外,人们在对混沌理论的研究中发现:一类混沌系统在一定条件下对小信号具有极端的敏感性而同时对噪声具有很强的免疫力,这使得此类混沌系统在信号检测领域中具有广阔的应用前景.著名的Duffing振子就具有这种特性,当其处于混沌分岔区域附近时,对临近频率的微弱信号敏感而对强噪声免疫.这使得该方法在强噪声背景下同样有效.文献[11]在对Duffing系统的研究过程中,发现了由于对频率参数的敏感性,导致系统状态时而周期,时而混沌,像呼吸一样时断时续,并且“呼”与“吸”都有固定的时间跨度,故称其“breather”.文献[12]提出了Duffing方程右侧的周期策动力幅值的敏感性,即存在某个临界值,幅值稍大于该临界值时就会呈现大尺度周期状态,小于该临界值时则呈现混沌状态,并把对于该敏感性的检测用于测定某已知频率的信号幅值.文献[13]通过计算周期状态时一段时间内的过零点次数来确定未知信号频率.文献[14]和[15]提出在Duffing方程右侧加上策动力频率附近的微弱信号后,会产生阵发混沌现象,策动力幅值会随时间变化.当幅值大于文献[12]中的临界值时,呈现大尺度周期运动状态,小于该临界值时,呈现混沌状态,这样就会出现时而周期,时而混沌的运动状态,把这种现象称为阵发混沌现象.通过观察该现象是否发生来判断待测信号中是否存在与策动力频率相近的周期信号.文献[16]对一个带有外部扰动的Duffing方程进行了详尽的分析,包括其最大Lyapunov指数、分叉情况和阵发特性等.

在前人研究的基础上,本文在Duffing方程的右端施加谐振扰动,产生阵发混沌现象,分析了扰动后的系统状态的变化情况.通过对该现象的观测和测定,确定施加的待测信号的各项参数,包括频率、幅值和相位.仿真实验得到了较为精确的结果,证明了该方案的可行性,且对多频信号同样适用.

2.阵发混沌现象的理论分析和仿真

经典的Duffing方程为

式中x(t)为状态变量,k为阻尼比,fcos(t)为周期策动力,-x(t)+x3(t)为非线性恢复力.改写为如下动力学方程:

当k=0.3时,改变参数f的取值,观察系统的x-y相图,会发生如图1所示的状态变化.使系统从混沌状态转变到大尺度周期状态的参数f的临界值Fc∈(0.53492,0.53493).从图2所示分叉图中也可看出状态变化的临界值范围.

令t=ω τ,f=F0,(1)式可改写为

在(3)式右侧加上频率ω附近的待测周期信号

对(4)式右侧化简

其中

图1 参数f改变引起的系统x-y相图变化 k=0.3.(a)f=0.2,(b)f=0.35,(c)f=0.53492,(d)f=0.53493

图2 Duffing方程x变量随参数f变化的分叉图 k=0.3.(a)f∈(0,10),(b)f∈(0.2,0.65)

观察(5)式与(1)式,当ω=1时,二者差别在于(5)式中多出了一个参数相位角.文献[17]讨论了相位为随机噪声时对系统状态变化的影响.这里的相位是与幅值相关的函数,考虑其影响系统状态的情况.

由(6)式得到

把(8)和(9)式代入(7)式中,令θ(τ)=g(F(τ)),得

图3 F(τ)引起的x-y相图变化 (a)F(τ)=0.5413((9)式取正号),(b)F(τ)=0.5414((9)式取正号),(c)F(τ) =0.537((9)式取负号),(d)F(τ)=0.5371((9)式取负号)

取F0=0.52,A=0.1,k=0.3,ω=1,由(6)式考虑到F(τ)∈(F0-A,F0+A)=(0.42,0.62),观察系统状态变化情况(见图3和图4).发现存在临界值Fc,当(9)式取正号时,Fc∈(0.5413, 0.5414),当(9)式取负号时,Fc∈(0.537, 0.5371),当F(τ)>Fc时,系统处于周期状态;当F(τ)

由于F(t)的取值会随着时间t变化,时而大于Fc,时而小于Fc,(10)式所描述的系统就会呈现混沌状态和周期状态交替出现的情况,这就是阵发混沌现象.令F0=0.52,A=0.04,ω=1,Δ ω=0.01,通过时域信号图可以清晰地观察到阵发混沌现象(见图5(a)).

图4 (10)式中变量x随F(τ)变化的分叉图 (a)(9)式取正号,(b)(9)式取负号

图5 阵发混沌现象 (a)Δ ω=0.01,(b)Δ ω=0.05

这里需要注意的是Δ ω不可以过大或过小.过小时,F(t)变化缓慢,导致阵发混沌现象的每个混沌状态和周期状态跨度过大,要检测出相对频差就需要大大增加采样和计算的工作量;另外,由于系统完成状态变化通常需要1个周期以上的时间,当Δ ω过大时,F(t)变化过快,系统完成相变的过程中没有足够长时间的激励,系统无法很好地响应如此快的变化,阵发混沌现象失去原有的规律性.通过仿真实验,我们一般取Δ ω∈(0.005,0.04).

3.利用阵发混沌现象测定待测信号参数

3.1.测定信号频率

加上待测信号后的策动力幅值f在临界值Fc附近随时间变化

若要检测频率ω′∈(ω1,ω2)的待测信号,先选取合适的公比1+Δ ω,从ω1开始以公比1+Δ ω取值作为周期策动力的频率,检测阵发混沌现象是否发生,如果发生就用上述方法算出ω′,若一直没有发生则说明不存在该频率区间内的周期信号.

3.2.测定信号相位

设当cos(ω Δ ωt+φ)=a,a∈0,1时,F(t) =Fc,则当cos(ω Δ ωt+φ)>a时,系统处于大尺度周期状态,对应图6中的T1时间段;当cos(ω Δ ωt+ φ)

设(4)式中的周期策动力为F0cos(ωt+φ1),待测信号Acos(ω(1+Δ ω)t+φ2),二者相位差Δ φ =φ2-φ1,但F(t)和θ(t)均有变化

假设相位差Δ φ=0,当t=0时,cos(ω Δ ωt+ Δ φ)=1,F(t)=F0+A>Fc,随着时间变化,t= T1时,cos(ω Δ ωt+Δ φ)减小到a,F(t)从F0+A减小到Fc,系统从大尺度周期状态进入混沌状态.同样的,如果cos(ω Δ ωt+φ)从a增加到1,相应的F(t)从Fc增加到F0+A,对应的,也需要时间T1.这样,我们可以得到如下结论:当时,cos(ω Δ ωt+Δ φ)=1,系统处于大尺度周期状态且正处于该状态的时间中点,即该状态恰好完成了一半.图7所示的仿真结果,以时间零点观察,图7(b)图中零点恰好为周期状态进行到一半的时间点,证明了该结论.

图6 系统分别处于大尺度周期和混沌状态的时间段

找出系统由混沌进入大尺度周期状态的时间t1和由大尺度周期状态进入混沌状态的时间t2,,

待测信号相位φ2=φ1+Δ φ.

3.3.测定信号幅值

由图6所示关系得

其中,T1,T2分别是大尺度周期状态和混沌状态的持续时间.

对于(11)式描述的幅值随时间变化的函数, F(t)=Fc为临界值t为取到临界值的时间点,所以得到函数关系

通过仿真结果找到T1,T2的值,代入(16)式中,解该方程即可求得待测信号幅值A.

4.仿真实验

先预估信号频率ω′∈(45,55),从ω=45开始以公比1.03取值作为周期策动力的频率,直到观测到阵发混沌现象,若没有发生阵发混沌现象,则不存在该频率范围的周期信号.图8所示的是周期策动力频率取不同值时(周期策动力相位取为0,幅值F0=0.53),x(t)的时域信号.

图7 Δ φ不同时的时域信号图

图8 不同ω时的时域信号图 (a)ω=45,(b)ω=47.74,(c)ω=49.17

通过时域信号图可以看出,ω=47.74时阵发混沌现象周期尚不明显,ω=49.17时发生了明显且有规则的周期性阵发混沌现象.以ω=49.17作为周期策动力频率,通过Matlab仿真结果(见表1)找出由混沌状态向大尺度周期状态转变的时间点.

表1 混沌状态向大周期状态转变的仿真结果

通过比较定向过零点间距,可以看出自t= 3.29后过零点间距变化很小,且渐渐趋于稳定,故认为t=3.29时F(t)开始超过Fc,判断t=3.29为混沌状态向大尺度周期状态转变的时间点.

同样的方法找出t=11.08,18.81,26.10, 33.76,41.18,48.84,56.38,64.17,71.90也是混沌状态向大尺度周期状态转变的时间点.求大周期的平均值T=7.62333,ω Δ ω=0.82420,待测信号频率ω′=ω(1+Δ ω)=49.994.

通过仿真结果比较过零点间距(见表2),发现自t=8.49后过零点间距开始逐渐失稳,故认为t=8.49时F(t)开始小于Fc,判断t=8.49为大尺度周期状态向混沌状态转变的时间点,则t=5.89为大尺度周期状态的时间中点,因此按照3.2节中所述方法,待测信号相位

表2 大尺度周期状态向混沌状态转变的仿真结果

同样的方法找出t=16.03,23.50,31.17, 38.71,46.38,53.72,61.39,69.06,76.60,84.07也是大尺度周期状态向混沌状态转变的时间点,求出各个大尺度周期状态的时间中点,然后测定相位

求得平均值φ′=0.4931π.

大尺度周期状态结束的时间(见表2)减去开始的时间(见表1),得到大周期状态的持续时间T1,多次求T1再求其平均值可得T1=4.994,混沌状态的持续时间T2=T-T1=2.629,代入(16)式中,解得待测信号幅值A=0.5070.其中临界值Fc= 0.537,周期策动力幅值F0=0.53.

实例2 待测信号f′(t)=0.8cos(100t+π)

预估待测信号频率ω′∈(90,110),周期策动力频率从ω=90开始以公比1.03递增,发现ω= 98.35时出现规则的阵发混沌现象.设定周期策动力频率为ω=98.35,步长Δt选取为0.005,通过仿真结果找出t=4.335,8.105,12.09,15.71,19.6, 23.30,27.27,30.975,34.745,38.67为混沌状态向大尺度周期状态转变的时间点.求大周期的平均值待测信号频率ω=ω(1+Δ ω)=99.997.

通过仿真结果找出t=7.335,11.07,14.97, 18.705,22.54,26.34,30.175,33.975,37.71, 41.61.多次求相位再求平均值得φ′=0.9087π.

大尺度周期状态持续时间T1=2.94,混沌状态的持续时间T2=T-T1=0.875,代入(16)式中,解得待测信号幅值A=0.8032.其中临界值Fc= 0.537.

对于该多频信号,先检测ω=50附近,设置周期策动力频率ω=49,周期策动力幅值F0=0.53,通过仿真结果求得大周期平均值T=6.2844,ω Δ ω =0.9998,待测信号频率ω′=ω(1+Δ ω)= 49.9998.多次求相位再求平均值φ′=1.4294π.大尺度周期状态持续时间T1=3.871,混沌状态的持续时间T2=T-T1=2.4134,代入(16)式中,解得待测信号幅值A=0.3923.

为了提高其他频率分量的检测精度,把ω=50的分量从原信号中减去,即f″(t)=f′(t)-0.3923cos(49.9998t+1.4294π).然后再检测ω= 100附近,设置周期策动力频率ω=98,周期策动力幅值F0=0.53,通过仿真结果求得大周期平均值T =3.1422,ω Δ ω=1.9996,待测信号频率ω′=ω(1 +Δ ω)=99.9996.多次求相位再求平均值φ′= 0.9310π.大尺度周期状态持续时间T1=1.8585,混沌状态的持续时间T2=T-T1=1.2837,代入(16)式中,解得待测信号幅值A=0.3177.

f‴(t)=f″(t)-0.3177cos(99.9996t+ 0.9310π),再设置周期策动力频率ω=195,周期策动力幅值F0=0.53,通过仿真结果求得大周期平均值T=1.2442,ω Δ ω=5.0500,待测信号频率ω′=ω(1+Δ ω)=200.0500.多次求相位再求平均值φ′=0.3044π.大尺度周期状态持续时间T1= 0.7443,混沌状态的持续时间T2=T-T1= 0.4999,代入(16)式中,解得待测信号幅值A= 0.2303.

5.方法改进

5.1.相位测量的误差分析及改进

通过实例仿真发现该方法测量相位的结果误差较大,其原因如下:

1)系统从某个状态向另一个状态转变的过程中,需要一定时间的激励,一般情况下需要超过一个周期的稳定激励才能发生状态变化,使得测量结果较实际值有一定滞后.

2)该方法通过比较过零点间距判断状态,决定了系统状态的起始点只能在过零点,而实际状态变化则会在任意时间点完成,这也使得结果较实际值滞后.

3)(7)式中ω Δ ω决定了误差大小,因此减小Δ ω,即减小每次搜索频率的范围,可以减小误差.

表3是对实例2的参数进行改动后得到的实验数据,表明了相位滞后的大小主要由ω Δ ω决定,受待测信号幅值和相位影响不大,因此在ω Δ ω测定后,可通过仿真实验估算出滞后相位的大小.测量结果加上该滞后相位,即可得到较精确的相位值.如实例3,可对滞后相位的大小求平均值Δ φ= 0.0994π,改进后的相位测量值φ″=φ′+Δ φ= 1.0081π,大幅度提高了相位测量的精度.

表3 滞后相位与待测信号参数的关系

5.2.幅值测量方法的改进

对待测信号幅值的测定,是建立在精确测得大尺度周期状态和混沌状态的起始点的基础上的,这样才可以准确得到大尺度周期状态和混沌状态的持续时间,代入(16)式中才能精确地求解出待测信号幅值.因此如何精确测出系统状态起始点,是提高幅值测量精度的关键.

1)离散化时尽量选择较小的步长,如果步长过大,则过零点间距的变化对系统状态变化不够敏感,影响状态起始点的测量精度.但考虑到步长过小会影响算法的效率,因此如何选择步长,要视实际情况而定.

2)当待测信号幅值较小时,通过过零点间距来判断系统状态变化的方法可能会失效,因为系统在F(t)稍大于Fc时,混沌状态向周期状态转变时存在较长时间的过渡过程(见图9,观察过零点间距,趋于等距要在两个周期以后),而系统在该时间点停留时间内无法完成状态的变化,按照过零点间距的判定方法,误认为该状态仍是混沌状态.另外,当A-F0>Fc时,(6)式中F(t)∈(A-F0,A+F0),系统将一直处于周期状态,阵发混沌现象消失.

表4是对实例3的幅值参数进行改动得到的实验数据,它反映了幅值过小或过大使得该方法失效的现象.因此使用该方法前,需保证待测信号幅值处于有效区间内,如果不满足,则需先预估幅值,再放大或缩小使之处于有效区间.该区间由信号频率ω′和周期策动力频率ω决定,可在测出频率后通过仿真实验确定.表4中ω=98,ω′=100.

图9 F(t)=0.55时的时域信号图

表4 幅值测量误差与实际幅值的关系

6.结 论

本文提出了一种利用阵发混沌现象来检测未知信号各项参数的方法,包括频率、幅值和相位,从而完整地还原待测信号.该方法的具体步骤如下:

1)将待测信号加到Duffing方程右侧,取合适的步长对周期策动力和待测信号之和进行采样.

2)在待测信号可能存在的频率区间内取合适的公比,按等比区间进行遍历搜索,若检测到阵发混沌现象发生,则按3.1节所述方法确定信号频率.

3)分别按照3.2节和3.3节所述方法测定待测信号幅值和相位.

4)对于多频信号,每次将所得前一频率分量的信号从方程右侧减去,再在下一个频率分量可能存在的区间进行搜索,如此重复上述步骤,确定所有频率分量信号的各项参数.

理论分析和实例仿真都证明了该方法的可行性.Duffing方程的性质决定了该方法具有以下特点:

1)对信号中各频率分量要先预估,可能存在频率的区间要取到合适的公比进行遍历搜索,才能防止频率的遗漏.相对频率差一般取0.005<Δ ω< 0.04,过大或过小都会使阵发混沌现象消失,所以搜索时的公比也要取到该区间.

2)测定幅值之前也需要先进行预估,然后对待测信号放大或缩小,使幅值处于有效区间内,该区间由周期策动力频率ω和待测信号频率ω′决定,二者间相对频差Δ ω越大,该区间也越大.幅值测定后,可以把测量结果从原信号中减去,判断阵发混沌现象是否消失,若消失则证明结果正确.

3)由于系统状态变化需要一定时间的激励才能完成,导致相位的测量结果较实际值有不同程度的滞后,该滞后相位的大小主要由ω Δ ω决定,因此可以在测定频率后通过仿真确定滞后相位的大小,加上测量结果,即得到较精确的相位值.

相较于混沌同步完成参数辨识的思想,本文的研究方法有以下优点:

1)速度快.一般在一个大周期后,如实例1,在t =11 s左右即可完成参数测定,而前者一般要40 s以上(见文献[1]).

2)抗噪性好.一般可在信噪比达-100 dB以上的条件下检出信号参数(见文献[9]).

3)该方法能完整地测出待测信号的各项参数,这是目前前者所无法完成的,但在幅值和相位的测量精度方面还有待提高.

[1]LüJ H,Lu J A 2003Chaos Soliton.Fract.17 127

[2]Zhou J,Lu J A,LüJ H 2006IEEE Transactions on Automatic Control51 652

[3]Wang X Y,Wu X J 2006Acta Phys.Sin55 605(in Chinese) [王兴元、武相军2006物理学报55 605]

[4]Wang X Y,Wu X J 2008Int.J.M od.Phys.B 22 1015

[5]Wang X Y,Wang Y 2008Int.J.M od.Phys.B 22 2453

[6]Wang X Y,Li X G 2008Int.J.M od.Phys.B 22 4069

[7]LüL,Li Y,Guo Z A 2009Sci.Chin.Ser.G 39 77(in Chinese)[吕 翎、李 义、郭治安2009中国科学,G辑39 77]

[8]Kim J H,Park C W,Kim E,Park M 2005Phys.Lett.A 334 295

[9]Hu J,ZhangQ J 2008Chin.Phys.B 17 503

[10]Ren H P,Han C Z,Liu D 2008Chin.Phys.B 17 1202

[11]Yang J G,Qu Z L,Hu G 1996Phys.Rev.E 53 4402

[12]Li Y,Yang B J,Shi Y W 2003Acta Phys.Sin.52 526(in Chinese)[李 月、杨宝俊、石要武2003物理学报52 526]

[13]Li Y,YangB J,Lin H B,Liu X H 2005Acta Phys.Sin.54 1994(in Chinese)[李 月、杨宝俊、林红波、刘晓华2005物理学报54 1994]

[14]LiA,Wang S 2004Electr.Quality06 06(in Chinese)[李安、王 殊2004电子质量06 06]

[15]Zhao X Y,Wei M X,Liu J H 2004Journal of Test and Measurement Technology18 16(in Chinese)[赵向阳、魏明献、刘君华2004测试技术学报18 16]

[16]Huang J C,Jing Z J 2007Chaos Soliton.Fract.40 1449

[17]Li S,Xu W,Li R H 2005Acta Phys.Sin.55 1049(in Chinese)[李 爽、徐 伟、李瑞红2005物理学报55 1049]

PACC:0545

Deter m ination of the parameters of unknown signals based on inter m ittent chaos＊

ZhaiDu-Qing†Liu Chong-Xin Liu Yao Xu Zhe

(Institute of Electrical Engineering,Xi’an Jiaotong University,Xi’an 710049,China)

12 May 2009;revised manuscript

17 June 2009)

We use the extreme frequency sensitivity ofDuffing’s equation to produce intermittent chaos(we call it“breather”) and propose a new method to quantitatively detect the parameters of unknown weak periodic signals.The theoretical analysis and instance simulation have proved its feasibility.We also put forward a way to improve the detection results and enhance the accuracy.

Duffing’s equation,intermittent chaos,detecting the parameters of signals

†E-mail:zdq850831@sina.com

†E-mail:zdq850831@sina.com