Duffing

Duffing（共5篇）

Duffing 篇1

摘要：在研究Duffing系统的基础上,为了利用Duffing系统进行弱信号检测,必须实现一种可靠的Duffing系统模型。对3种基于欧拉方程的模型实现方式进行了对比,通过分析欧拉方程的几何意义,并结合DSP Builder的特点,提出了一种易于工程实现的改进欧拉形式。通过对四种算法进行仿真实验得出,运用不同的状态模型建立Duffing系统会改变振子相变的临界策动力幅值,但不会对系统固有存在的状态产生影响。

关键词：Duffing系统,欧拉方程,DSP Builder,振子相变

0 引 言

混沌理论的飞速发展为自然科学领域提供了一种全新的研究手段。混沌理论使人们认识到随机输入并不是造成系统输出不规则性的惟一根源,某些确定的非线性系统也能产生非常不规则的运动。由于混沌理论在信息科学、医学、生物、机械等领域具有很大的应用潜力及发展前景,再加上计算机技术的日益发达,使得混沌已成为学术研究热点,而Duffing系统作为混沌学研究中经典的非线性系统,因其在弱信号检测[1,2,3]的优越性,得到了航空航天、通信、机械、电力等各个领域的广泛应用。

可靠的Duffing系统模型是信号检测能否顺利实现的前提和基础,本文通过分析欧拉方程的几何意义,构建了几种不同Euler算法下的Duffing振子模型并分别进行了探讨。提出了一种改进的Euler算法建立的Duffing振子模型,主要目的是为Duffing系统的FPGA实现创造前提和基础,用到的主要实验工具是DSP Builder软件。

1 Duffing混沌振子的欧拉方程实现

1.1 Duffing方程

Duffing方程是描述软弹簧振子弱阻尼运动的非线性系统。随着外界激励的变化,Duffing系统的状态会在大周期态和混沌状态之间变化[1]。Holmes型Duffing方程的标准形式如式(1)所示:

$X^{″}+k{X}^{\prime }-X(t)+X{}^3(t)=\gamma \cos t(1)$

式中:γcos t是驱动力,k是阻尼比,-X(t)+X3(t) 是非线性回复力。

在k固定的情况下 (k=0.5),系统状态将随γ 的变化出现有规律的变化。

1.2 欧拉方程实现Duffing系统的几何意义

运用欧拉方法[4]实现Duffing振子的状态模型具有鲜明的几何意义,首先,Duffing振子可以写成式(2)所示的微分方程形式,

$y^{\prime }=f(t,y),a\le t\le b(2)$

图1给出了欧拉公式几何解释的示意图,由式(2)中右端函数f(t,y)在(t,y)平面上,一个带形区域R(a≤t≤b,-∞<y<∞)内可确定一个方向场。因此,求解Duffing方程的初值问题,在几何上相当于需要确定一个通过初始点(t0,y0)的曲线y=y(t),使曲线上每点的切线均与已知向量场在该点具有一致的方向。由Duffing系统的初始点(t0,y0)出发求取(t1,y1)(其中t1是已知的等分点值)点的过程,实质上就是从(t0,y0)点出发,沿着方向场在此点的方向f(t0,y0)绘一线段,它与直线t=t1的交点便是所求的(t1,y1)点。同样,从(t1,y1)点出发,沿方向场在这一点的方向f(t1,y1)绘一线段,其与直线t=t2的交点便是所求的(t2,y2)点。依此类推,从(tk-1,yk-1)点出发,最终可达到所需的(tk,yk)点。若以A0,A1,A2,…,Ak,…,An分别表示点(t0,y0),(t1,y1),(t2,y2),…,(tk,yk),…,(tn,yn),其横坐标间隔为h,然后将这些点分别用直线段连接起来,便得到折线$\stackrel{-}{{\displaystyle A{}_{0}A{}_1}},\stackrel{-}{{\displaystyle A{}_{1}A{}_2}},\cdots ,\stackrel{-}{{\displaystyle A{}_{k-1}A{}_k}},\cdots ,\stackrel{-}{{\displaystyle A{}_{n-1}A{}_n}}‚$并记作P0。因此,欧拉方程实现Duffing振子的几何意义就是以折线P0近似代替通过初始点(t0,y0)的解函数曲线[5]y=y(t)。

同理,以h/2为步长,按照上面的方法,从(t0,y0)出发,便可得到点列(t${}_1^{(1)}$,y${}_1^{(1)}$),(t${}_2^{(1)}$,y${}_2^{(1)}$),…,(t${}_k^{(1)}$,y${}_k^{(1)}$),…,(t${}_n^{(1)}$,y${}_n^{(1)}$),用点集A${}_0^{(1)}$,A${}_1^{(1)}$,…,A${}_k^{(1)}$,…,A${}_n^{(1)}$表示。然后将点集用直线段连接起来,便可以得到一条折线P1,为$\overline{A{}_0^{(1)}A{}_1^{(1)}},\overline{A{}_1^{(1)}A{}_2^{(1)}},\cdots ,\overline{A{}_{k-1}^{(1)}A{}_k^{(1)}},\cdots ,\overline{A{}_{n-1}^{(1)}A{}_n^{(1)}}$。依此类推,以h/2k为步长,便可以得到折线Pk为$\overline{A{}_0^{(k)}A{}_1^{(k)}},\overline{A{}_1^{(k)}A{}_2^{(k)}},\cdots ,\overline{A{}_{k-1}^{(k)}A{}_k^{(k)}},\cdots ,\overline{A{}_{n-1}^{(k)}A{}_n^{(k)}}$。其中A${}_0^{(k)}$(k=1,2,…)代表图中的A0。这样便得到一个折线序列{Pk(k=0,1,2,…)},当k越大(即步长越小)时,折线越接近于Duffing方程解函数y=y(t)的曲线[6]。

1.3 前向欧拉方程构建Duffing系统

由上面的分析可以构建Duffing系统欧拉模型的状态方程,通过导数化差商、数值积分法中的矩形公式,以及运用泰勒方法展开等途径,可以将随机微分方程化成一个统一的公式[7]:

$\{\begin{array}{l}y{}_0=\eta \\ y{}_{k+1}=y{}_k+hf(t{}_k,y{}_k),k=0,1,2,\cdots ,n-1\end{array}(3)$

式(3)称为前向欧拉方程。由式可以看出,根据前一个点上的函数值yk可以直接计算出后一个点上的函数值yk+1。因此,式(3)是一个显式公式,只要初值y0已知,就可以递推计算出微分方程的数值解序列{yk(k=0,1,2,…,n)}。

由式(1)和式(3)经过推导,可以得到Duffing系统的前向欧拉方程如式(4)所示:

$\{\begin{array}{l}x(n+1)=x(n)+\omega hy(n)\\ y(n+1)=(1-0.5\omega h)y(n)+\\ \omega h[x(n)-x{}^3(n)]+\gamma \omega h\text{c}\text{o}\text{s}(h\omega n)\end{array}(4)$

1.4 后向欧拉方程构建Duffing系统

在数值积分法中采用矩形公式,可以得到如式(5)所示的方程形式[8]:

$\{\begin{array}{l}y{}_0=\eta \\ y{}_{k+1}=y{}_k+hf(t{}_k,y{}_{k+1}),k=0,1,2,\cdots ,n-1\end{array}(5)$

由式(5)可以看出,根据前一个点上的函数值yk尚不能直接计算出后一个点的函数值yk+1,因为在公式右端的第二项中仍包含有需要求的函数值yk+1。因此,式(5)是一个隐式公式,若初值y0已知,由公式逐次递推地解出其数值解的序列{yk(k=0,1,2,…,n)}。

后向欧拉公式写成(6)的形式:

$\{\begin{array}{l}y{}_0=\alpha \\ y{}_k=y{}_{k+1}-hf(t{}_{k+1},y{}_{k+1}),k=0,1,2,\cdots ,n-1\end{array}(6)$

据此可以根据k+1时刻的yk+1直接计算出k时刻的函数值yk。

由式(1)和式(6),得到Duffing系统的后向欧拉方程如式(7)所示:

$\{\begin{array}{l}x(n+1)=x(n)+\omega hy(n+1)\\ y(n+1)=(1-0.5\omega h)y(n)+\omega h\\ [x(n)-x{}^3(n)]+\gamma \omega h\text{c}\text{o}\text{s}[(h\omega n)]\end{array}(7)$

1.5 梯形构建Duffing系统

在数值积分法中采用梯形公式[9],就可以得到如式(8)的方程形式:

$\{\begin{array}{l}y{}_0=\eta \\ y{}_{k+1}=y{}_k+\frac{h}{2}[f(t{}_k,y{}_k)+f(t{}_{k+1},y{}_{k+1})],\\ k=0,1,2,\cdots ,n-1\end{array}(8)$

由梯形公式(8)可以看出,由k时刻的函数值yk也还不能直接计算出k+1时刻的函数值yk+1,因为在公式右端的第二项中依旧包含待求的函数值yk+1。因此,式(8)也是一个隐式公式,但它实际上已经是一个非线性方程,只要初值y0已知,即可根据上式,逐次递推地解出微分方程的数值解序列{yk(k=0,1,2,…,n)}。

由式(1)和式(8),得到Duffing系统的梯形公式如式(9)所示:

$\{\begin{array}{l}x(n+1)=x(n)+\omega (h/2)[y(n+1)+y(n)]\\ y(n+1)=(1-0.5\omega h)y(n)+\\ \omega h[x(n)-x{}^3(n)]+\gamma \omega h\text{c}\text{o}\text{s}[h(\omega n)]\end{array}(9)$

1.6 改进欧拉方程构建Duffing系统

针对Duffing振子的工程化模型建立展开研究,此模型建立的主要目的是为Duffing振子通信信号解调的FPGA实现创造前提和基础,用到的主要实验工具是DSP Builder软件[10]。

DSP Builder是嵌入在Matlab中的模块化处理软件,与Simulink模块软件不同,它可以通过相应的转换模块,转化成FPGA的编程VHDL语言,最终实现系统的硬件化集成。在运用DSP Builder的Duffing实现过程中,振子状态量x(n)只能够将x(n+1)通过延迟模块获得,因此为了在DSP Builder中搭建Duffing系统,一定要首先得到x(n+1)的表达式。因为在方程中包含x3(n)项和x(n)项,无法通过DSP Builder直接建立x(n+1)的确切表达式。为此,可将前向欧拉算法改进为式(10)形式,对式(4)中的第一个方程做延迟变形,得到改进后方程,如式(10)所示:

$x(n+2)=x(n+1)+\omega hy(n+1)(10)$

对式(10)进行变形,得到式(11):

$y(n+1)=\frac{x(n+2)-x(n+1)}{\omega h}(11)$

将式(11)代入式(4),得到式(12),

$\begin{array}{l}x(n+2)=x(n+1)+(\omega h-0.5\omega {}^{2}h{}^2)\\ y(n)+\omega {}^{2}h{}^{2}x(n)-\omega {}^{2}h{}^{2}x{}^3(n)+\\ \gamma \omega {}^{2}h{}^2\cos [h(\omega n)](12)\end{array}$

联立式(11)和式(12),得到改进的Duffing系统的欧拉方程,如式(13),

$\{\begin{array}{l}x(n+2)=x(n+1)+(\omega h-0.5\omega {}^{2}h{}^2)\\ y(n)+\omega {}^{2}h{}^{2}x(n)-\omega {}^{2}h{}^{2}x{}^3(n)+\\ \gamma \omega {}^{2}h{}^2\cos [h(\omega n)]\\ y(n+1)=\frac{x(n+2)-x(n+1)}{\omega h}\end{array}(13)$

2 仿真实验

为比较四种欧拉方法用于检测的性能,做如下的仿真实验:

设Duffing振子的周期策动力为正弦信号,并取系统内置周期策动力频率为50 Hz,内置策动力强度为γ=0.82,循环仿真步数设为4×104次,在未加噪声的前提下,借助Matlab仿真工具比较前向欧拉方程、后向欧拉方程、梯形公式和改进欧拉方程的Duffing振子状态相图,由图2～图5所示,由图2可以看出,当选用不同的模型建立Duffing系统时,使系统发生相变的内置策动力的临界值有显著的不同。

即γ=0.82时,在未加噪声的前提下,运用前向欧拉方程和改进欧拉方程,系统已经进入大尺度周期态,而后向欧拉方程和梯形公式产生的Duffing系统尚处于混沌状。

由此可见,运用不同的模型进行Duffing系统实现将会改变使Duffing振子发生相变的策动力临界值,但不会对系统用于检测的两个状态产生影响,即运用不同的欧拉方法构建Duffing振子模型,可以通过调整振子的临界策动力值,完成相应的检测过程。

3 结 语

(1) 运用不同的状态模型建立Duffing系统会改变振子相变的临界策动力幅值,但并不会对系统固有存在的状态产生影响,即通过改变理论和工程实验中Duffing振子的策动力临界值,即使模型不同,也并不影响信号检测的结果;

(2) 通过模块化设计DSP Builder的Duffing振子改进欧拉形式,可以实现从软件到硬件的过渡,为检测系统的工程实现奠定了基础。

参考文献

[1]王冠宇,陈大军,林建亚,等.Duffing振子微弱信号检测方法的统计特性研究[J].电子学报,1998,26(10):38-44.

[2]李士林,尹成群,尚秋峰,等.基于图像识别理论的混沌特性判别方法[J].中国电机工程学报,2003,23(10):47-50.

[3]赵华,尹成群,尚秋峰,等.双振子差分混沌特性判别方法[J].中国电机工程学报,2006,26(23):32-35.

[4]HU N Q,WEN X S.The application of Duffing oscillator in characteristic signal detection of early fault[J].Journal of Sound and Vibration,2003,68(5):917-931.

[5]田增锋,魏跃春,胡良剑.随机微分方程Euler法的均方稳定性和指数稳定性[J].自然杂志,2002,24(6):369-370.

[6]叶青华,黄海宁,张春华.用于微弱信号检测的随机共振系统设计[J].电子学报,2009,37(1):216-220.

[7]方同.工程随机振动[M].北京:国防工业出版社,1995.

[8]王永生,姜文志,赵建军,等.一种Duffing弱信号检测新方法及仿真研究[J].物理学报,2008,57(4):2053-2059.

[9]谢涛,魏学业.混沌振子在微弱信号检测中的可靠性研究[J].仪器仪表学报,2008,29(6):1265-1269.

[10]文忠,李立萍.基于Duffing振子的弱Chirp信号检测与参数估计[J].自动化学报,2007,33(5):236-539.

Duffing 篇2

强噪声背景下微弱周期信号的检测方法研究对于故障诊断、目标探测与识别、无线通信和生物医学信息提取等领域具有重要指导意义，同时也是现代信号处理领域的热点与难点之一。传统的功率谱分析、相干平均、线性滤波等方法在某些理想情况下可以检测出噪声掩盖下的微弱信号，但一般对背景噪声的统计特性有特殊要求，或因需要长时间的数据资料而效率较低。一般情况下，当带限噪声的频率范围覆盖信号频率且其能量远大于待测信号(信噪比低于-20dB)时，很难准确检测出信号是否存在[1]。利用混沌系统检测微弱周期信号是近年来提出的一种新型信号检测方法，具有常规方法所不具备的独特优势，可以利用较少的测量数据在任意色噪声背景下实现极低的信噪比工作门限[2,3]。近年来，在利用混沌系统的突变特性检测微弱谐和信号的研究中，国内学者取得了一定进展。最初是通过Duffing-Holmes方程的强迫激励分叉进行检测[4]，信噪比可以达到常规方法的下限。在此基础上提出的修正Duffing-Holmes方程系统[5]可进一步降低信噪比工作下限，增强检测效果。本文在对混沌弱信号检测法进行分析的基础之上，结合齿轮裂纹实验数据，研究了该方法在齿轮裂纹检测中应用，并对其结果与频谱分析法进行比较，通过比较可以看出，其检测性能得以大大提高。

1 修正Duffing弱信号检测原理

为了进一步提高Duffing振子检测系统对于微弱信号的敏感性，文献[6]提出了一种修正Duffing方程，可以获得更低的检测信噪比门限，尤其适合于检测微弱谐和信号，其模型为：

相比于Duffing-Holms方程，该模型对系统非线性项作了修改，并综合考虑了系统检测信噪比和系统混沌判据等方面的因素。以强迫激励频率ω为采样频率对系统变量x进行采样并记录采样值。修正Duffing方程出现混沌的参数阈值也可由Melnikov方法对其进行计算。

为了检测其它频率的周期信号，可令t=ωτ,上述系统可变为时间尺度τ上的动力学方程。写成状态方程形式为：

间歇混沌(又称为阵发混沌)是非线性系统在时间和空间上表现出的有序和无序交替出现的特殊动力学形态。在某些时空段，运动十分接近规则的周期运动;而在规则的运动段落之间，又夹杂着看起来很随机的跳跃。本文主要利用Duffing振子由于微小频差导致的间歇混沌现象，研究微弱信号检测方法。

带外加摄动力的修正Duffing振子表示为[7,8,9]：

式中，γrcosωt+dcos[(ω+△ω·ω)t+φ]为总驱动力;γrcosωt为内驱动信号;dcos[(ω+△ω·ω)t+φ]为外界小幅值的摄动信号，是对内驱动信号的摄动;△ω·ω为周期摄动信号与内驱动力间的频差，γr+d稍大于γc，γc为Duffing振子从混沌状态到大尺度周期状态产生相变的临界阈值，n(t)为噪声。

在上述系统中，总驱动力的幅值F(t)将在γr+d和γr-d之间周期性地交替变化。其中，γr与d之间的合成驱动力的幅值将在某一时段大于γc，系统因此过渡到大尺度周期运动状态;而在另一时段小于γc，系统则退化到以前的混沌状态。这样系统就出现时而混沌，时而周期的间歇混沌现象。

当△ω很小(△ω<0.03，△ω>0.04时将很难辨别出有规律的间歇混沌现象)时，F(t)的变化非常缓慢，远远慢于相变过程，因此系统对驱动力的缓变能够很好地响应，通过Duffing振子的时序图则可以看出周期和混沌的出现是泾渭分明的。作为本方法的关键之处，我们正是要通过Duffing振子时序图中这种泾渭分明的周期和混沌的交替出现来判断待检信号中是否包含所需的信号。

2 数值研究

以修正Duffing方程为基础，可以构造微弱信号检测系统如下：

其中输入信号u(t)=acos(ωt)+n(t),a表示待测谐和信号的幅值，n(t)表示背景噪声，γcos(ωt)是系统内置的激励信号。具体为：预先设定混沌系统周期策动力的频率为外界输入信号频率的下限值，并将杜芬振子的内策动力f调整为稍小于阈值fd，此时系统处于由混沌运动向周期状态转变的边沿。将带有强噪声的外界信号作为系统内部周期激励(内策动力信号)的补充引入系统，一旦带有与系统内策动力信号存在微小频差的周期信号时，尽管幅值很小，会导致杜芬振子出现混沌和周期交替出现的阵发混沌现象。同时根据阵发混沌的周期T的平均值(T1+T2+T3+…TN)/N,可求得外界输入弱信号频率为ω+2π/T,实验中将信号u(t)=acos(ωt)+n(t)并入系统。在系统(4)中，当a=0.01，白噪声标准差为σ=3,f=ω/2π分别取f1=19.04Hz、f2=20.52Hz、f3=170.03Hz、f4=179.38Hz时检测结果如图1所示:

从图1可以看出，本文分析的微弱信号检测方法是一种有效的检测方法。该方法对噪声免疫，经仿真分析可以求出，上述检测系统(4)的最低检测信噪比为SNR=-48dB，这相对于目前时域方法处理信号的最低信噪比门限SNR=-10dB来说有了很大的改进。

3 在齿轮裂纹检测中的应用实验

通过实验平台采集齿轮振动声信号(齿轮已经出现较大裂纹)，其信号波形如图2～3所示。利用频谱分析法进行检测分析，其结果如图4～5所示。

从图4和图5中可以看出，利用频谱分析法进行检测，由于噪声信号较强，在频率域内期望信号的幅值与噪声信号的幅值相比过小，几乎无法观测，从而导致期望信号不能正确地进行检测分析。

利用本文所提出的混沌弱信号检测法进行检测分析，其结果如图6～9所示。

从图6和图7中可以看出，信号作为外加激励加入到检测系统中，由于齿轮裂纹的声信号频率接近内驱动信号频率，所以尽管声信号幅值很小，并且噪声信号很强，但从系统的时序图中都可以观察到明显的间歇混沌现象发生。由于齿轮裂纹的声信号不是一个单频信号，而是一个窄带信号，所以观察到的间歇混沌现象不是规则的间歇混沌现象，但足以说明问题。图8和图9是齿轮无裂纹时采集的声信号送入检测系统后的检测结果，在本文中亦即纯噪声的检测结果。从这两幅图中可以看出，尽管噪声很强，但由于其频率与检测系统的内驱动信号频率相差较大，所以不会引起间歇混沌现象的发生。

4 结论

混沌系统对小信号的敏感性和对噪声的免疫力，使它在信号检测中非常具有潜力。基于此，本文详细分析了修正Duffing方程产生间歇混沌的原理，并且研究了利用间歇混沌原理检测强噪声中的微弱超声信号的方法。将方法与实验数据相结合，并与频谱分析法进行比较可得：混沌弱信号检测法原理简单、易于实现、具有很好的应用价值;与频谱分析法相比，具有较强的噪声免疫力和抗干扰能力。因此，混沌弱信号检测法是齿轮裂纹检测的一个有效检测方法，它为齿轮的安全工作提供了有效保障。

参考文献

[1]李月,杨宝俊.混沌振子检测引论[M].北京:电子工业出版社,2004.

[2]李月,杨宝俊.检测强噪声背景下周期信号的混沌系统[J].科学通报,2003,48(1):19~20.

[3]李月,路朋,杨宝俊等.用一类特定的双耦合Duffing振子系统检测强色噪声背景中的周期信号[J].物理学报,2006,55(4):1672~1677.

[4]王冠宇,陶国良,陈行等.混沌振子在强噪声背景信号检测中的应用[J].仪器仪表学报,1997,18(2):209-212.

[5]李月,路朋,杨宝俊等.用一类特定的双耦合Duffing振子系统检测强色噪声背景中的周期信号[J].物理学报,2006,55(4):1672~1677.

[6]李月,杨宝俊.混沌振子检测引论[M].北京:电子工业出版社,2004.

[7]Hu N Q,Wen X S.The application of Duffing oscillator in characteristic signal detection of early fault[J].Journal of Sound and Vibration,2003,268(5):917-931.

[8]Gan Chunbiao.Noise-induced chaos and basin erosion in softening Duffing oscillator[J].Chaos,Solitons and Fractals,2005,25(5):1069-1081.

Duffing 篇3

x″(t)+k·x′(t)-x(t)+x3(t)=f cos(t) (1)

f cos(t)为周期策动力,k为阻尼比,-x(t)+x3(t)为非线性恢复力。若k固定,则f超过一定阈值fc后继续增大,系统很快出现周期倍化分叉,随后进入混沌状态;随着f继续增大,当f大于另一个阈值fd时,系统进入大尺度的周期运动状态[1,2]。

文献[3]分析了如果要使系统运行轨迹发生改变,被测正弦信号的相位必须在一定范围内取值,而且相差过大时系统很容易收到噪声影响。而在实际的通信系统中,接受信号acos(ωt+θ)的相位不可能总是为0或接近0,而是随机分布的,因此本节讨论一种利用Duffing方程组检测随机相位的微弱信号的实现方法[4]。

1 原理分析

将2π均匀N等分,令N个Duffing方程的周期策动力分别预置这N个初相位[4],即Φi=(2πi/N),i=0,1,2,…,N-1。

得到如下N个方程构成的方程组

x″(t)+k·x′(t)-x(t)+x3(t)=f·cos(wt+Φi),i=0,1,2,…,N-1 (2)

仿照文献[5]中求混沌阈值的方法,容易推出无论Φi取何值,上述方程都存在混沌解,而且混沌阈值的表达式为[5]

当$\frac{f}{k}>0$时,解得$\frac{f}{k}＞\frac{4\cosh \left(\frac{\text{π}w}{2}\right)}{3\sqrt{2}\text{π}w}$,阈值为

$R(\omega )=\frac{4\cosh \left(\frac{\text{π}w}{2}\right)}{3\sqrt{2}\text{π}w}(3)$

当$\frac{f}{k}＜0$时,解得$\frac{f}{k}＜-\frac{4\cosh \left(\frac{\text{π}w}{2}\right)}{3\sqrt{2}\text{π}w}$,阈值为

$R(\omega )=-\frac{4\cosh \left(\frac{\text{π}w}{2}\right)}{3\sqrt{2}\text{π}w}(4)$

当k(k=0.5)固定时,各方程的阈值fd近似相等,但也有所不同。因此预先找出各fdi,并将各方程的策动力f 预设在fdi,此时各个方程的解皆处于混沌状态。

按照文献[3]中的参数设置,通过实验分别找出了16个阈值,如表1所示。

虽然已求出16个方程对应的阈值,但并不是16个方程同时用于检测信号。当待检测信号加入系统时,只有一个方程起作用,可以逐个验证每个方程,看是否有能够将信号检测出来的方程出现,将所有能检测信号的方程都作为有效,将检测到的信号求平均值。

待检测信号的未知相位θ可能会出现两种情况:

(1)待检测信号的相位θ与16个策动力的初相位θi,i=0,1,2,…,15之一相等。则当待测信号通过系统时,该θi对应的方程出现由混沌状态进入大尺度周期状态解的情况,可以准确地检测出待检测信号幅值。

(2)待测信号的相位θ处于两个θi之间。通过仿真实验发现,当待测微弱信号通过系统时,方程组中可能会有两个甚至多个出现临界状态的转变方程。为减小检测带来的误差,通过求平均的方法来减小误差。

2 仿真结果

上节求出16个方程对应的阈值,本节将用来检测信号。实验时选取k=0.5,ω=100π rad/s。

实验1 加入信号s=0.002 4cos(100πt+π)。

首先调整系统方程策动力的相位,使θ=θ0=0,然后调整幅值f0,使得f0=fd0=0.825 9,此时系统处于临界状态,然后将信号s加入系统,此时发现系统回到混沌状态,这说明待测信号的相位不在0附近。

然后继续改变系统策动力的幅值与相位,使θ=θ1= 0.125π,f1=fd1=0.825 9,此时系统又处于临界状态,然后继续将信号s加入系统,此时发现系统又回到混沌状态,说明待测信号的相位同样不在0.125π附近,然后继续将幅值与相位改变。

实验发现,当θ=θ8=π,f8=fd8=0.825 8时系统出现大尺度周期运动状态,利用信号检测方法调整f8,实验发现当f8=0.823 5时,系统重新回到临界状态。因此待测信号的幅值a=fd8-f8=0.825 8-0.823 5=0.002 3。误差为0.002 4-0.002 3=0.000 1。

实验继续下去,发现没有出现大尺度周期运动状态的方程,因此最终检测的信号幅度为0.002 3,相位为π。

实验2 加入信号s=0.002 4cos(100πt+15π/16)。

仿照实验一的步骤,当θ=θ7=0.875π,f7=fd7=0.826 0时系统出现大尺度周期运动状态,调整f7,实验发现当f7=0.824 4时,系统重新回到临界状态。因此待测信号的幅值a1=fd7-f7=0.826 0-0.824 4=0.001 6。误差为0.002 4-0.001 6=0.000 8。

实验继续,当θ=θ8=π,f8=fd8=0.825 8时系统出现大尺度周期运动状态,调整f8,实验发现当f8=0.824 3时,系统重新回到临界状态。因此待测信号的幅值a2=fd8-f8=0.825 8-0.824 3=0.001 5。误差为0.002 4-0.001 5=0.000 9。

实验继续下去,发现没有出现大尺度周期运动状态的方程。为了提高检测精度,采用时域经常用到的平均求值方法。将检测到的a1,a2求平均,可以得到

a=(a1+a2)/2=(0.001 6+0.001 5)/2=0.001 55 (5)

θ=(θ7+θ8)/2=(7π/8+π)/2=15π /16 (6)

3 结束语

对于相位未知的微弱正弦信号,可以采用16个方程的方法将带有未知相位微弱正弦信号检测出来,利用这种方法可以检测到的相位误差<π/16。上述仿真实验是在没有噪声输入时进行的,幅值误差主要为系统误差和相位差引起的误差。如果实际系统中有噪声存在,那么测量误差会增加,但是只要能检测出任意相位的输入信号,这种检测方法则是有效的。

参考文献

[1]王冠宇.混沌振子在强噪声背景信号检测中的应用[J].仪器仪表学报,1997,18(2):209-212.

[2]聂春燕.基于混沌相平面变化的弱信号检测方法研究[J].长春大学学报,1999,9(4):1-4.

[3]刘立,孙军.基于混沌振子的微弱信号检测方法研究[J].沈阳农业大学学报,2005,36(6):667-670.

[4]李健,何坤,乔强,等.应用混沌系统实现弱信号的检测[J].四川大学学报:自然科学版,2004,41(6):1180-1183.

Duffing 篇4

电力系统基波频率不仅是电网信号的重要指标参数,也是影响电力系统安全稳定运行的重要因素［１］。基波频率的波动会导致电网电能质量发生变化,从而影响电力系统的分析结果,特别是基波频率波动过大时,很可能影响某些电力系统的检测装置精度［２］,因而基波频率的精确检测对电力系统谐波分析与补偿具有重要意义。目前常用的频率测量方法主要有以下几种:1加窗相位差校正方法［３］,其实时性不够好;2跟踪频率采样算法,该算法虽然可以实现高精度频率测量［４］,但是采样频率也需跟踪调整,不便进行实际应用;3 强跟踪滤波器(STF)算法［５］,该算法与跟踪频率采样算法［４］相比具有更快的跟踪能力,但其采样频率要求高达10kHz,且涉及大量矩阵运算,特别是在频率突变的情况下,与鲁棒型扩展复数卡尔曼滤波(RECKF)算法［６］同样存在较大的跟踪超调;4双正交数字滤波器与加权平滑相位差分(DODF-WSPD)算法,以上2种算法对稳态信号具有较好的跟踪效果［７］,但对暂态信号的跟踪较慢,且存在较大的跟踪超调;5改进的插值离散傅里叶变换(DFT)算法,该算法的缺点在于计算量较大［８］,因而难以进行实际应用;6过零点测频算法［９］,该算法对谐波和白噪声干扰较为敏感;7其他方法,如最小均方误差法［１０］、牛顿递归法［１１］、自适应滤波法［１２］以及正交分量法［１３］等,这些方法在高谐波和噪声干扰下,通常不能同时兼顾静态输出精度和动态影响特性。此外,文献[14]提出了基于幅值线性变化模型的频率测量新算法;文献[15]提出了一种自适应基波提取与频率跟踪算法;文献[16]在传统傅里叶算法的基础上对其进行了改进,并提出了迭代算法与均值最小误差算法,取得了较好的效果。

本文提出了基于扩展型Duffing振子测量基波频率的方法,该方法有效避免了传统Duffing振子只能检测低频信号的局限性,具有良好的噪声免疫特性,仅用一组确定的参数条件即可对电力系统基波频率进行高精度检测。

1 扩展型Duffing振子模型与混沌检测特性

传统Duffing系统方程为［１７－１８］:

式中:k为阻尼系数;-x+x３为非线性恢复力;A为激励幅值;ω 为激励频率,其值为1rad/s;τ 为广义时间变量;F(ωτ)=Acosωτ为周期驱动力。

由于非线性项的存在,Duffing方程具有丰富的非线性动力学特性,是研究混沌的常用模型之一。混沌系统的特性之一是对初始状态的敏感性,初始状态的微小变化都将导致系统状态的显著变化。这一初始状态的敏感性等价于对系统参数A的敏感性［１７－１８］。

1.1 Duffing振子在微弱信号检测中的局限性

尽管基于传统Duffing振子的检测方法可以检测含有噪声的微弱周期或准周期信号,但该方法仍存在以下几方面的局限性。

1)系统临界混沌幅值的多变性

系统临界混沌幅值Aｄ的取值依赖于Duffing系统的各参数条件,即完全由阻尼系数、初始状态x(0)、初始状态的一阶导数x·(0)、周期驱动力角频率以及Runge-Kutta法的计算步长等因素共同决定。上述任何一个因素发生改变都会引发Aｄ的改变。

2)系统临界混沌幅值的不确定性

在保持其他参数不变的情况下,如果系统的激励频率不同,相应的Aｄ也不同。当待测信号的频率未知时,系统的周期驱动频率难以确定,因此Aｄ也就难以确定。

3)不良的动态响应品质与检测效果

Duffing系统仅在小频率(1rad/s左右)参数条件下具有良好的动态特性和检测效果,实际工程中则较难满足小频率成分的条件。通常情况下,实际信号可能为宽频带信号。在该条件下,系统的动态响应特性会恶化,将给信号检测带来较大的困难和误差。

4)噪声免疫功能失效

在保持其他参数不变的情况下,如果Runge-Kutta法的计算步长(即待测信号采样步长)改变,会引发Aｄ的改变,甚至会使噪声免疫特性明显变差,文献[17]对该观点进行了仿真验证。

基于上述局限性,本文提出了扩展型Duffing振子模型。

1.2 扩展型Duffing振子模型

对方程(1)进行广义时间尺度变换,设t=ωτ,则有= dx/dt = (1/ω)dx/dτ。 同理, =(1/ω２)d２x/d２τ,将τ换为t,代入方程(1)并整理,可得:

式中:ω=2πf,f为激励频率。

设x１=x,x２==１,代入方程(2),整理可得:

由此可见,式(1)经广义时间尺度变换后得式(3),变换后方程的性质和式(1)是完全一致的,但激励频率已由原先的1rad/s变为ω,表明经过广义时间尺度变换后,由原来激励角频率1rad/s的Duffing振子模型(1)变成了任意激励角频率为 ω的扩展型Duffing振子模型(3),因而获得了具有频率尺度变换的扩展型Duffing振子。由于变换不改变方程的原有性质,因此,在所有参数条件包括采样频率都相同的情况下,变换后的激励角频率为ω 的系统混沌阈值和原系统激励角频率为1rad/s的阈值应完全相同。

1.3 扩展型Duffing振子混沌特性

取k=0.5,初始状态x(0)=x·(0)=0,设采样频率fｓ=4 000Hz,则计算步长Tｓ=1/fｓ。令f=2fｅ=100Hz,其中fｅ为50Hz的标准激励频率,则可得系统的激励角频率为 ω =2πf =4πfｅ≈628rad/s。用4阶Runge-Kutta方法进行求解,计算点数为M=4 000,并由后2 000个点绘出系统的相轨迹。由仿真可计算确定扩展型Duffing系统的临界混沌幅值为Aｄ=0.825 9,大尺度周期状态临界幅值为Aｐ=0.826 0。临界混沌状态和大周期状态的相轨迹分别如图1和图2所示。

图1和图2 表明,扩展型Duffing振子在周期驱动力幅值A=Aｄ=0.825 9 时处于临界混沌状态,而在A=Aｐ=0.826 0时则处于大尺度周期运动状态。很显然,扩展型Duffing振子在临界混沌状态时,周期驱动力幅值只需发生0.000 1的微小变化就会使扩展型Duffing振子从临界混沌状态过渡到大尺度周期状态。周期驱动力幅值的这一敏感特性为微弱信号的检测奠定了理论基础。同时也表明,扩展型Duffing振子与传统Duffing振子相比具有对任意频率微弱信号的检测潜力。

1.4 扩展型Duffing振子噪声免疫特性测试

对式(3)取A=Ad=0.825 9,即扩展型Duffing振子处于临界混沌状态。如果此时只输入50dB的Gauss白噪声信号Nb(t),则噪声免疫特性测试系统将变为:

式中:F(t)=Aｄcos t;ω=2πf=4πfｅ。

在保证其他所有参数不变的情况下,Gauss白噪声输入到系统中后,扩展型Duffing振子仍然处于混沌状态,且具有良好的动态品质,如图3所示。

图3表明,扩展型Duffing振子具有较强的噪声免疫特性。

2 电力系统基波频率检测原理

设电力系统的待测基波信号为s(t),噪声信号为N(t),在系统(3)中取A=Aｄ,即扩展型Duffing振子处于临界混沌状态,将计入噪声的待测基波信号f(t)=s(t)+N(t)输入到系统(3),可得基于扩展型Duffing振子的基波信号检测系统,即

式中:F(t)=Aｄcos t;α为待测信号放大系数。

如果将周期信号视为二维直角坐标系上旋转的向量,则由于电力系统基波和各次谐波的初始相位角不为零,因此相同频率的系统驱动信号与电力系统待测信号在任何时刻都无法同步,通过对50 Hz标准激励信号频率的微调难以实现对电力系统基波信号的高精度检测。为此,本文将式(5)中驱动信号的激励频率修改为如下形式:

则系统的激励角频率为:

式中:n为满足1<n<4的实数。

经过式(6)的处理后,由于系统激励频率大于待测信号的基波频率,因此,在一周期(2π弧度)内,驱动信号与待测基波信号f(t)至少同相叠加1次,从而使得系统周期驱动力幅值增加。当且仅当驱动力幅值达到或超过大尺度周期运动状态的临界幅值时,系统(5)就会从临界混沌状态转化为大尺度周期状态,从而检测出f(t)的基波频率。

由于待测基波信号的频率未知,仅已知其在50Hz附近波动,因此,具体检测方法如下:在50Hz附近微调系统驱动信号的标准激励频率,使其与待测基波信号频率f１在误差范围内相近,即fｅ≈f１时,若取n=2.01,则系统激励频率为:

此时系统的激励角频率为:

因此,扩展型Duffing振子在t时刻的瞬时相位角为:

则系统周期驱动力信号为:

设待测电力系统基波信号为:

式中:ф１为待测基波信号的初始相位。

比较式(11)和式(12)可知,系统驱动力信号与待测信号的形式相同,且都包含相同的频率成分f１,在一周期内的某一时刻tｋ,一定存在如式(13)所示的待测基波信号初始相位瞬时匹配等式:

式中:m为正整数。

因此,在时刻tｋ,系统驱动力F(tｋ)与待测基波信号f(tｋ)同相叠加,从而使周期驱动力的幅值由原来的Aｄ增加到Aｄ+A１,当且仅当Aｄ+A１大于大尺度周期运动的幅值Aｐ时,扩展型Duffing振子由临界混沌状态转化为大尺度周期运动状态,从而检测出待测信号的基波频率f１。以上即为扩展型Duffing振子的基波频率检测原理。

3 仿真实验

目前智能变电站中最常用的采样频率为4kHz,国标GB/T 15945—2008规定的小系统允许的频率偏差范围为50±0.5 Hz。由于系统在不同激励频率下的临界混沌幅值略有偏差,而且待测基波信号的放大系数也略有偏差,因此,为了有效提高待测基波信号频率的检测精度,首先应求解Aｄ和α在标准测试信号情况下各频率点的具体数值。本文仿真实验中,对检测系统(5)的参数作如下设置:k=0.5,x(0)=x·(0)=0,fｓ=4 000 Hz,系统激励频率f=2.01f１。

标准测试信号为:

首先应确定系统激励频率下的临界混沌幅值Aｄ,然后在临界混沌状态时输入频率为f１的标准信号,通过调整α 使系统处于大尺度周期状态。所得的结果如表1所示。

3.1 仿真算例1

取以下待测基波信号［１６］:

仿真实验1 中,f１分别取49.6,49.7,49.8,49.9,50.0,50.1,50.2,50.3,50.4 Hz,噪声则取为50dB的Gauss白噪声信号nｂ(t),计算用数据长度为2 000点。临界混沌幅值和放大系数根据表1对应确定,由于式(15)的待测信号幅度为100,因此,表1的放大系数在仿真实验中应相应缩小100倍。所得的仿真结果如表2所示。

3.2 仿真算例2

设仿真使用的信号如式(16)所示［３］:

式中:A０为直流分量;Aｎ为n次谐波的幅值;фｎ为n次谐波的初始相位。

仿真实验2 中,f１分别取49.6,49.7,49.8,49.9,50.0,50.1,50.2,50.3,50.4 Hz。采样频率取4kHz,各次谐波的幅值和相位如文献[3]中的表1所示。由于文献[3]中谐波的最大幅值为100V,因此表1中α相应缩小100倍。相应的临界混沌幅值根据表1确定。取n=2.01,则扩展型Duffing振子的激励频率为f=2.01fｅ,计算用数据长度为2 000点,通过对49.6~50.4 Hz共9 个频率值进行频率扫描,并在相应频率值处对fｅ进行微调,可获得待测基波频率的近似值。此外,仿真实验中,将待测信号的直流以及各次谐波分量作为随机噪声看待,则可得待测基波频率的测试结果如表3所示。

4 结语

本文提出了一种基于扩展型Duffing振子的频率检测方法,该方法在非同步采样下具有较高的频率测量精度,仿真实验结果验证了该方法的有效性。不管是独立基波信号还是混有直流与各次谐波成分的信号,本文方法都能高精度地获取待测基波信号的频率,且频率检测的绝对误差小于0.001 Hz,比国标GB/T 14549—1993《电能质量公用电网谐波》中规定的误差范围小1个数量级。

Duffing 篇5

一、基于Duffing-Holms模型序参量的混沌模型

在经典力学当中, Duffing-Holms模型作为具有摆动的非线性方程, 同时还具有混沌现象的典型特征。一般认为Duffing-Holms模型的基本形式是:

其中, x表示系统的状态;ẍ表示系统变化的加速度;p1表示阻尼系数;q cos (ωt) 表示系统外力项;ω表示外力项频率;p1、p3和q一般大于0, p2<0 (p2取-1) 。为方便对本模型系统进行分析, 并实现降价的目的, 可以按照以下方式描述系统:

假设各参数取值分别为:p1=0.154, p2=-1, p3=4, ω=1.1, q=0.034, 同时系统初始值取x1 (0) =0, x2 (0) =ε (ε是一个计算机可识别的极小常数, 为10-10) , 由此可推断出系统处于周期振荡状态, 见图1。

在p1=0.154, p2=-1, p3=4, ω=1.1, q=0.088时, 系统由周期振荡状态转化为混沌状态, 此时其根轨迹和混沌相轨迹如图2 (a) 、 (b) 所示。

由此可知, 在频率ω不变的情况下, 模型系统的特性会随着系统外力q cos (ωt) 的外力系数q的改变而发生明显变化。当q取值0.088时, 系统处于临界状态, 此时系统由周期振荡模型转为混沌模型。

二、基于Duffing-Holms模型的混沌控制方式

(一) 变量反馈控制方式

在外力参数q取值0.088时, 系统表现为混沌模型。此时可通过引入合适的反馈变量来实现对系统混沌的控制。例如在 (1) 式中引入反馈项-kx1, 那么可得到:

设置k=0.06, 且其他参数保持不变, 此时可实现对混沌模型的控制, 并转化为周期振荡状态。见图3。

由上述分析可知, 通过引入合适的反馈变量-kx1, 并对其反馈系数k进行调节, 实现对混沌模型的特性的控制, 将其转化为周期振荡状态。

(二) 参数微扰控制方式

基于Duffing-Holms模型序参量的混沌模型, 除了可以通过变量反馈法加以控制, 以及自适应控制方式、外力反馈控制方式和延时反馈控制方式等, 对系统参数随时间的连续微扰进行控制之外, 还可通过参数共振微扰的方式来进行控制。

考虑到混沌模型受参数变化的影响较为剧烈, 以及Duffing-Holms模型当中的x13的倍增性质, 本研究将参数p3设置为p3 (1+c cosΩt) , 其中, c表示参数微扰幅度;Ω表示微扰频率。在微扰频率Ω和外力项频率ω产生共振的过程中, 基于Duffing-Holms模型序参量的混沌模型的相关特性会被抑制, 同时也可让该模型回到期望的周期振荡状态。此时 (1) 式可重新表达如下:

三、结束语

本研究通过利用变量反馈方法和参数微扰方法对基于Duffing-Holms模型序参量的混沌模型进行控制, 可以看出, 变量反馈方法只需要引入并调节反馈系数k, 便可实现对混沌模型的特性的控制, 甚至将其转化为周期振荡状态;而参数微扰方法则可通过外力项频率和微扰频率的共振, 来改善混沌模型, 并使其能够达到所期望的特性。在实际生活当中, 可将该模型广泛应用于金融市场, 并通过变量反馈方法和参数微扰方法来实现对系统的有效控制, 以确保系统的稳定和可靠。

参考文献

[1]马超群, 邹琳, 李红权.股票市场的非线性结构与混沌效应检验:基于BDS方法与CR方法[J].湖南大学学报, 2008 (5)