采用粒子群算法的自适应变步长随机共振研究

张仲海，王 多，王太勇，林锦州，蒋永翔

(1.天津大学 机构理论与装备设计教育部重点实验室，天津 300072;2.航天晨光股份有限公司研究院，南京 211100)

自从Benzi等［1］研究古气象冰川问题提出随机共振(Stochastic Resonance，SR)概念以来，SR现象受到了广泛的关注。随机共振现象是一种非线性现象，它在一定条件下，将部分高频噪声能量转移到低频信号上，在降低噪声的同时能够使淹没于噪声中的弱信号得到共振加强，极大地提高输出信噪比，从而实现从强烈噪声干扰中检测微弱信号的目的。传统的随机共振受到绝热近似理论［2］的限制，只适用于小参数信号(信号幅值、信号频率、噪声强度远小于1)，这极大地制约了随机共振在工程实际中的应用。

为了突破随机共振只能适应于小参数的限制，许多学者进行了有益探索，取得不少阶段性研究成果［3－6］，为随机共振应用于工程实际提供了可能。例如，文献［7］提出了二次采样随机共振(TSSR)方法，即按照一定的变换尺度将较高采样信号变换为较低的二次采样信号，以满足小参数条件，随机共振输出后，再按变换尺度恢复至较高采样信号。文献［8］提出一种移频变尺度随机共振(FRSR)方法，即通过移频和尺度变换等手段压缩分析频率，从而使之满足绝热近似小参数条件。文献［9］则通过双稳系统参数的调节，从而达到共振状态。以上方法为随机共振应用于工程实测信号提供了理论基础，但在实际应用中，如何实现参数的自适应选取是一个难题。针对自适应随机共振的问题，文献［10］提出了一种基于近似熵测度的自适应随机共振方法，在固定双稳系统结构参数b的条件下，对变步长(或二次采样)随机共振［8，10］系统参数a和计算步长h进行优化，当参数b选取不合适时，则无法达到真正意义上最优的共振状态。文献［11］提出了一种基于遗传算法的自适应随机共振算法，能对双稳系统的结构参数a、b进行优化，但在处理大参数信号时需要在优化前手动进行移频和尺度变换等操作，当调制频率fc和变尺度压缩率R选取不恰当时会影响随机共振的输出效果。对于某一确定的含噪信号，需要对随机共振系统进行多参数的同步调节，当且仅当双稳系统参数、信号频率与噪声相互之间达到最佳协同时，才能达到最佳的共振状态。

变步长随机共振能突破绝热近似理论对于小参数的限制，而应用于大参数条件下的工程实测信号，但如何对结构参数a，b和计算步长h进行自适应选取仍是一个难题。粒子群算法［12－13］作为一种多变量全局优化方法，没有遗传算法的交叉和变异，因此效率高、速度快，具有很强的工程应用价值。本文以变步长双稳随机共振系统为研究对象，利用粒子群算法的全局并行搜索优化能力，以双稳系统输出的信噪比作为粒子群算法的适应度函数，对变步长双稳随机共振系统的结构参数a，b和计算步长h进行同步优化，使a，b和h等三个参数达到最佳的协同，从而最优地实现对大参数条件下微弱信号的自适应检测。

1 变步长随机共振的基本原理

1.1 双稳随机共振概述

在随机共振的研究中，非线性双稳系统通常由势函数表示:

其中，a和b是双稳系统的结构参数。

在不考虑噪声的情况下得到一维非线性双稳系统的确定性动力学系统方程:

由定态方程dx/dt=0得到双稳系统的三个解为一个不稳定的定态解x1=0和两个稳定的定态解x2，3=从物理意义上理解，双稳系统的响应输出可以表现为一个粒子在双稳势阱内的运动，如图1所示。

由图1可以看出，系统参数a和b不仅调节双稳系统势垒ΔU=a2/4b，而且还改变着双稳系统两个势阱之间的距离

图1 非线性双稳系统势函数Fig.1 Non-linear bi-stable system potential function

当以单周期正弦信号Asin(2πf0t)和白噪声n(t)为输入信号时，对应的Langevin方程为:

这里，

其中D是噪声强度。

从粒子动力学角度来考察，式(3)描述了一个过阻尼的质点布朗运动。当A和n(t)均为0，即没有调制和噪声作用时，根据系统不同的初始状态，质点将处于双稳势阱中的某一个。当A＞0且足够小时，在无噪情况下，由于外加信号的驱动作用，整个系统不再处于平衡状态，势阱按信号频率发生周期性倾斜，但是粒子的运动将被限制在某一势阱内。然而，在引入噪声的情况下，即使A＜Ac，甚至A≪Ac时，质点将越过势垒而从当前势阱跃迁到另一势阱。如果外加噪声的强度合适，这种大幅度跃迁与周期驱动力达到很好的协同，那么原来两个势阱之间的随机跃迁运动就会变成与周期调制信号频率相一致的有序跃迁运动，随机共振现象就发生了。

1.2 变步长随机共振及其存在的问题

式(3)是一种非线性随机微分方程，可通过四阶Runge-Kutta法进行数值求解，具体算法如下:

其中，n=1，2，…，N。

式(4)中，Sn和xn分别是双稳系统输入S(t)=Asin(2πf0t)+n(t)和输出X(t)的第n个采样值，h=1/fs(fs)为数值计算步长。

当待测信号为小参数信号，满足绝热近似理论时，无需改变计算步长仅通过调节双稳系统结构参数就达到很好的共振状态。而当待测信号为大参数信号时，由于绝热近似理论的限制，仅通过调节系统参数无法产生共振。为了克服绝热近似理论对小参数的限制，冷永刚等［7］提出了二次采样随机共振思想，即对于大参数的待测信号，通过对其进行二次采样，将其转化为小参数，经过双稳系统处理得到共振输出后，按相应的变换尺度恢复还原到实际信号频率。设有一大参数的含噪信号，信号频率为f，采样频率为fs，二次采样频率为fsr，通过改变式(4)中的计算步长h，使其等于二次采样频率的倒数，即h=1/fsr，把原信号频率f变换成f0=f·fsr/fs，通过式(4)求解，就可得出频率f0下的共振输出。在实际数值求解过程中，可先调节计算步长h使双稳系统达到共振状态，得到变换后的信号频率f0，再按变换尺度R=hfs还原恢复实际信号频率f=Rf0。以上便是变步长随机共振(SCSR)思想。

大量分析表明，对大参数信号进行变步长随机共振处理时，通过适当的改变计算步长h，能拓宽随机共振系统输出的频带宽度。不过仅改变计算步长h不能获得很好的输出，还需要结合非线性系统结构参数a和b的调节，以改善非线性系统特性，进而得到待检测的微弱信号。但是，如何对结构参数a，b和计算步长h进行联合调节，进而使a，b，h三个参数达到最佳协同，获得双稳系统的最优输出，目前还没有一个定量的规律可循，只能凭经验探索性地选取，给随机共振在工程实际中的应用带来很大的不便。

2 基于粒子群算法的自适应变步长随机共振

如何对变步长随机共振系统的结构参数a，b和计算步长h进行自适应选取是一个难题，本文利用粒子群算法来实现结构参数 a，b和计算步长 h的同步优化。

2.1 粒子群优化算法

Kennedy等［12－13］研究鸟群觅食行为时，受到启发而提出了粒子群优化算法(Particle Swarm Optimization，PSO)。该算法实现方便，与遗传算法相比需要设置的参数少，是一种高效、实用的搜索优化算法。

图2 粒子群算法优化搜索示意图Fig.2 Searching schematic diagram of PSO

在粒子群优化算法中，每个粒子表征某一确定的待优化问题函数的一个可能解。首先，初始化一组粒子种群，粒子的初始速度和位置随机产生，之后，粒子群追随当前最优粒子，不断更新自己的速度和位置，在多维解空间中搜索，经过若干次迭代找出最优解。在每次迭代中，粒子通过跟踪其自身的当前最优解(即个体极值Pbest)和整个粒子群的当前最优解(即全局极值Nbest)来更新自己的速度和位置，直到达到设定的最大迭代次数或找到最优解。粒子的优化运动轨迹如图2所示。

PSO算法通常的数学描述为:设在某D维空间中，种群 X=(x1，…，xi，…，xm)由 m 个粒子组成。xi(t)=(xi1，xi2，…，xiD)T和 vi(t)=(vi1，vi2，…，viD)T分别表示t时刻种群中第i个粒子的位置和速度;Pbesti(t)=(Pbesti1，Pbesti2，…，PbestiD)T表示 t时刻第 i个粒子的个体极值，t时刻整个种群的全局极值则由Nbest(t)=(Nbest1，Nbest2，…，NbestD)T表示。那么，粒子 xi将按式(5)来更新其速度和位置。

式中，j=1，2，…，D;t为当前进化代数;r1，r2为均匀分布于［0，1］之间的随机数;c1，c2为学习因子，通常取c1=c2=2;w为惯性权重，即保持原来速度的系数。

惯性权重的计算公式如下:

式中;wmax为惯性权重上限;wmin为惯性权重下限;t为当前进化代数;Tmax为最大进化代数。

粒子速度和位置更新方式如图3所示。

图3 粒子速度和位置更新方式示意图Fig.3 Concept of modification of a searching point by PSO

图3中，vPbesti表示粒子在解空间中朝着个体极值进行搜索的速度，对应于c1r1(t)［Pbesti(t)－xi(t)］速度部分;vNbesti表示粒子在解空间中朝着全局极值进行搜索的速度，对应于c2r2(t)［Nbest(t)－xi(t)］速度部分。

2.2 目标优化函数

粒子群算法是根据各个粒子的适应度大小来调整进化搜索能力的。而适应度函数和目标优化函数是相关的。此处采用的适应度函数即为目标优化函数——随机共振输出信噪比。本文的适应度函数为:

式中:sr(a，b，h)为变步长随机共振的输出结果;SNRout(sr(a，b，h))表示随机共振输出的信噪比。

已知随机共振输出信噪比定义［14］如下:

式中:F0为信号频率;S(F0)为信号功率;P为系统总功率，包括信号功率和噪声功率;P－S(F0)即为噪声功率。

设输入信号为Asin(2πF0t)+n(t)，该含噪信号经采样频率为Fs的采样得到长度为L的离散序列Zl。Zl经过二次采样频率为Fsr的变步长随机共振，输出信号sr(a，b，h)中频率分量F'0=F0Fsr/Fs对应于输入频率F0。设F'0分量的单边谱幅值为 X(k0)，且有 k0=LF'0/Fsr=LF0/Fs，由式(8)可得:

2.3 基于粒子群算法的自适应变步长随机共振流程

粒子群优化算法容易实现且具有全局并行搜索优化能力，本文利用粒子群优化算法以2.2节所述的随机共振输出信噪比为适应度函数，对变步长随机共振系统的结构参数a、b和计算步长h等三个参数进行同步优化。算法的流程如图4所示。

下面给出其实现的具体步骤:

步骤1:种群初始化。设置种群数量、参数a、b和h的搜索范围以及最大进化代数Tmax，最大搜索速度取最大调整步长的10%～20%，这里的最大调整步长指的是所设定的粒子位置范围的上限值减去粒子位置范围的下限值所得的差值。随机初始化搜索点的位置xi(0)和速度vi(0)，设置各个粒子的Pbesti坐标为其当前位置xi(0)，并计算出其相应的个体极值，记录整个粒子群中个体极值最大的粒子序号，设置Nbest为该最大粒子的当前位置。

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

步骤3:粒子的更新。根据式(5)更新所有粒子的速度和位置。

图4 粒子群优化算法流程图Fig.4 Flowchart of PSO

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

步骤5:检测结果。根据对a、b和h优化输出的最优解，对原始信号进行变步长随机共振，尺度恢复后得到最终的微弱信号检测结果。

3 仿真数据分析

基于粒子群算法的自适应变步长随机共振通过自动联合调节结构参数a，b和步长h，得到了自适应条件下双稳系统的最优输出。下面利用仿真数据证明该自适应方法的有效性。

设输入信号为u(t)=A0sin(2πf0t)+n(t)，其中，A0=0.2，f0=20 Hz，添加均值为0、方差 D=3.1 的高斯白噪声 n(t)，采样频率 fs=2 048 Hz，采样点数 n=2 048。该输入信号的原始时域波形及其频谱分别如图5(a)和图5(b)所示，由于强噪声的加入，从图5(a)上很难发现周期成分，在图5(b)上也很难辨认出60 Hz的频率分量。

图5 大参数数据的仿真结果Fig.5 Simulation results of large parameters data

现利用本文提出的自适应变步长随机共振方法对该信号进行处理。初始化种群:设置种群数量为40，a、b、h的搜索范围分别为［0.1，10］、［0.1，1 000］和［0.02，0.2］，最大搜索速度取最大调整步长的20%，即a、b、h 的最大搜索速度分别为 1.98、199.98 和 0.036，最大进化代数为150，最小错误标准为1×10－4。从图5c的收敛曲线可以看出，经过62次迭代，算法收敛，输出的最优参数分别为 a=6.33、b=217.39、h=0.041 3。将最优参数代入变步长随机共振系统，对原始信号进行随机共振处理后分别得到图5(d)和图5(e)所示的时域波形和频谱。从图5d可以看出，噪声已经被极大地削弱了，在图5(e)上可以发现频率为0.236 7 Hz的分量非常突出。按变换尺度 R=hfs=0.041 3×2 048=84.582 4还原恢复后，可以得到f0=Rf=84.582 4×0.236 7=20 Hz，其频率正好对应于原始信号中的20 Hz频率成分。仿真结果表明，将粒子群优化算法用于变步长随机共振系统能实现对参数a、b和h的自适应最优选取，从而快速有效的检测出大参数条件下的微弱信号。

4 工程应用

实验选用型号为6205－2RS JEM SKF的深沟球轴承，该轴承的尺寸和故障频率如表1和表2如示。

表1 滚动轴承6205－2RS的尺寸参数(厘米)Tab.1 Size parameters of 6205 －2RS

表2 滚动轴承6205－2RS的故障频率(转频的倍数)Tab.2 Defect frequencies of 6205 －2RS

使用电火花加工技术在该轴承内圈上布置了单点故障，故障直径为0.017 78 cm，根据表2可以算出该滚动轴承内圈故障的特征频率为156.14 Hz。试验中，该轴承用于支承电机轴，电机转速为1 730 r/min，使用加速度传感器采集振动信号，采样频率为Fs=12 kHz，采样点数n=8 192。图6(a)和图6(b)所示的分别为原始采样信号的时域波形和幅值谱，在图6(b)的频谱图上根本无法辨认156.14 Hz的滚动轴承内圈故障频率。由于原始信号中冲击成分较为明显，对其进行包络解调处理，所得的时域波形和幅值谱分别如图6(c)和图6(d)所示，从图6(d)的包络谱中可以看到156 Hz的故障频率，但谱线很不明显。

现对该滚动轴承故障信号的包络进行基于粒子群算法的自适应变步长随机共振处理。首先初始化粒子种群:设置种群数量为80，a、b和h的搜索范围分别为［0.01，30］、［0.01，15 000］和［0.002，0.8］，最大搜索速度取最大调整步长的20%，即a、b、h的最大搜索速度分别为 5.998、2 999.998 和 0.159 6，最大进化代数为200，最小错误标准为1×10－4。

从图7的收敛曲线可以看出，经过103次迭代，算法收敛，输出的最优参数分别为 a=8.15、b=2 220.6、h=0.048 7。将最优参数代入变步长随机共振系统，对包络信号进行随机共振处理后分别得到图6(e)和图6(f)所示的时域波形和频谱。对比图6(c)和图6(e)后可以发现，随机共振处理后的时域波形中的噪声成分被极大的削弱了。在图6(f)中可以非常清楚的看到0.266 8 Hz的频率分量及其二倍频，按变换尺度R=hFs=0.048 5×12 000=584.4 还原恢复后，可以得到F0=RF=584.4 ×0.266 8≈155.9 Hz，考虑到存在计算舍入误差，即该频率正好是轴承内圈的故障特征频率，这与滚动轴承存在内圈故障的事实相吻合。

图6 滚动轴承内圈故障的诊断结果Fig.6 Diagnosis results of the fault of inner ring

图7 PSO算法的最优收敛曲线Fig.7 The optimal convergence curve of PSO algorithm

5 结论

本文提出了一种采用粒子群优化算法的自适应变步长随机共振方法，选用共振输出的信噪比作为粒子群算法的适应度函数，利用粒子群算法的全局并行搜索优化能力，对变步长随机共振系统的多个参数进行同步优化，实现了结构参数a、b和计算步长h等三个参数的自适应选取，从而最优检测出强噪声背景下的高频微弱信号。该方法克服了单参数优化及步长选取依赖经验的缺点，充分体现了联合调参的思想。将该方法用于仿真微弱信号的检测及滚动轴承的故障诊断，结果表明所提方法能快速有效地检测出大参数条件下的微弱信号。

［1］Benzi R，Parisi G，VulpianiA. Theoryofstochastic resonance in climatic chaner［J］.SIAM Journal on Applied Mathematics，1983，43(3):565 －578.

［2］Mcnamara B， Wiesenfeld K， RoyR. Observation of Stochastic Resonance in a Ring Laser［J］.Physical Review Letters.1988，60(25):2626 －2629.

［3］夏均忠，刘远宏，马宗坡，等.基于调制随机共振的微弱信号检测研究［J］.振动与冲击，2012，31(3):132 －135，140.XIA Jun-zhong，LIU Yuan-hong，MA Zong-po，et al.Weak signal detction based on the modulated stochastic resonance［J］.Journal of Vibration and Shock.2012，31(3):132 －135，140.

［4］李晓龙，冷永刚，范胜波，等.基于非均匀周期采样的随机共振研究［J］.振动与冲击，2011，30(12):78 －84.LI Xiao-long， LENG Yong-gong，FAN Sheng-bo，et al.Stochastic resonance based on periodic non-uniform sampling［J］.Journal of Vibration and Shock.2011，30(12):78 －84.

［5］郝 研，王太勇，万 剑，等.基于级联双稳随机共振和多重分形的机械故障诊断方法研究［J］.振动与冲击，2012，31(8):181－185.HAO Yan，WANG Tai-yong，WAN Jian，et al.Mechanical faultdiagnosis based on cascaded bistable stochastic resonance and multi-fractal［J］.Journal of Vibration and Shock，2012，31(8):181 －185.

［6］赵艳菊，王太勇，徐 跃，等.双稳随机共振降噪下的经验模式分解研究［J］.振动与冲击，2009，28(3):149－151.ZHAO Yan-ju，WANG Tai-yong，XU Yue，et al.Empirical mode decomposition based on bistable stochastic resonance denoising［J］.Journal of Vibration and Shock，2009 ，28(3):149－151.

［7］冷永刚，王太勇.二次采样用于随机共振从强噪声中提取弱信号的数值研究［J］.物理学报，2003，52(10):2432－2437.LENG Yong-gang，WANG Tai-yong.Numerical research of twice sampling stochastic resonance for the detection of a weak signal submerged in a heavy Noise［J］.Acta Physica Sinica，2003，52(10):2432 －2437.

［8］Tan J，Chen X，Wang J，et al.Study of frequency-shifted and re-scaling stochastic resonance and its application to fault diagnosis［J］.Mechanical Systems and Signal Processing，2009，23(3):811－822.

［9］陈 敏，胡茑庆，秦国军，等.参数调节随机共振在机械系统早期故障检测中的应用［J］.机械工程学报，2009，45(4):131－135.CHEN Min，HU Niao-qing，QIN Guo-jun，et al.Application of parameter-tuning stochastic resonance for detecting early mechanical faults［J］.Journal of Mechanical Engineering，2009，45(4):131－135.

［10］Qiang L，Taiyong W，Yonggang L，et al.Engineering signal processing based on adaptive step－changed stochastic resonance［J］.Mechanical Systems and Signal Processing，2007，21(5):2267－2279.

［11］Wang J，Zhang Q，Xu G H.Genetic Stochastic Resonance:A new fault diagnosis method to detect weak signals in mechanical systems［J］.Advanced Science Letters，2011，4(6－7):2508－2512.

［12］Kennedy J，Eberhart R.Particle swarm optimization［C］.Perth，Aust:IEEE，1995.

［13］Eberhart R，Kennedy J.A new optimizer using particle swarm theory［C］.New York，NY，USA:IEEE，1995.

［14］万 频，詹宜巨，李学聪，等.一种单稳随机共振系统信噪比增益的数值研究［J］.物理学报，2011，60(4):60－66.WAN Pin，ZHAN Yi-ju，LI Xue-cong，et al.Numerical research ofsignal-to-noise ratio gain on a monostable stochastic resonance［J］.Acta Physica Sinica，2011，60(4):60－66.