指数型随机共振微弱振动信号检测方法

2019-05-27贺利芳

振动与冲击 2019年9期

张刚,曹莉,贺利芳,易甜

(1.重庆邮电大学通信学院，重庆 400065;2.信号与信息处理重庆市重点实验室，重庆 400065)

微弱信号指自身强度很小并且淹没在噪声中一类信号，微弱信号检测即为采用各种技术手段提高信噪比将淹没在强噪声中的微弱信号提取出来。一种技术手段针对噪声，通过噪声与信号之间差异抑制噪声，例如滤波技术、时频分析技术、相关检测技术，其应用涉及到实际生产的各个方面，如故障信号、地震信号、生物电信号，这些技术已经成为广大学者的研究热点[1]。另一种技术为增强信号，将噪声与信号通过非线性方程进行能量转化，将噪声能量转化为信号能量使得信号突出被检测。随机共振(Stochastic Resonance,SR)最早为1981年物理学家Benzi等[2]在对地球冰川期气候变化的研究中首次提出；1983年Fauve等[3]在施密特触发器实验中观察到了随机共振现象。随着各领域间随机共振发现，人们对随机共振现象机理研究越来越火热，并取得丰富研究成果[4-5]。

随机共振微弱信号检测技术在微弱信号检测及提取方面应用广泛，相比较于传统的抑制噪声技术，在保留信号信息的基础上，最大化将噪声能量转化为信号能量。常采用一维Langevin方程模型即传统双稳[6](Classic Bistable Stochastic Resonance,CBSR)系统，在参数满足绝热近似理论条件下实现信号检测，为了更好地利用噪声，一些新的模型陆续被提出并被研究，Qiao等[7]通过构造分段双稳从势函数势阱峭度研究势阱峭度对随机共振输出信号饱和度影响；周玉飞等[8]使用级联双稳研究能量多次转换问题；陆思良等[9]采用两个一维随机共振组合而成的二维互补随机共振应用于轴承信号检测中，这些模型大多基于一维非线性方程。随着研究的深入发现二阶非线性方程也能产生随机共振现象。Duffing振子方程是一类典型的二阶非线性方程不仅能通过混沌理论检测微弱信号，也能够产生随机共振现象。Duffing振子随机共振最早提出在20世纪80年代[10-11]，随后Almog等[12-13]搭建硬件平台实现Duffing振子随机共振信号检测，Leng等[14]研究不同参数调节Duffing振子随机共振对输出信号影响；Lai等[15]研究在Duffing方程下噪声诱导随机共振发生；Li[16]将Duffing方程随机共振应用在MSK信号中；Zhang等[17]提出的耦合双稳Duffing随机共振系统不仅能提高输出信噪比，有效检测正弦信号及方波信号。

Duffing振子随机共振扩展了随机共振理论与实践研究，本文提出指数型双稳作用在Duffing方程随机共振系统，推导其逃逸率并研究阻尼参数对随机共振输出信噪比影响；采用衰减冲击信号和振动谐波信号进行验证，证明所提系统有增强信号作用；并且提出随机共振与经验模态分解结合下信号检测方法，应用于轴承故障信号检测，研究结果表明有效的检测到故障频率以及故障的二倍频信号、三倍频信号。

1 势函数模型

1.1 简谐势阱模型

简谐势阱为线性随机共振系统原始模型，后来演化为幂函数型势阱。系统结构简单参数单一，其表达式[18]为

(1)

式中，a为系统参数，不同系统参数a下势函数图形如图1所示。a越大势阱越陡峭、势阱底部中心宽度越窄粒子运动穿过势阱底部时间越短。调节系统参数、信号参数、噪声使三者协同产生随机共振。简谐势阱只有一个系统参数不用考虑系统参数间的相互影响，使之成为研究热点之一。

图1 简谐势函数Fig.1 Harmonic potential function

1.2 Gaussian potential势阱模型

在核物理学中，常用GP模型来描述复核散射，径向GP模型[19]可以表示为

(2)

式中，V表示势阱深度，R表示势阱宽度。图2表示GP势阱的势函数，由图可知，势阱的两端迅速收敛于0。由图2(a)可知，固定V，单独调节R可以改变势阱的宽度，并且随着R的减小，势阱宽度逐渐减小，势阱壁逐渐变得陡峭。图2(b)显示，GP势的势阱深度由V唯一表征。

(a)V=1

(b)R=1

1.3 组合模型

基于以上两种模型，将简谐势阱模型和GP势阱模型相结合，提出一种指数型双势阱模型，势函数如下

(3)

图3为结构参数为：a=1，V=1，R=0.5的指数型双稳势函数。势阱函数的两个势阱是对称的，每个势阱的宽度和势阱壁的陡峭度均可通过参数a，V，R进行调节。在简谐势阱模型中加入GP势阱使其变成双势阱，使振荡粒子由在单一势阱中运动变成在两势阱间进行跃迁，提高噪声的利用率从而改善输出信噪比，得到较好的随机共振效果。从该新型系统的势函数来看，通过调节系统参数，可以实现模型在单稳态势阱和双稳态势阱之间的切换。因此，该新型模型兼具单稳势阱模型与双稳势阱模型的优点。

图3 指数型双稳势函数Fig.3 The exponential bistable potential function

2 随机共振

2.1 随机共振机理

随机共振模型为

(4)

粒子逃逸率反映布朗粒子运动规律，噪声引起布朗粒子跃迁Kramers逃逸率[20]为

(5)

(6)

布朗粒子在某一势阱中的平均驻留时间TK=1/rk与势函数变化时间相等才使得随机共振发生

TK=T/2

(7)

(8)

由于噪声太大信号被淹没,从波形图和频谱图中无法识别，如图4(a)和4(b)所示。经过随机共振后输出波形接近方波如图4(c)所示，0.01 Hz信号频谱幅度为87，如图4(d)所示，相对于周围噪声幅度信号突出，说明指数型双稳随机共振使信号能量大幅度提高。

2.2 系统参数对双稳随机共振影响

传统Duffing双稳随机共振在噪声系数过大下，增大阻尼比k，达到调节粒子跃迁速度使之与势函数周期变化趋势同步的目的，产生共振现象。指数型Duffing方程随机共振条件为式(8)，在满足绝热近似理论条件下定义函数F表达式如下

F(a,u,V,D,k)=

(a)输入信号波形

(b)输入频谱

(c)输出信号波形

(d)输出频谱

(9)

固定一组参数采用数值仿真，信号采用幅度0.1，频率0.01 Hz正弦周期信号，在采样频率5 Hz下采样5 000点，a=1.74，V=1，R=0.8不同k值下仿真500次取平均SNR曲线图如图5(a)所示。由图5(a)可知随着D的增大输出信噪比出现先增大后减小趋势，不同阻尼系数k在同一噪声强度对应输出信噪比值一样大，D=0.31的为随机共振输出最高点，此时输出信噪比为-9.825 dB，相对于输入信号(信噪比-31.416 9 dB)信噪比提高21 dB。固定信号参数、采样频率、采样点、噪声强度D=0.31、V=1不变。R=0.8时参数a对输出信噪比的影响如图5(b)所示，随着a的增加，输出信噪比出现先增后减最后趋于平稳趋势。在较小a的情况下，双稳势函数势阱边缘平坦，粒子跃迁无法通过，即使存在噪声也无法驱动粒子完成跃迁。随着a的增大布朗粒子运动在一个势阱内局部运动，此时系统没达到跃迁条件，部分信号能量转化为噪声能量导致输出信噪比降低。a继续增大，系统条件转化为对信号有利条件，噪声能量逐渐转化为信号能量，当a达到某一值时达到跃迁条件随机共振现象发生。固定a=1.75，图5(c)为参数R对输出信噪比影响曲线，随着R的逐渐增加，输出信噪比先减后增大到极大值后再减小，随着R的增大势垒高度逐渐降低，势垒高度达到布朗粒子运动跃迁高度，达到随机共振发生条件此时三者协同作用下噪声能量最大化转化为信号能量。由此可通过调节噪声强度以及系统参数实现随机共振。

2.3 随机共振参数优化

当信号噪声固定时，需要调节系统参数达到最佳共振效果，单一的调节会忽略参数之间影响。目前随机共振的衡量指标主要有功率谱放大系数、信噪比增益、相关系数等等。其中，信噪比增益更能直观地反应随机共振效果，本文采用噪比增益作为衡量指标，其公式为[21]

(10)

式中，SNRout为输出信噪比，SNRin为输入信噪比。SNRgain大于1时表明随机共振系统对输入待测弱信号有明显的改善和增强作用。本文采用粒子群算法以信噪比增益为适应度函数，对每一组实现数据进行参数选择。

(a)D

(b)a

(c)R

步骤1 种群初始化。固定信号幅度，噪声强度，及阻尼比固定为0.5，设置种群数量，参数a，V，R搜索范围以及最大进化代数Tmax，参数最大搜索速度取最大调整步长的10%～25%，这里的最大调整步长指的是所设定的粒子位置范围的上限值减去粒子位置范围的下限值所得的差值。随机初始化搜索点的位置，并计算出其相应的个体极值，记录整个粒子群中个体极值最大的粒子序号，设置Nbest为该最大粒子的当前位置。

步骤2 评价每个粒子。适应度函数采用SNRgain公式，计算粒子适应度值。与该粒子当前个体极值进行比较，若大于后者，则更新粒子个体极值，若在该粒子的邻域内所有粒子的个体极值中最大的大于当前的Nbest，则设置Nbest为该粒子的位置，记录该粒子的序号，并更新Nbest的函数值。

步骤3 粒子的更新,采用粒子群算法迭代公式[22]更新粒子速度及位置。

步骤4 检验是否符合结束条件。判断当前的迭代次数是否达到最大进化代数Tmax，满足条件则停止迭代，并输出最优参数，否则转至步骤2。

步骤5 检测结果。根据对a，V，R优化输出的最优解，计算适应度值，输出最优a，V，R。

3 经验模态分解

经验模态分解(Empirical Mode Decomposition，EMD)是用来实现降低噪声的方法之一，通过数据特征时间尺度来获得固有波动模式，从而自适应选择函数实现信号分解。经验模态分解后可以明显的得到频率大、小的信号，频率大的信号优先被分解得出，频率低的随后被分解，所有的IMF分量叠加得到的信号为最原始信号。分解过程中产生分解量为模态函数(IMF)。IMF信号分量中，信号极值点和过零点数目相等或相差1，为了得到满足条件模态函数，需要多次筛选，具体算法如下[23]

步骤1 计算信号局部极大值极小值，通过3次样条插值法拟合出上下包络线平均值m1(t)，并且认为h1(t)=x(t)-m1(t)为残余分量。

步骤2 理想情况下h1为第一个IMF分量，判断h1是否为满足IMF分量，不满足需要反复筛分，接下来将h1作为新信号，重复上述步骤，循环k次后，得到IMF条件h1k(t)。其中筛选次数约束满足柯西准则[24]

(11)

式中，T为信号时间长度，ε为门限值范围为0.2～0.3。

步骤3 得到第一阶IMF的c1(t)即为h1k(t)，r1(t)=x(t)-c1(t)将c1(t)作为原始信号，反复重复上面两个步骤，得到c2(t),c3(t),…,cn(t)，和剩余分量rn(t)，分解结束条件为rn(t)单调，由此可以将信号分解为n个模态分量。

经验模态分解对于小波分解优点在于不需要预先设定基函数，只需根据信号自身特性进行平滑处理，最后达到降低噪声，检测出信号的效果。

4 仿真信号分析

为了验证本文所提指数型双稳随机共振系统在微弱信号检测中的效果，分别构造冲击衰减信号和谐波振动信号进行模拟实际工程仿真，更加贴近实际情况。

4.1 冲击衰减信号

机械运动轴承信号产生都是具有冲击性，这种信号淹没在噪声中更加难以发现，为模拟冲击震荡信号构造冲击衰减信号

sn(t)=H(t)×exp[-100(t-t/T/T)]×sin[1 200π(t-t/T/T)]+n(t)

(12)

式中，H(t)为阶跃函数，T=0.128 s是冲击周期，n(t)为均值为0方差为0.6的高斯噪声，采样点为2 048，采样频率fs=2 000 Hz，特征信号频率为f=1/T=7.8 Hz，固定H(t)=0.2构造周期冲击衰减信号如图6所示。图中显示周期冲击衰减信号加入噪声后无法分辨出信号，将含噪信号送入指数型双稳随机共振系统中，特征信号频率为大参数需要预处理变为小参数，进行二次采样处理大频率信号，二次采样频率为fsr=5 Hz，采样压缩比为R=fs/fsr=400，相当于待测信号压缩为0.018 75 Hz满足小参数条件。

(a)波形图

(b)频谱图

经过参数为(a,b)=(1,1)经典双稳随机共振后特征信号7.8 Hz频谱幅度为27.96，周围噪声依然很大，信号不能识别如图7所示。由于阻尼系数对随机共振影响不大，固定k=0.5使用3.3节寻优参数找到(a,V,R)=(1.528 3,3.512 7,0.151 8)进行指数型随机共振，将输出信号进行大频恢复，特征信号7.8 Hz频谱幅度为44.9，如图8所示，相对于周围噪声信号成分突出明显，证明此系统在振动信号中检测有用性。

图7 经典双稳随机共振输出Fig.7 The output signal spectrum after CBSR

4.2 谐波振动信号分析

谐波振动信号是一种典型的多频振动信号，通过分析这种振动信号，可以使我们准确了解机器的运行状态。构造一种包含两种谐波振动成分的含噪信号

图8 指数型随机共振输出Fig.8 The output signal spectrum after exponential SR

sn(t)=A1sin(2π×ft)+A2sin(2π×2ft)+n(t)

(13)

图9 含噪信号Fig.9 Signal with noise

将二次采样后含噪信号送入指数型双稳系统，固定k=0.5使用3.3节寻优参数找到(a,V,R)=(1.687 4,0.470 4,0.050 8)进行指数型随机共振，随机共振后信号频谱如图11所示，高频噪声降低，低频成分升高，振动信号成分能量增强。图11可知振动成分f频谱幅度为203.1、2f频谱幅度为176.8。上述两实验分别验证两种信号下指数型随机共振在信号检测方面优于经典双稳系统。

含噪信号在指数型双稳随机共振后信号能量得到较大增加，但周围依然存在大量噪声。经验模态分解是很好的降噪技术手段，能将信号的成分从高频到低频逐层分解，为了有效的检测到待测的多频信号，将指数型双稳随机共振输出信号进行EMD分解得到IMF成分如图12所示，IMF1～IMF3含有幅度太小的噪声，IMF4～IMF5含有f信号成分，IMF6含有2f信号成分IMF7～IMF8为剩余分量。将IMF4～IMF6含有信号部分叠加得到降低噪声后的合成信号如图13所示，可以明显的观察到一次谐波与二次谐波信号，这相对于图11更加准确的检测到信号。

图10 经典双稳随机共振输出Fig.10 The output signal spectrum after CBSR

图11 指数型随机共振输出Fig.11 The output signal spectrum after exponential SR

图12 EMD分解后IMFs频谱分量Fig.12 The spectrum of IMFs decomposed by EMD

因此在多频信号大噪声环境下，可以先通过随机共振进行部分能量转移，再将输出信号进行降噪处理提取待测信号成分进行合成，能得到比随机共振后更加明显信号成分。

图13 合成信号频谱Fig.13 The synthetic signal spectrum

5 工程应用

滚动轴承故障信号具有非平稳性、调制性、微弱性，由于常常淹没在强大的背景噪声中难以发现和提取[25]。轴承故障出现对工程应用有不可估量的损失常常发生在内圈、外圈、滚动体等。轴承故障信号存在往往以倍频形式，故障越严重倍频成分越多。本文提出基于指数型双稳随机共振系统与EMD联合检测微弱信号方法并对轴承的内圈故障信号进行验证。首先，用二次采样对不满足绝热近似理论的故障信号进行预处理；其次，在将信号送入随机共振系统后，参数寻优后找到最佳参数对信号进行随机共振；再次，将随机共振后信号通过经验模态分解，观察故障频率、二倍频信号和三倍频信号；最后，含有故障信号及其倍频信号的模态分量进行合成输出，实现故障信号检测。本文故障信号采用Case Western Reserve University电气工程实验室[26]公开数据，选用型号为6205-2RS JEM SKF的深沟球轴承为实验对象，其主要参数如表1所示。轴承转速r=1 796 r/min(29.93 Hz)，采样频率fs=12 kHz，二次采样频率5 Hz，采样点数N=6 000，轴承的故障特征频率如表2所示。外圈的故障特征频率为162.2 Hz，其二倍频为324.4 Hz。三倍频率为486.6 Hz。

表1 滚动轴承主要计算参数Tab.1 The main computation parameters of the rolling bearing

表2 滚动轴承故障特征频率Tab.2 Rolling bearing fault feature frequency

故障信号时域图与频谱图如图14(a)所示，轴承转动由于故障之处产生周期性冲击信号，2 000～4 000频段存在干扰信号，故障信号存在低频段。将信号二次采样预处理后送入参数为(a,V,R,k)=(2.530 1,3.524 2,1.758 1,0.5)指数型双稳随机共振系统，随机共振后输出信号频谱如图14(b)所示，特征160 Hz显现为故障频率(误差在允许范围内)，但其二倍频、三倍频特征没有显现。将随机共振后输出信号进行EMD分解如图15所示，IMF4～IMF3成分含有特征信号，IMF3～IMF2含有二倍频信号，IMF2含有三倍频信号。图16为IMF2～IMF4分量进行叠后合成信号，由此可见大噪声存在时，随机共振都无法识别多频成分，结合经验模态分解后可使没被检测出的信号得到检测。