CN101561314A - 随机共振-混沌微弱信号检测方法 - Google Patents

Abstract

本发明公开了一种随机共振－混沌微弱信号检测方法，首先通过一种自适应参数调节随机共振方法实现微弱特征频率信号是否存在的检测，得到实际待检信号最可能的信噪比大小；在此基础上，设计出相应的Duffing混沌振子，并计算得出不同信噪比条件下混沌振子的输入周期激励与输出大尺度周期运动幅值间的变化关系，再根据随机共振方法估计得出的待检信号信噪比大小即可估计得出微弱特征频率信号的幅值大小。该方法适合于物理、化学、生物、机械故障早期检测等领域低信噪比条件下微弱特征频率信号的检测。

Description

随机共振-混沌微弱信号检测方法

技术领域

本发明涉及一种强噪声背景下的微弱信号检测方法，特别适用于低信噪比条件下，机械系统早期故障诊断中对相关故障是否发生以及发生程度的全面检测。

背景技术

随机共振技术于上世纪八十年代由意大利学者巴兹等人在研究地球古气候变化时提出。随机共振利用噪声增强微弱信号传输的机制，使其与其它微弱信号检测方法相比具有独特的优势，因而在生物信号处理、视觉图像与听觉识别、电磁系统及光信号处理等领域得到广泛重视。研究表明，利用随机共振方法能把微弱特征频率信号是否存在有效检测出来，但由于随机共振发生机制中包含噪声的能量向信号能量的转移，且系统产生随机共振时噪声不是某个确定的值，而是某个区间，故无法直接用随机共振方法进行微弱信号幅值估计。机器实际运行过程中，需要全面掌握故障是否发生以及发生的程度，以便进行合理维修。而传统的微弱信号幅值估计方法(如最大似然法等)，在信噪比较低时，估计效果不理想。由此，关于随机共振在机械系统早期故障微弱特征频率信号检测中的应用存在以下亟待解决的问题：低信噪比条件下，如何在利用随机共振方法实现微弱特征频率信号是否存在的基础上，对其幅值进行估计？

发明内容

本发明目的在于解决随机共振技术在实际微弱特征频率信号检测中无法估计信号幅值大小的问题，提供一种低信噪比条件下，既能快速、准确判断出特征频率信号有无，又能对其幅值大小进行估计的微弱信号检测方法。

为了实现上述发明目的，本发明的随机共振-混沌微弱信号检测方法包括如下步骤：首先通过一种自适应参数调节随机共振方法实现微弱特征频率信号是否存在的快速、准确检测，得到实际待检信号最可能的信噪比大小；在此基础上，设计出相应的Duffing混沌振子，并计算得出不同信噪比条件下混沌振子的输入周期激励与输出大尺度周期运动幅值间的变化关系，再根据随机共振方法估计得出的待检信号信噪比大小即可估计得出微弱特征频率信号的幅值大小。从而实现微弱特征频率信号的“有无”和“大小”的全面检测。

根据上述方法，本发明进一步的技术方案是采用如下步骤：

1)根据随机共振理论，对周期信号与噪声共同作用的双稳系统模型：

$\stackrel{\·}{x}=\mathrm{ax}-{\mathrm{bx}}^3+u\left(t\right)+\mathrm{\Γ}\left(t\right)---\left(1\right)$

<Γ(t)>＝0，<Γ(t)，Γ(t′)>＝2Dδ(t-t′)

进行数值仿真，得出模型(1)中a＝b＝1、周期信号u(t)频率f₀＝0.1Hz时，非线性双稳系统最佳匹配随机共振的噪声强度D₀(即噪声均方根值为σ₀， ${\mathrm{\&sigma;}}_0=\sqrt{2D_0}$ )与输入信噪比 ${\mathrm{SNR}}_{\mathrm{in}}({\mathrm{SNR}}_{\mathrm{in}}=20\log \left(\frac{A}{\mathrm{\&sigma;}}\right))$ 间的对应关系；

(1)式中x为双稳系统输出，a、b为系统结构参数，t是时间变量，Γ(t)是均值为0、噪声强度为D(即噪声均方根值为σ， $\mathrm{\&sigma;}=\sqrt{2D}$ )的高斯分布白噪声，u(t)是幅值为A、频率为f₁的微弱周期信号；

2)滤掉待检信号中的工频成分，估计其噪声均方根值σ₁；

3)改变待检信号信噪比值，同时对双稳系统模型(1)进行数值仿真，得到实际待检信号最可能的信噪比大小；根据步骤1)所得结果，即可得到该信噪比条件下，频率f₀＝0.1Hz所对应的最佳匹配随机共振噪声均方根值σ₀，再按照归一化尺度反变换

$a=\frac{f_1}{f_0},$ $b=\frac{{{\mathrm{\&sigma;}}_0}^2{f_1}^2}{{{\mathrm{\&sigma;}}_1}^2{f_0}^2}---\left(2\right)$

计算出模型(1)中的a，b值；

4)将a，b值代入模型(1)，得到对应于该输入信号的随机共振检测模型，根据该模型输出频谱图即可判断输入信号中是否存在特征频率为f₁(角频率为ω₁)的信号；

5)如果待检信号中存在特征频率为f₁的信号，即可针对f₁(即ω₁)，设计出相应的Duffing混沌振子

$\frac{1}{{{\mathrm{\&omega;}}_1}^2}\frac{{\mathrm{dx}}^2}{{\mathrm{dt}}^2}+\frac{\mathrm{\&delta;}}{{\mathrm{\&omega;}}_1}\frac{\mathrm{dx}}{\mathrm{dt}}-x+x^3={\mathrm{\&gamma;}}_c\cos \left({\mathrm{\&omega;}}_{1}t\right)+u\left(t\right)+\mathrm{\&Gamma;}\left(t\right)---\left(3\right)$

(3)式中，δ为阻尼系数，γ_c为根据Melnikov方法计算得出的混沌振子混沌解与周期解的分叉阈值，u(t)和Γ(t)在实际信号检测中往往混和在一起，构成实际待检信号；

6)采用四阶Rung-Kutta算法对混沌振子(3)进行数值求解，在此基础上估计其大尺度周期解的幅值A₁；

7)通过对混沌振子(3)进行数值仿真，得到待检信号信噪比在某区间内变化时，混沌振子(3)的激励幅值大小γ与相应的振子大尺度周期解幅值A₁间的变化关系为：

A₁(γ)＝p₁γ+p₂                (4)

(4)式中，γ为混沌振子(3)中分叉阈值γ_c与周期信号u(t)幅值A之和，A₁由步骤6)得出，p₁和p₂分别为线性拟合系数；

8)根据步骤3)所估计的待检信号信噪比值，结合步骤7)所得到的该信噪比条件下的振子大尺度周期解幅值A₁和振子激励幅值γ间的拟合经验公式(4)，即可求出γ值；

9)用γ减去混沌振子(3)中的分叉阈值γ_c，即可估计出所检测的微弱特征频率信号幅值。

本发明是一种非线性信号检测方法，其特点在于将随机共振与混沌两种方法巧妙地结合起来：所设计的参数调节随机共振方法(前4个步骤)在快速、准确检测出微弱特征频率信号是否存在的同时，还为基于混沌振子的微弱信号幅值估计方法(后5个步骤)提供所需的待检信号信噪比大小，再设计相应的Duffing混沌振子估计出微弱信号幅值大小，估计效果优于传统的最大似然法。本发明将以上微弱信号检测方法命名为随机共振-混沌方法。该方法回避了混沌振子微弱信号检测中系统是否发生相变的判断这一难点，弥补了随机共振无法检测信号大小的不足，为微弱信号检测提供了一条新的途径。该方法适合于物理、化学、生物、机械故障早期检测等领域低信噪比条件下微弱特征频率信号的检测。

下面结合附图进一步说明本发明的技术方案。

附图说明

图1是原始碰摩信号时域波形。

图2是原始碰摩信号频谱图。

图3是处理后的碰摩信号时域波形。

图4是处理后的碰摩信号频谱图。

图5是σ/A＝10时，信噪比增益与信号幅值关系曲线。

图6是自适应参数调节随机共振模型输出时域波形。

图7是自适应参数调节随机共振模型输出局部频谱图。

具体实施方式

下面用一个具体实例来说明本发明用于微弱信号检测的具体过程与有益效果：

动静件碰摩是转子系统中较常见且难以检测和捕捉的一类故障现象。主要发生在机组动静叶片密封、转子轴封以及滑动轴承等处。碰摩发生的时间短且有时位置不确定，如果碰摩现象经常发生，可能引发其它严重故障的发生，造成重大经济损失。因此，设法在转子系统碰摩故障发生早期将其检测出来，对于转子系统避免碰摩及继发性故障以及机组的健康运行具有现实意义。

通过试验研究转子系统中动静件间尖锐碰摩时的振动特征规律可知在早期碰摩阶段，发现有工频的1/3X、2/3X等分量稳定存在，这一特性可以为这类故障的早期诊断提供依据。下面据此特性，采用随机共振-混沌方法，对碰摩故障进行早期检测。

图1是工频为60Hz时的一组典型数据，采样频率为f_s＝1000Hz，数据长度为N＝2000。上述获得的试验数据是在实验理想环境中获得的，表征早期碰摩故障的各种微弱特征还能在谱图中隐约可见，如果在现场的转子系统运行环境中，这样微弱的特征信号将被背景噪声所完全淹没。为分析方便，假设背景噪声为高斯白噪声，将图1中数据进行滤波处理(此操作同时完成了步骤2)中的滤除工频信号影响)，仅考虑包含被检测频率的窄带宽的情况，所采用的带通数字滤波器(Chebyshev I型6阶滤波器)的通带截止频率分别为16Hz、30Hz，阻带截止频率分别为5Hz、40Hz，通带波纹0.5dB，阻带衰减-45dB。再加入同样长度且强度为D＝2.5的高斯白噪声数据，结果如图3所示。图2、图4分别为上述两种信号的功率谱图。由图1可知，图1所包含的数据存在频率为20Hz、40Hz两个明显的分量。由图3、图4根本分辨不出该信号中包含有弱周期信号。首先采用参数调节随机共振方法对其特征频率信号进行检测。

假设针对20Hz信号设计出的双稳系统随机共振检测模型为：

$\stackrel{\&CenterDot;}{x}=\mathrm{ax}-{\mathrm{bx}}^3+s\left(t\right)---\left(5\right)$

其中，s(t)为图3描述的待检信号，相当于模型(1)中的周期信号u(t)和噪声信号Γ(t)之和。

步骤1)：通过对模型(1)进行数值仿真，得到a＝b＝1、周期信号u(t)频率f₀＝0.1Hz时，非线性双稳系统最佳匹配随机共振的噪声强度D₀(即噪声均方根值为σ₀)与输入信噪比间的对应关系。

模型(1)中a＝b＝1、f₀＝0.1Hz时，对模型(1)进行数值仿真，采用四阶Rung-Kutta(参见“黄明游编.数值计算方法.北京：科学出版社，2003”)算法进行数值求解，取计算步长h＝0.05，数据长度N＝2000。实际操作中采用噪声均方根值与信号幅值的比值σ/A来描述输入信噪比SNR_in。下面以σ/A＝10时双稳系统最佳匹配随机共振噪声(以下简称最佳噪声)的计算为例，来说明不同σ/A时最佳噪声的计算过程。

令模型(1)中周期信号u(t)幅值A以0.05的步长在[0.12.5]区间内变化，保持σ/A＝10不变，对模型(1)进行数值仿真，容易得到信噪比增益(输出信噪比与输入信噪比的比值)与信号幅值之间的变化曲线如图5。由图可知σ/A＝10时，双稳系统最佳匹配随机共振(图5中信噪比增益最大值处)时周期信号幅值A₀＝0.4，即最佳噪声均方根值σ₀＝4。

为进一步获得各个信噪比条件下的最佳噪声，让输入信噪比σ/A值在5到20之间(大量实验结果表明，随机共振方法在该区间内检测效果较好)，以0.5为步长递增，按照上述σ/A＝10时的计算方法，依次计算得出各个信噪比下的最佳噪声，计算结果见表1。

表1最佳匹配随机共振噪声大小与输入信噪比间的对应关系

步骤2)：滤掉待检信号中的工频成分，估计其噪声均方根值σ₁。

首先对待检信号进行滤波去除工频信号的影响(此操作在对信号进行预处理时已经完成)，再将待检信号近似为纯噪声估计出其噪声均方根σ₁＝0.2505。

步骤3)：计算实际待检信号最可能的输入信噪比，根据归一化尺度反变换关系式(2)计算实际检测模型(5)中的a，b值。

步骤2)已经估计得出了待检信号的噪声均方根值σ₁，假设待检信号的输入信噪比按照表1所列各值变化，并对模型(5)进行数值仿真，依次计算各个输入信噪比条件下的信噪比增益，找出信噪比增益最大值所对应的输入信噪比大小，即为该实际待检信号最可能的信噪比。

针对实际转子碰摩信号，计算得出该信号最可能的信噪比为σ/A＝10。由表1，对应模型(1)中a＝b＝1时，双稳系统的最佳噪声σ₀＝4.0。至此，式(2)中的各个参数值均以得到：f₁＝20Hz、f₀＝0.1Hz、σ₁＝0.2505、σ₀＝4.0，代入式(2)，即可计算得出模型(5)中的参数a＝200，b＝2.048×10⁷。

步骤4)：将a，b值代入实际检测模型(5)，得到对应于该输入信号的随机共振检测模型，根据该模型输出频谱图即可判断输入信号中是否存在特征频率为20Hz的信号。

将a＝200，b＝2.048×10⁷代入模型(5)，采用四阶Rung-Kutta算法对其进行数值计算，计算步长h＝0.001、数据长度为N＝2000。图6、图7是其输出的时域和频域谱图，从图7中可以清晰分辨出20Hz信号，即可确定待检信号中1/3分频信号的存在。

在此基础上，进一步采用混沌方法对其故障特征频率信号幅值大小进行估计，以便对其故障发生程度进行分析。

步骤5)：针对20Hz信号，设计出相应Duffing混沌振子为：

$\frac{1}{{{\mathrm{\&omega;}}_1}^2}\frac{{\mathrm{dx}}^2}{{\mathrm{dt}}^2}+\frac{\mathrm{\&delta;}}{{\mathrm{\&omega;}}_1}\frac{\mathrm{dx}}{\mathrm{dt}}-x+x^3={\mathrm{\&gamma;}}_c\cos \left({\mathrm{\&omega;}}_{1}t\right)+s\left(t\right)---\left(6\right)$

其中，δ＝0.5，ω₁＝40π，γ_c＝0.8260，γ_c通过Melnikov方法(参见“刘曾荣编.混沌的微扰判据.上海：上海科技教育出版社，1994”)计算得到，s(t)为待检信号，相当于模型(3)中的周期信号u(t)和噪声信号Γ(t)之和。

步骤6)：应用四阶Rung-Kutta算法对混沌振子(6)进行数值求解，在此基础上估计其大尺度周期解的幅值A₁；

采用四阶Rung-Kutta算法对混沌振子(6)进行数值计算，步长h＝0.005，计算长度N＝2000，得到其大周期解序列，在此基础上，根据最大似然法估计出其大周期解幅值A₁＝1.4202。

步骤7)：数值计算得出待检信号信噪比在某区间内变化时，混沌振子(3)的激励幅值大小γ(γ_c与γ₁之和)与相应的振子大尺度周期解幅值A₁间的变化关系。

讨论ω＝1时的Duffing振子为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}+\mathrm{\&delta;}\stackrel{\&CenterDot;}{x}-x+x^3={\mathrm{\&gamma;}}_c\cos t+u\left(t\right)+\mathrm{\&Gamma;}\left(t\right)---\left(7\right)$

其中，δ＝0.5；u(t)为待检微弱特征频率信号，且u(t)＝γ₁cost；γ_c为根据Melnikov方法计算得出的混沌振子混沌解与周期解的分叉阈值，且γ_c＝0.8260。由基于混沌的微弱信号检测相关理论可知，当待检信号中存在ω＝1微弱周期信号时，振子处于大尺度周期运动状态。

设已经处于大尺度周期状态的Duffing振子周期运动振幅大小为A₁，振子激励的振幅大小为γ(γ_c与γ₁之和)。大量仿真计算结果表明，当γ₁以0.002为步长在区间[0.002，0.11]内取值，即γ在区间[0.8280，0.9360]内变化，待检信号信噪比σ/γ₁值在区间[5，20]内变化时，A₁与γ基本保持正比关系。用

A₁(γ)＝P₁γ+p₂                    (8)

表示各个σ/γ₁时拟合的A₁随γ变化的经验公式，表2列出了σ/γ₁以0.5为步长在区间[5，20]内变化时的各个p₁，p₂值。仿真实验结果表明，以下表格对混沌振子(3)中的频率ω≠1时都能成立。

表2不同信噪比对应的A₁与γ间经验公式拟合系数

σ/γ1	5	5.5	6	6.5	7	7.5	8	8.5	9
										p1	1.235	1.253	1.269	1.286	1.302	1.317	1.332	1.347	1.361
p2	0.3187	0.3049	0.2915	0.2786	0.266	0.2538	0.2419	0.2305	0.2193
										σ/γ1	9.5	10	10.5	11	11.5	12	12.5	13	13.5
p1	1.375	1.388	1.401	1.414	1.426	1.438	1.45	1.461	1.472
										p2	0.2085	0.198	0.1879	0.178	0.1684	0.1591	0.1501	0.1414	0.1329
σ/γ1	14	14.5	15	15.5	16	16.5	17	17.5	18
										p1	1.483	1.493	1.503	1.513	1.523	1.532	1.541	1.55	1.559
p2	0.1246	0.1166	0.1088	0.1013	0.094	0.0869	0.08	0.0733	0.0668
										σ/γ1	18.5	19	19.5	20
p1	1.567	1.575	1.583	1.591
										p2	0.0605	0.0544	0.0484	0.0427

步骤8)：计算混沌振子(7)中频率为20Hz周期信号的激励幅值大小γ。

根据步骤3)所估计的信号输入信噪比σ/A＝10，步骤6)计算得出的A₁＝1.4202，以及表2所列出的σ/A＝10时，p₁＝1.388，p₂＝0.198，将以上数值代入公式(8)，即可计算得出混沌振子(7)中频率为20Hz的激励幅值大小γ＝0.88055。

步骤9)：估计所检测的微弱特征频率信号幅值。

由于混沌振子(7)中的激励幅值γ相当于分叉阈值γ_c和实际待检信号微弱周期信号幅值

之和，故碰摩信号中20Hz成分幅值大小 ${\widehat{\mathrm{\&gamma;}}}_{20}=\mathrm{\&gamma;}-{\mathrm{\&gamma;}}_c=0.05455,$ 这与其实际大小γ₂₀＝0.05538非常接近。

按照同样方法可以实现故障信号中40Hz和80Hz频率成分的检测，在应用随机共振方法检测出存在基础上，估计出幅值大小分别为： ${\widehat{\mathrm{\&gamma;}}}_{40}=0.0636,$ ${\widehat{\mathrm{\&gamma;}}}_{80}=0.093,$ 实际幅值为γ₄₀＝0.0645，γ₈₀＝0.0958。可见，应用随机共振-混沌方法在成功检测出碰摩故障特征信号的同时，还实现了其特征频率信号幅值的估计，能够为实际故障检测中故障发生程度分析提供一定依据。

Claims (2)

1、一种随机共振-混沌微弱信号检测方法，包括如下步骤：

首先通过一种自适应参数调节随机共振方法实现微弱特征频率信号是否存在的检测，得到实际待检信号最可能的信噪比大小；在此基础上，设计出相应的Duffing混沌振子，并计算得出不同信噪比条件下混沌振子的输入周期激励与输出大尺度周期运动幅值间的变化关系，再根据随机共振方法估计得出的待检信号信噪比大小即可估计得出微弱特征频率信号的幅值大小。

2、根据权利要求1所述的随机共振-混沌微弱信号检测方法，其特征是包括以下步骤；

1)根据随机共振理论，对周期信号与噪声共同作用的双稳系统模型：

$\stackrel{\&CenterDot;}{x}=\mathrm{ax}-{\mathrm{bx}}^3+u\left(t\right)+\mathrm{\&Gamma;}\left(t\right)$

                                                (1)

＜Γ(t)＞＝0，＜Γ(t)，Γ(t′)＞＝2Dδ(t-t′)

进行数值仿真，得出模型(1)中a＝b＝1、周期信号u(t)频率f₀＝0.1Hz时，非线性双稳系统最佳匹配随机共振的噪声强度D₀(即噪声均方根值为σ₀， ${\mathrm{\&sigma;}}_0=\sqrt{2D_0}$ )与输入信噪比SNR_in( ${\mathrm{SNR}}_{\mathrm{in}}=20\log \left(\frac{A}{\mathrm{\&sigma;}}\right)$ )间的对应关系；

(1)式中x为双稳系统输出，a、b为系统结构参数，t是时间变量，Γ(t)是均值为0、噪声强度为D(即噪声均方根值为σ， $\mathrm{\&sigma;}=\sqrt{2D}$ )的高斯分布白噪声，u(t)是幅值为A、频率为f₁的微弱周期信号；

2)滤掉待检信号中的工频成分，估计其噪声均方根值σ₁；

3)改变待检信号信噪比值，同时对双稳系统模型(1)进行数值仿真，得到实际待检信号最可能的信噪比大小；根据步骤1)所得结果，即可得到该信噪比条件下，频率f₀＝0.1Hz所对应的最佳匹配随机共振噪声均方根值σ₀，再按照归一化尺度反变换

$a=\frac{f_1}{f_0},$ $b=\frac{{{\mathrm{\&sigma;}}_0}^2{f_1}^2}{{{\mathrm{\&sigma;}}_1}^2{f_0}^2}---\left(2\right)$

计算出模型(1)中的a，b值；

4)将a，b值代入模型(1)，得到对应于该输入信号的随机共振检测模型，根据该模型输出频谱图即可判断输入信号中是否存在特征频率为f₁(角频率为ω₁)的信号；

5)如果待检信号中存在特征频率为f₁的信号，即可针对f₁(即ω₁)，设计出相应的Duffing混沌振子

$\frac{1}{{{\mathrm{\&omega;}}_1}^2}\frac{{\mathrm{dx}}^2}{{\mathrm{dt}}^2}+\frac{\mathrm{\&delta;}}{{\mathrm{\&omega;}}_1}\frac{\mathrm{dx}}{\mathrm{dt}}-x+x^3={\mathrm{\&gamma;}}_c\cos \left({\mathrm{\&omega;}}_{1}t\right)+u\left(t\right)+\mathrm{\&Gamma;}\left(t\right)---\left(3\right)$

(3)式中，δ为阻尼系数，γ_c为根据Melnikov方法计算得出的混沌振子混沌解与周期解的分叉阈值，u(t)和Γ(t)在实际信号检测中往往混和在一起，构成实际待检信号；

6)采用四阶Rung-Kutta算法对混沌振子(3)进行数值求解，在此基础上估计其大尺度周期解的幅值A₁；

7)通过对混沌振子(3)进行数值仿真，得到待检信号信噪比在某区间内变化时，混沌振子(3)的激励幅值大小γ与相应的振子大尺度周期解幅值A₁间的变化关系为：

A₁(γ)＝p₁γ+p₂                    (4)

(4)式中，γ为混沌振子(3)中分叉阈值γ_c与周期信号u(t)幅值A之和，A₁由步骤6)得出，p₁和p₂分别为线性拟合系数；

8)根据步骤3)所估计的待检信号信噪比值，结合步骤7)所得到的该信噪比条件下的振子大尺度周期解幅值A₁和振子激励幅值γ间的拟合经验公式(4)，即可求出γ值；

9)用γ减去混沌振子(3)中的分叉阈值γ_c，即可估计出所检测的微弱特征频率信号幅值。

Images