CN102735330A - 基于粒子群优化算法的自适应随机共振微弱信号检测方法 - Google Patents

Abstract

一种基于粒子群优化算法的自适应随机共振微弱信号检测方法：1）粒子种群的初始化；2）变步长随机共振；3）个体适应度评价；4）更新粒子的速度和位置；5）终止条件判定；6）输出检测结果。本发明简单易行，适用范围广，收敛速度快，能有效地检测出强噪声背景下的高频微弱信号，为随机共振参数的自适应选取和在工程中的实际应用提供了一种新的方法。

Description

基于粒子群优化算法的自适应随机共振微弱信号检测方法

技术领域

本发明涉及一种微弱信号检测方法。特别是涉及一种可以自适应地选取最优的变步长随机共振系统结构参数和计算步长，并能有效地检测出大参数条件下的微弱信号的基于粒子群优化算法的自适应随机共振微弱信号检测方法。

背景技术

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

目前，针对绝热近似理论大参数条件下的随机共振研究很多，取得了不少阶段性研究成果，大大拓宽了随机共振在工程实测信号处理中的应用，例如，具有代表性有二次采样随机共振(TSSR)方法、移频变尺度随机共振（FRSR）方法等。以上方法为随机共振应用于工程信号提供了理论基础，但在实际应用中，如何实现参数的自适应选取是一个难题。针对自适应随机共振的问题，现有的一些方法虽能在一定程度上实现参数的自适应选取，但这些方法要么是在固定某一结构参数的条件下实现对另一结构参数的优化，要么就是在事先设定好尺度变换等参数的前提下对结构参数进行优化，往往无法达到真正意义上的共振状态。而对于某一确定的含噪信号，需要对随机共振系统进行多参数同步调节，当且仅当双稳系统参数、信号频率与噪声相互之间达到最佳协同时，才能达到最佳的共振状态。

针对强噪声背景下、大参数条件下的微弱信号检测，变步长随机共振能突绝热近似理论对小参数的限制，从而应用于大参数条件下的工程实际信号，但如何对结构参数a、b和计算步长h进行自适应选取仍是一个难题，这限制了随机共振在工程实际中的进一步应用。

发明内容

本发明所要解决的技术问题是，提供一种将粒子群优化算法应用于变步长随机共振系统，为大参数条件特别是强噪声背景下高频微弱信号的检测提供了一种简单、快速、有效解决方案的基于粒子群优化算法的自适应随机共振微弱信号检测方法。

本发明所采用的技术方案是：一种基于粒子群优化算法的自适应随机共振微弱信号检测方法，包括如下阶段：

1）粒子种群的初始化，包括如下步骤：

（1）设置种群维数、种群数量、双稳系统结构参数a、b和计算步长h的搜索范围以及最大迭代次数T_max；

（2）设置最大搜索速度为最大调整步长的10%～20%；

（3）搜索点的初始位置x_i(0)和初始速度υ_i(0)在设定的范围内随机产生；

（4）将矩阵Pbest的各列依次设置为该列所对应粒子的当前位置x_i(0)，并计算出各粒子的个体极值，记录整个粒子群中个体极值最大的粒子序号，将向量Gbest设置为该最大粒子的当前位置。

2）变步长随机共振

将待测数据作为双稳系统输入S(t)，依次将当前各个粒子所对应的一组双稳系统结构参数a、b和计算步长h的具体值代入如下公式：

$\left\{\begin{array}{c}{x}_{n+1}=x_n+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)\\ k_1=h({\mathrm{ax}}_n-{\mathrm{bx}}_n^3+S_n)\\ k_2=h[a(x_n+\frac{k_1}{2})-b{(x_n+\frac{k_1}{2})}^3+S_n]\\ k_3=h[a(x_n+\frac{k_2}{2})-b{(x_n+\frac{k_2}{2})}^3+S_{n+1}]\\ k_4=h[a(x_n+k_3)-b{(x_n+k_3)}^3+S_{n+1}]\end{array}\right.$ 分别求解出各个粒子所对应的共振输出，

公式中：n=1,2，…，N；S_n和x_n分别是双稳系统输入S(t)=u(t)+n(t)和输出X(t)的第n个采样值；h=1/f_s为数值计算步长，其中f_s为采样频率；

3）个体适应度评价

根据适应度函数，计算每个粒子的适应度值，如果该适应度值大于该粒子当前的个体极值，则将该粒子所对应的Pbest列设置为该粒子的位置，并更新个体极值；若在该粒子的邻域内所有粒子的个体极值中最大的大于当前Gbest所对应的全局极值，则将Gbest设置为该粒子的位置，记录该粒子的序号，并更新Gbest的函数值。

4）更新粒子的速度和位置

根据下式对每一个粒子的速度和位置进行更新：

$\left\{\begin{array}{c}{v}_{\mathrm{ij}}(t+1)={\mathrm{wv}}_{\mathrm{ij}}\left(t\right)+c_1{r}_1\left(t\right)(p_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right))+c_2{r}_2\left(t\right)(p_{\mathrm{gj}}\left(t\right)-x_{\mathrm{ij}}\left(t\right))\\ x_{\mathrm{ij}}(t+1)=x_{\mathrm{ij}}\left(t\right)+v_{\mathrm{ij}}(t+1)\end{array}\right.$

式中，j=1,2，…，D；i=1,2，…，m；t为当前进化代数；r₁，r₂为均匀分布于[0,1]之间的随机数；c₁，c₂为学习因子，通常取c₁=c₂=2；υ_ij(t)为第i粒子的速度矢量中第j个元素的当前值；x_ij为第i粒子的位置矢量中第j个元素的当前值；υ_ij(t+1)为第i粒子的速度矢量中第j个元素更新后的值；x_ij(t+1)为第i粒子的位置矢量中第j个元素更新后的值；w为惯性权重，即保持原来速度的系数；

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

$w\left(t\right)=w_{\max }-\frac{(w_{\max }-w_{\min })\×t}{T_{\max }}$

式中，w_max为惯性权重上限；w_min为惯性权重下限；t为当前进化代数；T_max为最大进化代数；

5）终止条件判定

如果当前的迭代次数达到了预先设定的最大次数或达到最小错误要求，则停止迭代，输出最优解，否则转到步骤2。

6）输出检测结果

根据优化输出的一组最优双稳系统结构参数a、b值和计算步长h值，对待测数据进行变步长随机共振处理，尺度恢复后得到最终的微弱信号检测结果。

由于优化参数有分别是双稳系统结构参数a、b和计算步长h，故种群维数固定为3；所述的双稳系统结构参数a、b和计算步长h依次对应于粒子位置向量的三个元素。

第1）阶段中所述的设置种群数量和参数a、b、h的搜索范围时，需要根据优化要求以及待测数据的特点进行设置。

第1）阶段中所述的最大调整步长指的是设定的粒子位置范围的上限值减去粒子位置范围的下限值所得的差值。

第1）阶段中所述的矩阵Pbest为一D×m矩阵，Pbest的每一列依次存放第1个粒子至第m个粒子各自个体极值的位置，列向量Gbest为一D维列向量，Gbest存放整个种群全局极值的位置，其中，D是种群维数，m是种群数量。

第3）阶段中是以随机共振输出的信噪比为适应度函数，此处采用的适应度函数即为目标优化函数——随机共振输出信噪比，适应度函数具体为：

F(a，b，h)=SNR_out(sr(a，b，h))

式中：sr(a，b，h)为变步长随机共振的输出结果；SNR_out(sr(a，b，h))表示随机共振输出的信噪比；

已知随机共振输出信噪比定义如下：

$\mathrm{SNR}=10\log \frac{S\left(F_0\right)}{P-S\left(F_0\right)}$

式中：F₀为信号频率；S(F₀)为信号功率；P为系统总功率，包括信号功率和噪声功率；P-S(F₀)即为噪声功率；

设输入信号为Asin(2πF₀t)+n(t)，该含噪信号经采样频率为F_s的采样得到长度为L的离散序列Z_l。Z_l经过二次采样频率为F_sr的变步长随机共振，输出信号sr(a，b，h)中频率分量F′₀=F₀F_sr/F_s对应于输入频率F₀，设F′₀分量的单边谱幅值为X(k₀)，且有k₀=LF′₀/F_sr=LF₀/F_s，由公式 $\mathrm{SNR}=10\log \frac{S\left(F_0\right)}{P-S\left(F_0\right)}$ 得到：

${\mathrm{SNR}}_{\mathrm{out}}\left(\mathrm{sr}(a,b,h)\right)=10\log \frac{2{\left|X\left(k_0\right)\right|}^2}{{\mathrm{\Σ}}_{k=0}^{L-1}{\left|X\left(k\right)\right|}^2-2{\left|X\left(k_0\right)\right|}^2}.$

本发明的基于粒子群优化算法的自适应随机共振微弱信号检测方法，简单易行，适用范围广，收敛速度快，能有效地检测出强噪声背景下的高频微弱信号，为随机共振参数的自适应选取和在工程中的实际应用提供了一种新的方法。

附图说明

图1是粒子群优化算法流程图；

图2是粒子群优化算法优化搜索示意图；

图3是粒子速度和位置更新方式示意图；

图4是滚动轴承内圈故障的原始信号时域波形；

图5是滚动轴承内圈故障的原始信号频谱；

图6是针对滚动轴承内圈故障信号的粒子群优化算法的最优收敛曲线；

图7是滚动轴承内圈故障的随机共振输出时域波形；

图8是滚动轴承内圈故障的随机共振输出频谱。

具体实施方式

下面结合实施例和附图对本发明的基于粒子群优化算法的自适应随机共振微弱信号检测方法做出详细说明。

图1所示的是本发明的基本流程图，下面结合附图对本发明做进一步详细描述：本发明的基于粒子群优化算法的自适应随机共振微弱信号检测方法，包括如下阶段：

1）粒子种群的初始化

将粒子群优化算法应用于变步长随机共振，首先需要初始化一组粒子种群。

粒子群优化算法（Particle Swarm Optimization,PSO）是Kennedy和Eberhart于1995年研究鸟群觅食行为受到启发而提出的。该算法实现方便，与遗传算法相比需要设置的参数少，是一种高效、实用的搜索优化算法。

在粒子群优化算法中，每个粒子表征某一确定的待优化问题函数的一个潜在解。初始化一组粒子种群，粒子的初始速度和位置随机产生，之后，粒子群追随当前最优粒子在多维解空间中搜索，经过若干次迭代找出最优解。在每次迭代中，粒子通过跟踪两个“极值”来更新自己，第一个就是粒子本身所找到的最优解，这个解叫做个体极值，另一个极值是整个种群目前找到的最优解，这个极值是全局极值。粒子始终跟随这两个极值变更自己的位置和速度直到找到最优解。粒子的优化运动轨迹如图2所示。

粒子群优化算法数学描述为：

设在某D维空间中，种群X=(x₁，…，x_i，…x_m)由m个粒子所组成，其中第i个粒子位置为x_i=(x_i1，x_i2，…，x_iD)^T，其速度为v_i=(v_i1，v_i2…，v_id，…v_iD)^T，它的个体极值为p_i=(p_i1，p_i2，…p_iD)^T，种群的全局极值为p_g=(p_g1，p_g2，…p_gD)^T。Pbest为一D×m矩阵，Pbest的每一列依次存放第1个粒子至第m个粒子各自个体极值的位置；Gbest为一D维列向量，Gbest存放整个种群全局极值的位置。

粒子群的初始化操作包括如下步骤：

（1）设置种群维数、种群数量、双稳系统结构参数a、b和计算步长h的搜索范围以及最大迭代次数T_max；

（2）设置最大搜索速度为最大调整步长的10%～20%；

（3）搜索点的初始位置x_i(0)和初始速度υ_i(0)在设定的范围内随机产生；

（4）将矩阵Pbest的各列依次设置为该列所对应粒子的当前位置x_i(0)，并计算出各粒子的个体极值，记录整个粒子群中个体极值最大的粒子序号，将向量Gbest设置为Nbest_i该最大粒子的当前位置。

由于本发明涉及到的优化参数有3个，分别是双稳系统结构参数a、b和计算步长h，故种群维数固定为3；双稳系统结构参数a、b和计算步长h依次对应于粒子位置向量的三个元素。则对于第i个粒子，其位置为x_i=(x_i1，x_i2，x_i3)^T，其中，x_i1对应于a，x_i1对应于b，x_i1对应于h。

在设置种群数量和双稳系统结构参数a、b和计算步长h的搜索范围时，需要根据优化要求以及待测数据的特点综合考虑。当种群数量和双稳系统结构参数a、b和计算步长h的搜索范围设置的比较大时，有利于得到全局最优解，但收敛速度较慢；当种群数量和双稳系统结构参数a、b和计算步长h的搜索范围设置的比较小时，算法能很快收敛，但容易陷入局部最优。

最大调整步长指的是设定的粒子位置范围的上限值减去粒子位置范围的下限值所得的差值。

2）变步长随机共振

对当前粒子群中的各个粒子进行变步长随机共振处理，得到各个粒子的共振输出。

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

$\frac{\mathrm{dx}}{\mathrm{dt}}=\mathrm{ax}-{\mathrm{bx}}^3+u\left(t\right)+n\left(t\right)---\left(1\right)$

式中，a和b是双稳系统的结构参数；u(t)为微弱周期信号；n(t)是一定噪声强度下均值为0的高斯分布白噪声。

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

$\left\{\begin{array}{c}{x}_{n+1}=x_n+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)\\ k_1=h({\mathrm{ax}}_n-{\mathrm{bx}}_n^3+S_n)\\ k_2=h[a(x_n+\frac{k_1}{2})-b{(x_n+\frac{k_1}{2})}^3+S_n]\\ k_3=h[a(x_n+\frac{k_2}{2})-b{(x_n+\frac{k_2}{2})}^3+S_{n+1}]\\ k_4=h[a(x_n+k_3)-b{(x_n+k_3)}^3+S_{n+1}]\end{array}\right.---\left(2\right)$

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

式(2)中，S_n和x_n分别是双稳系统输入S(t)=u(t)+n(t)和输出X(t)的第n个采样值，h=1/f_s（f_s为采样频率）为数值计算步长。

当待测信号满足绝热近似理论小参数条件时，计算步长直接取采样频率的倒数，通过调节双稳系统结构参数就能达到很好的共振状态。而当待测信号为大参数信号时，仅通过调节系统参数无法产生共振。变步长随机共振能突破绝热近似理论仅适用于小参数的限制，而应用于工程实际中的大参数信号。

变步长随机共振（SCSR）思想为：设有一大参数的含噪信号，信号频率为f，采样频率为f_s，二次采样频率为f_sr，通过改变(4)式中的计算步长h，使其等于二次采样频率的倒数，即h=1/f_sr，把原信号频率f变换成f₀=f·f_sr/f_s=f/hf_s，通过（2）式求解，就可得出频率f₀下的共振输出。在实际数值求解过程中，可先调节计算步长h使双稳系统达到共振状态，得到变换后的信号频率f₀，再按变换尺度R=hf_s还原恢复实际信号频率f=Rf₀。

每个粒子对应一组双稳系统结构参数a、b和计算步长h，对各个粒子进行变步长随机共振处理的具体是：将待测数据作为双稳系统的输入S(t)，依次将各个粒子所对应的一组a、b、h值代入式(2)，分别求解出各个粒子所对应的共振输出。

3）个体适应度评价

根据适应度函数，计算每个粒子的适应度值，如果该适应度值大于该粒子当前的个体极值，则将该粒子所对应的Pbest列设置为该粒子的位置，并更新个体极值；若在该粒子的邻域内所有粒子的个体极值中最大的大于当前Gbest所对应的全局极值，则将Gbest设置为该粒子的位置，记录该粒子的序号，并更新Gbest的函数值。

本发明以随机共振输出的信噪比为适应度函数，适应度函数具体表述如下：

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

F(a，b，h)=SNR_out(sr(a，b，h))                                           （3）

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

已知随机共振输出信噪比定义如下：

$\mathrm{SNR}=10\log \frac{S\left(F_0\right)}{P-S\left(F_0\right)}---\left(4\right)$

式中：F₀为信号频率；S(F₀)为信号功率；P为系统总功率，包括信号功率和噪声功率；P-S(F₀)即为噪声功率。

设输入信号为Asin(2πF₀t)+n(t)，该含噪信号经采样频率为F_s的采样得到长度为L的离散序列Z_l。Z_l经过二次采样频率为F_sr的变步长随机共振，输出信号sr(a，b，h)中频率分量F′₀=F₀F_sr/F_s对应于输入频率F₀。设F′₀分量的单边谱幅值为X(k₀)，且有k₀=LF′₀/F_sr=LF₀/F_s，由(4)式可得：

${\mathrm{SNR}}_{\mathrm{out}}\left(\mathrm{sr}(a,b,h)\right)=10\log \frac{2{\left|X\left(k_0\right)\right|}^2}{{\mathrm{\Σ}}_{k=0}^{L-1}{\left|X\left(k\right)\right|}^2-2{\left|X\left(k_0\right)\right|}^2}---\left(5\right)$

4）更新粒子的速度和位置

按(6)式更新所有粒子的速度和位置，更新方式示意图如3所示。

$\left\{\begin{array}{c}{v}_{\mathrm{ij}}(t+1)={\mathrm{wv}}_{\mathrm{ij}}\left(t\right)+c_1{r}_1\left(t\right)(p_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right))+c_2{r}_2\left(t\right)(p_{\mathrm{gj}}\left(t\right)-x_{\mathrm{ij}}\left(t\right))\\ x_{\mathrm{ij}}(t+1)=x_{\mathrm{ij}}\left(t\right)+v_{\mathrm{ij}}(t+1)\end{array}\right.---\left(6\right)$

式中，j=1,2，…，D；i=1,2，…，m；t为当前进化代数；r₁，r₂为均匀分布于[0,1]之间的随机数；c₁，c₂为学习因子，通常取c₁=c₂=2；υ_ij(t)为第i粒子的速度矢量中第j个元素的当前值；x_ij为第i粒子的位置矢量中第j个元素的当前值；υ_ij(t+1)为第i粒子的速度矢量中第j个元素更新后的值；x_ij(t+1)为第i粒子的位置矢量中第j个元素更新后的值；w为惯性权重，即保持原来速度的系数。

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

$w\left(t\right)=w_{\max }-\frac{(w_{\max }-w_{\min })\×t}{T_{\max }}---\left(7\right)$

式中，w_max为惯性权重上限；w_min为惯性权重下限；t为当前进化代数；T_max为最大进化代数；

5）终止条件判定

如果当前的迭代次数达到了预先设定的最大次数或达到最小错误要求，则停止迭代，输出最优解，否则转到步骤2。

6）输出检测结果

根据优化输出的一组最优双稳系统结构参数a、b和计算步长h的具体数值，对待测数据进行变步长随机共振处理，尺度恢复后得到最终的微弱信号检测结果。

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

选用美国Case Western Reserve大学轴承数据中心的型号为6205-2RS JEM SKF的深沟球轴承，该轴承的尺寸和故障频率如表1和表2如示。使用电火花加工技术在该轴承内圈上布置了单点故障，故障直径为0.007英寸（1英寸=2.54厘米），根据表2可以算出该滚动轴承内圈故障的特征频率为156.14Hz。试验中，该轴承用于支承电机轴，电机转速为1730RPM，使用加速度传感器采集振动信号，采样频率为F_s=12kHz，采样点数n=2048。图4和图5所示的分别为原始采样信号的时域波形和频谱，在图5的频谱图上根本无法辨认156.14Hz的滚动轴承内圈故障频率。

表1滚动轴承6205-2RS的尺寸参数（英寸）

对该滚动轴承的故障数据进行基于粒子群优化算法的自适应变步长随机共振处理。初始化种群数量为40，a、b和h的搜索范围分别为[0.01,30]、[0.01,15000]和[0.002,0.8]，最大搜索速度为最大调整步长的20%，最大进行代数为300。从图6的收敛曲线可以看出，经过76次迭代，算法收敛，输出的最优参数分别为a=20、b=699.36、h=0.064。将最优参数代入变步长随机共振系统，对原始信号处理后分别得到图7和图8所示的输出时域波形和频谱。对比图4和图7后可以发现，处理后的时域波形中的冲击成分已相当明显，且噪声被极大的削弱了。在图8中可以非常清楚的看到0.204Hz的频率分量及其二倍频，按变换尺度R=hF_s=0.064×12000=768还原恢复后，可以得到F₀=RF=768×0.204=156.67Hz，即轴承内圈故障特征频率，这与滚动轴承存在内圈故障的事实相吻合。

Claims (6)

1.一种基于粒子群优化算法的自适应随机共振微弱信号检测方法，其特征在于，包括如下阶段：

1）粒子种群的初始化，包括如下步骤：

（1）设置种群维数、种群数量、双稳系统结构参数a、b和计算步长h的搜索范围以及最大迭代次数T_max；

（2）设置最大搜索速度为最大调整步长的10%～20%；

（3）搜索点的初始位置x_i(0)和初始速度υ_i(0)在设定的范围内随机产生；

（4）将矩阵Pbest的各列依次设置为该列所对应粒子的当前位置x_i(0)，并计算出各粒子的个体极值，记录整个粒子群中个体极值最大的粒子序号，将向量Gbest设置为该最大粒子的当前位置；

2）变步长随机共振

将待测数据作为双稳系统输入S(t)，依次将当前各个粒子所对应的一组双稳系统结构参数a、b和计算步长h的具体值代入如下公式：

$\left\{\begin{array}{c}{x}_{n+1}=x_n+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)\\ k_1=h({\mathrm{ax}}_n-{\mathrm{bx}}_n^3+S_n)\\ k_2=h[a(x_n+\frac{k_1}{2})-b{(x_n+\frac{k_1}{2})}^3+S_n]\\ k_3=h[a(x_n+\frac{k_2}{2})-b{(x_n+\frac{k_2}{2})}^3+S_{n+1}]\\ k_4=h[a(x_n+k_3)-b{(x_n+k_3)}^3+S_{n+1}]\end{array}\right.$ 分别求解出各个粒子所对应的共振输出，

公式中：n=1,2，…，N；S_n和x_n分别是双稳系统输入S(t)=u(t)+n(t)和输出X(t)的第n个采样值；h=1/f_s为数值计算步长，其中f_s为采样频率；

3）个体适应度评价

根据适应度函数，计算每个粒子的适应度值，如果该适应度值大于该粒子当前的个体极值，则将该粒子所对应的Pbest列设置为该粒子的位置，并更新个体极值；若在该粒子的邻域内所有粒子的个体极值中最大的大于当前Gbest所对应的全局极值，则将Gbest设置为该粒子的位置，记录该粒子的序号，并更新Gbest的函数值；

4）更新粒子的速度和位置

根据下式对每一个粒子的速度和位置进行更新：

$\left\{\begin{array}{c}{v}_{\mathrm{ij}}(t+1)={\mathrm{wv}}_{\mathrm{ij}}\left(t\right)+c_1{r}_1\left(t\right)(p_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right))+c_2{r}_2\left(t\right)(p_{\mathrm{gj}}\left(t\right)-x_{\mathrm{ij}}\left(t\right))\\ x_{\mathrm{ij}}(t+1)=x_{\mathrm{ij}}\left(t\right)+v_{\mathrm{ij}}(t+1)\end{array}\right.$

式中，j=1,2，…，D；i=1,2，…，m；t为当前进化代数；r₁，r₂为均匀分布于[0,1]之间的随机数；c₁，c₂为学习因子，通常取c₁=c₂=2；υ_ij(t)为第i粒子的速度矢量中第j个元素的当前值；x_ij为第i粒子的位置矢量中第j个元素的当前值；υ_ij(t+1)为第i粒子的速度矢量中第j个元素更新后的值；x_ij(t+1)为第i粒子的位置矢量中第j个元素更新后的值；w为惯性权重，即保持原来速度的系数；

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

$w\left(t\right)=w_{\max }-\frac{(w_{\max }-w_{\min })\×t}{T_{\max }}$

式中，w_max为惯性权重上限；w_min为惯性权重下限；t为当前进化代数；T_max为最大进化代数；

5）终止条件判定

如果当前的迭代次数达到了预先设定的最大次数或达到最小错误要求，则停止迭代，输出最优解，否则转到步骤2；

6）输出检测结果

根据优化输出的一组最优双稳系统结构参数a、b值和计算步长h值，对待测数据进行变步长随机共振处理，尺度恢复后得到最终的微弱信号检测结果。

2.根据权利要求1所述的基于粒子群优化算法的自适应随机共振微弱信号检测方法，其特征在于，由于优化参数有分别是双稳系统结构参数a、b和计算步长h，故种群维数固定为3；所述的双稳系统结构参数a、b和计算步长h依次对应于粒子位置向量的三个元素。

3.根据权利要求1所述的基于粒子群优化算法的自适应随机共振微弱信号检测方法，其特征在于，第1）阶段中所述的设置种群数量和参数a、b、h的搜索范围时，需要根据优化要求以及待测数据的特点进行设置。

4.根据权利要求1所述的基于粒子群优化算法的自适应随机共振微弱信号检测方法，其特征在于，第1）阶段中所述的最大调整步长指的是设定的粒子位置范围的上限值减去粒子位置范围的下限值所得的差值。

5.根据权利要求1所述的基于粒子群优化算法的自适应随机共振微弱信号检测方法，其特征在于，第1）阶段中所述的矩阵Pbest为一D×m矩阵，Pbest的每一列依次存放第1个粒子至第m个粒子各自个体极值的位置，列向量Gbest为一D维列向量，Gbest存放整个种群全局极值的位置，其中，D是种群维数，m是种群数量。

6.根据权利要求1所述的基于粒子群优化算法的自适应随机共振微弱信号检测方法，其特征在于，第3）阶段中是以随机共振输出的信噪比为适应度函数，此处采用的适应度函数即为目标优化函数——随机共振输出信噪比，适应度函数具体为：

F(a，b，h)=SNR_out(sr(a，b，h))

式中：sr(a，b，h)为变步长随机共振的输出结果；SNR_out(sr(a，b，h))表示随机共振输出的信噪比；

已知随机共振输出信噪比定义如下：

$\mathrm{SNR}=10\log \frac{S\left(F_0\right)}{P-S\left(F_0\right)}$

式中：F₀为信号频率；S(F₀)为信号功率；P为系统总功率，包括信号功率和噪声功率；P-S(F₀)即为噪声功率；

设输入信号为Asin(2πF₀t)+n(t)，该含噪信号经采样频率为F_s的采样得到长度为L的离散序列Z_l。Z_l经过二次采样频率为F_sr的变步长随机共振，输出信号sr(a，b，h)中频率分量F′₀=F₀F_sr/F_s对应于输入频率F₀，设F′₀分量的单边谱幅值为X(k₀)，且有k₀=LF′₀/F_sr=LF₀/F_s，由公式 $\mathrm{SNR}=10\log \frac{S\left(F_0\right)}{P-S\left(F_0\right)}$ 得到：

${\mathrm{SNR}}_{\mathrm{out}}\left(\mathrm{sr}(a,b,h)\right)=10\log \frac{2{\left|X\left(k_0\right)\right|}^2}{{\mathrm{\Σ}}_{k=0}^{L-1}{\left|X\left(k\right)\right|}^2-2{\left|X\left(k_0\right)\right|}^2}.$

Images