CN108645505A - 一种随机共振微弱信号检测方法 - Google Patents

Abstract

本发明公开了一种随机共振微弱信号检测方法，该随机共振微弱信号检测方法通过追尾行为的混沌变步长萤火虫优化算法寻找随机共振系统的结构最优参数，代入二维Duffing振子系统中，实现随机共振，检测出待检测信号中的微弱周期信号，输出最大信噪比，检测出非高斯噪声背景下的微弱信号。本发明将鱼群追尾行为、混沌和变步长操作引入到萤火虫算法中，解决萤火虫算法在寻优时出现拥挤的问题，在较小的规模内快速收敛到全局最优解，具有寻优能力强、收敛速度快、收敛精度高等优点。

Description

一种随机共振微弱信号检测方法

技术领域

本发明属于微弱信号检测领域，具体涉及一种基于追尾行为的混沌变步长萤火虫优化算法的随机共振微弱信号检测方法。

背景技术

实际工程应用领域中，强噪声背景下的微弱信号检测向来是国内外研究者的研究热点，所谓微弱信号检测，就是通过一系列信号处理方法有效滤除或抑制噪声，其目的是提高输出信号信噪比，检测出强噪声背景下的弱特征信号的技术，在通信、物理、化学、生物医学、遥感和材料学等领域，都存在处理微弱信号和噪声的问题。因此，有效提高输出信噪比是改善微弱信号检测侧的首要任务。

以往的随机共振微弱信号检测研究大多都建立在高斯白噪声模型背景下，但在研究高斯白噪声背景下的随机共振时，通常把高斯白噪声理想化，然而在实际应用中，经常会碰到具有拖尾特性和显著尖峰脉冲特性的非高斯噪声(如：海杂波)。因此，非高斯背景下的随机共振研究显得尤为重要，实现非高斯噪声背景下的微弱信号检测是急需解决的问题。

混沌变步长萤火虫优化算法是一种改进的萤火虫算法，但其在迭代后期会出现寻优精度低、易陷入局部最优、早熟收敛以及收敛速度不够快等缺点。

发明内容

本发明的目的是为了解决现有混沌变步长萤火虫优化算法中存在的缺陷，提供一种能够解决萤火虫算法在寻优时出现拥挤问题的改进算法。

为了达到上述目的，本发明提供了一种随机共振微弱信号检测方法，该随机共振微弱信号检测方法通过追尾行为的混沌变步长萤火虫优化算法寻找随机共振系统的结构最优参数，代入二维Duffing振子系统中，实现随机共振，检测出待检测信号中的微弱周期信号，输出最大信噪比，检测出非高斯噪声背景下的微弱信号。

具体检测方法包括以下步骤：

(1)PCVSGSO参数初始化；

(2)混沌初始化及最优适应度值初始化：按照二维混沌公式对各只萤火虫完成混沌初始化，得到初始质量较优和分布均匀的萤火虫，然后获得混沌优化后的萤火虫适应度值，将混沌优化后的萤火虫适应度值赋值为每只萤火虫的局部最佳适应度值，把所有局部最佳适应度值中的最大值赋值为萤火虫的全局最佳适应度值；

(3)荧光素更新；

(4)计算拥挤度：计算要移向目标萤火虫周围的拥挤度；

(5)更新位置和决策域半径：得到下一次萤火虫迭代的位置；更新萤火虫位置后，完成更新位置后萤火虫的动态决策域更新；

(6)更新最优适应度值：根据步骤(5)更新的萤火虫位置后，得到更新后各只萤火虫的适应度值；如果获得的单只萤火虫局部最佳适应度值或者全局最佳适应度值优于步骤(2)中得到的上一代萤火虫的单只局部最佳适应度值或者全局最佳适应度值，则进行更新；

(7)输出最优随机共振：当迭代次数达到最大迭代次数时，依据最终输出的全局最佳适应度值所对应的萤火虫位置获得的最终二维Duffing振子随机共振系统参数的寻优结果，对原始信号完成随机共振输出。

进一步的，本发明基于PCVSGSO的随机共振微弱信号检测方法具体步骤如下：

第一步:PCVSGSO参数初始化。给各只萤火虫分配同样大小的荧光素l₀与感应范围γ₀，形成萤火虫初始的种群。并且将萤火虫的初始移动步长设为s(0)、最大迭代次数T_max、维度D与萤火虫种群规模m。在搜索空间内随机产生初始的萤火虫位置向量：

x_i＝(x_i1,x_i2,...,x_iD)(i＝1,2,...,m) (1)

式(1)中：x_i是萤火虫种群内第i只萤火虫所在的位置

第二步：混沌初始化及最优适应度值初始化。选用一种二维混沌式(2)作为本文的混沌系统。对式(1)中产生的各只萤火虫位置按照二维混沌公式完成混沌初始化，可得到初始质量较优和分布均匀的萤火虫，此时所处位置为x_i+1。然后使用式(3)可获得混沌优化后的萤火虫适应度值，将混沌优化后的萤火虫适应度值赋值为每只萤火虫的局部最佳适应度值Pbest_i(0)(i＝1,2,...,m)，把Pbest_i(0)中的最大值赋值为萤火虫的全局最佳适应度值gbest_i(0)。其中，Pbest_i(t)＝(Pbest_i1,Pbest_i2,...,Pbest_iD)^T是第i只萤火虫在t时刻的个体极值，gbest_i(t)＝(gbest_i1,gbest_i2,...,gbest_iD)^T是所有种群在t时刻的全局极值。

x_i+1＝tan(αx_i)exp(-βx_i) (2)

式(2)中：α,β是控制参量，取α＝5，β＝3,假定0＜x₀＜1，则上面公式达到了混沌状态。

F(a,b,k)＝SNR_out(sr(a,b,k)) (3)

式(3)中：sr(a,b,k)是PCVSGSO的随机共振输出结果，SNR_out(sr(a,b,k))是SR的输出信噪比，也为适应度函数。

式(5)中:X(f₀)是输出信号的功率谱；N是背景的噪声谱，即ω₀频率周围段内的平均值。

第三步：荧光素更新。萤火虫种群中的每只萤火虫根据式(6)来更新荧光素。

l_i(t)＝(1-p)l_i(t-1)+γf(x_i(t)) (6)

式(6)中：l_i(t)表示在第t代第i只萤火虫的荧光素含量；p表示为荧光素挥发因子，取p＝0.4；γ表示为荧光素更新率，取γ＝0.6；f(x_i(t))表示第i只萤火虫在第t代时，处于x_i位置的目标函数，也就是适应度值函数。

第四步：计算拥挤度。通过公式(7)计算出要移向目标萤火虫周围的拥挤度，当q_i(t)＜δ(t)时，说明要移向的目标萤火虫周围不拥挤，则萤火虫能够直接飞向要移向的目标萤火虫，当q_i(t)＞δ(t)时，说明要移向的目标萤火虫周围拥挤，则该萤火虫在自己的感知半径内以概率P_ij(t)随机挑选一只萤火虫飞行。

式(7)中：表示在t时刻第i只萤火虫的动态决策域。

δ(t)＝l-e^-ct (8)

式(8)中：c为阈值变化系数。

式(9)中：P_ij(t)为第i只萤火虫向邻居集内第j只萤火虫转移概率，l_i(t)表示在t时刻第i只萤火虫的荧光素含量，N_i(t)表示在t时刻第i只萤火虫的邻居。

第五步：更新位置和决策域半径。根据式(10)得到下一次萤火虫迭代的位置，即目标萤火虫位置，式(10)中的s(t)是变化的步长，采用式(11)的关系变化，随着迭代次数的增加呈曲线下降，在迭代初期时，由于步长较大，萤火虫群体拥有较高的全局寻优能力，在迭代后期，步长减小，可增加群体的局部寻优能力。在更新萤火虫位置后，根据式(12)完成更新位置后萤火虫的动态决策域更新。

式(10)中：x_i(t+1)表示在t+1时刻第i只萤火虫的位置，s(t)表示变化步长。

式(11)中：s_min是s的最小值，取s_min＝10^-4，s_max为s的最大值，取s_max＝1，初始步长s(0)＝0.03，t_max为迭代次数的最大值，取t_max＝100，g(t)为当前迭代次数。

式(12)中：为迭代到t+1代时，第i只萤火虫的动态决策域半径；r_s表示单只萤火虫的感知半径；β表示动态决策域的更新率，取值为0.08；N_i(t)表示迭代到第t代时，在其动态决策域内比第i只萤火虫亮的萤火虫数目。

第六步：更新最优适应度值。按照式(10)迭代更新萤火虫的位置，利用式(4)得出更新后各只萤火虫的适应度值，如果获得的单只萤火虫局部最佳适应度值Pbest_i(t+1)或全局最佳适应度值gbest_i(t+1)好于上一代萤火虫的Pbest_i(t)和gbest_i(t)，那么更新Pbest_i(t)和gbest_i(t)。

第七步：输出最优随机共振。当迭代次数达到最大迭代次数T_max时，此时依据最终输出的全局最佳适应度值gbest_i(t)所对应的萤火虫位置获得最终的二阶Duffing振子随机共振系统参数k、a和b的寻优结果，且利用所获得的参数对原始信号完成随机共振输出。

其中，二维Duffing振子随机共振系统模型为：

对于一个受噪声信号n(t)与微弱周期信号u(t)作用的二阶Duffing振子随机共振方程：

式(13)中k是阻尼比；x是系统输出；V(x)是势函数；a和b是大于0的非线性系统结构参数,系统势函数V(x)在相面上有三个极值点，分别为稳定点 以及不稳定点(0,0)，势垒高度为ΔV＝a²/(4b)；u(t)是二阶Duffing振子SR的输入信号，A是输入信号的幅度；δ是单位脉冲信号；是输入信号的相位角；xDξ(t)为α稳定噪声，D代表α稳定噪声的强度放大系数。

其中，α稳定噪声数据通常由α稳定分布序列产生，α稳定分布通常用特征函数表达：

Φ(t)＝exp{jut-γ|t|^α[1+jBsign(t)ω(t,α)]} (14)

式(15)中，α∈(0,2]为特征指数，决定分布的脉冲特性和拖尾特性；γ为尺度参数；u为位置参数，表明分布的中心；B∈[-1,1]为对称参数，B＝0为完全对称，α稳定分布为特殊高斯模型，B≠0时为非对称，α稳定分布为非高斯模型。

本发明相比现有技术具有以下优点：

1、本发明将鱼群追尾行为、混沌和变步长操作引入到萤火虫算法中，解决萤火虫算法在寻优时出现拥挤的问题，在较小的规模内快速收敛到全局最优解，具有寻优能力强、收敛速度快、收敛精度高等优点。

2、本发明利用追尾行为的混沌变步长萤火虫优化算法寻找随机共振系统的结构最优参数速度快、准确度高，能达到系统参数高精度匹配的要求。将寻找到的最优参数代入二维Duffing振子随机共振系统中，实现随机共振，此时，输入微弱信号、非高斯白噪声、二阶Duffing非线性系统产生某种协同效应时，在低频处将噪声的部分能量转移到微弱周期信号中，检测出非高斯噪声背景下的微弱周期信号，此时提高了输出信噪比，增强微弱信号的检测精度。

附图说明

图1为本发明随机共振微弱信号检测方法的流程图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。

本发明基于追尾行为的混沌变步长萤火虫优化算法的随机共振微弱信号检测方法，是将经过追尾行为的混沌变步长萤火虫优化算法运算得到的优化参数代入到由α稳定噪声驱动的二维Duffing振子随机共振系统中，实现随机共振，进行微弱信号的检测。其中，PCVSGSO是一种利用混沌搜索改进的萤火虫优化算法，根据每一只萤火虫的适应度来更新荧光素，从而调节进化搜索能力，完成性能最优萤火虫的寻找。基于PCVSGSO的随机共振微弱信号检测方法流程图如图1所示，其具体步骤如下：

第一步:PCVSGSO参数初始化。给各只萤火虫分配同样大小的荧光素l₀与感应范围γ₀，形成萤火虫初始的种群。并且将萤火虫的初始移动步长设为s(0)、最大迭代次数T_max、维度D与萤火虫种群规模m。在搜索空间内随机产生初始的萤火虫位置向量：

x_i＝(x_i1,x_i2,...,x_iD)(i＝1,2,...,m) (1)

式(1)中：x_i是萤火虫种群内第i只萤火虫所在的位置

第二步：混沌初始化及最优适应度值初始化。选用一种二维混沌式(2)作为本文的混沌系统。对式(1)中产生的各只萤火虫位置按照二维混沌公式完成混沌初始化，可得到初始质量较优和分布均匀的萤火虫，此时所处位置为x_i+1。然后使用式(3)可获得混沌优化后的萤火虫适应度值，将混沌优化后的萤火虫适应度值赋值为每只萤火虫的局部最佳适应度值Pbest_i(0)(i＝1,2,...,m)，把Pbest_i(0)中的最大值赋值为萤火虫的全局最佳适应度值gbest_i(0)。其中，Pbest_i(t)＝(Pbest_i1,Pbest_i2,...,Pbest_iD)^T是第i只萤火虫在t时刻的个体极值，gbest_i(t)＝(gbest_i1,gbest_i2,...,gbest_iD)^T是所有种群在t时刻的全局极值。

x_i+1＝tan(αx_i)exp(-βx_i) (2)

式(2)中：α,β是控制参量，取α＝5，β＝3,假定0＜x₀＜1，则上面公式达到了混沌状态。

F(a,b,k)＝SNR_out(sr(a,b,k)) (3)

式(3)中：sr(a,b,k)是PCVSGSO的随机共振输出结果，SNR_out(sr(a,b,k))是SR的输出信噪比，也为适应度函数。

式(5)中:X(f₀)是输出信号的功率谱；N是背景的噪声谱，即ω₀频率周围段内的平均值。

第三步：荧光素更新。萤火虫种群中的每只萤火虫根据式(6)来更新荧光素。

l_i(t)＝(1-p)l_i(t-1)+γf(x_i(t)) (6)

式(6)中：l_i(t)表示在第t代第i只萤火虫的荧光素含量；p表示为荧光素挥发因子，取p＝0.4；γ表示为荧光素更新率，取γ＝0.6；f(x_i(t))表示第i只萤火虫在第t代时，处于x_i位置的目标函数，也就是适应度值函数。

第四步：计算拥挤度。通过公式(7)计算出要移向目标萤火虫周围的拥挤度，当q_i(t)＜δ(t)时，说明要移向的目标萤火虫周围不拥挤，则萤火虫能够直接飞向要移向的目标萤火虫，当q_i(t)＞δ(t)时，说明要移向的目标萤火虫周围拥挤，则该萤火虫在自己的感知半径内以概率P_ij(t)随机挑选一只萤火虫飞行。

式(7)中：表示在t时刻第i只萤火虫的动态决策域。

δ(t)＝1-e^-ct (8)

式(8)中：c为阈值变化系数。

式(9)中：P_ij(t)为第i只萤火虫向邻居集内第j只萤火虫转移概率，l_i(t)表示在t时刻第i只萤火虫的荧光素含量，N_i(t)表示在t时刻第i只萤火虫的邻居。

第五步：更新位置和决策域半径。根据式(10)得到下一次萤火虫迭代的位置，即目标萤火虫位置，式(10)中的s(t)是变化的步长，采用式(11)的关系变化，随着迭代次数的增加呈曲线下降，在迭代初期时，由于步长较大，萤火虫群体拥有较高的全局寻优能力，在迭代后期，步长减小，可增加群体的局部寻优能力。在更新萤火虫位置后，根据式(12)完成更新位置后萤火虫的动态决策域更新。

式(10)中：x_i(t+1)表示在t+1时刻第i只萤火虫的位置，s(t)表示变化步长。

式(11)中：s_min是s的最小值，取s_min＝10^-4，s_max为s的最大值，取s_max＝1，初始步长s(0)＝0.03，t_max为迭代次数的最大值，取t_max＝100，g(t)为当前迭代次数。

式(12)中：为迭代到t+1代时，第i只萤火虫的动态决策域半径；r_s表示单只萤火虫的感知半径；β表示动态决策域的更新率，取值为0.08；N_i(t)表示迭代到第t代时，在其动态决策域内比第i只萤火虫亮的萤火虫数目。

第六步：更新最优适应度值。按照式(10)迭代更新萤火虫的位置，利用式(4)得出更新后各只萤火虫的适应度值，如果获得的单只萤火虫局部最佳适应度值Pbest_i(t+1)或全局最佳适应度值gbest_i(t+1)好于上一代萤火虫的Pbest_i(t)和gbest_i(t)，那么更新Pbest_i(t)和gbest_i(t)。

第七步：输出最优随机共振。当迭代次数达到最大迭代次数T_max时，此时依据最终输出的全局最佳适应度值gbest_i(t)所对应的萤火虫位置获得最终的二阶Duffing振子随机共振系统参数k、a和b的寻优结果，且利用所获得的参数对原始信号完成随机共振输出。

其中，二维Duffing振子随机共振系统模型为：

对于一个受噪声信号n(t)与微弱周期信号u(t)作用的二阶Duffing振子随机共振方程：

式(13)中k是阻尼比；x是系统输出；V(x)是势函数；a和b是大于0的非线性系统结构参数,系统势函数V(x)在相面上有三个极值点，分别为稳定点 以及不稳定点(0,0)，势垒高度为ΔV＝a²/(4b)；u(t)是二阶Duffing振子SR的输入信号，A是输入信号的幅度；δ是单位脉冲信号；是输入信号的相位角；xDξ(t)为α稳定噪声，D代表α稳定噪声的强度放大系数。

其中，α稳定噪声数据通常由α稳定分布序列产生，α稳定分布通常用特征函数表达：

Φ(t)＝exp{jut-γ|t|^α[1+jBsign(t)ω(t,α)]} (14)

式(15)中，α∈(0,2]为特征指数，决定分布的脉冲特性和拖尾特性；γ为尺度参数；u为位置参数，表明分布的中心；B∈[-1,1]为对称参数，B＝0为完全对称，α稳定分布为特殊高斯模型，B≠0时为非对称，α稳定分布为非高斯模型。

Claims (9)

1.一种随机共振微弱信号检测方法，其特征在于：所述随机共振微弱信号检测方法通过追尾行为的混沌变步长萤火虫优化算法寻找随机共振系统的结构最优参数，代入二维Duffing振子系统中，实现随机共振，检测出待检测信号中的微弱周期信号，输出最大信噪比，检测出非高斯噪声背景下的微弱信号。

2.根据权利要求1所述的随机共振微弱信号检测方法，其特征在于：所述检测方法包括以下步骤：

(1)PCVSGSO参数初始化；

(2)混沌初始化及最优适应度值初始化：按照二维混沌公式对各只萤火虫完成混沌初始化，得到初始质量较优和分布均匀的萤火虫，然后获得混沌优化后的萤火虫适应度值，将混沌优化后的萤火虫适应度值赋值为每只萤火虫的局部最佳适应度值，把所有局部最佳适应度值中的最大值赋值为萤火虫的全局最佳适应度值；

(3)荧光素更新；

(4)计算拥挤度：计算要移向目标萤火虫周围的拥挤度；

(5)更新位置和决策域半径：得到下一次萤火虫迭代的位置；更新萤火虫位置后，完成更新位置后萤火虫的动态决策域更新；

(6)更新最优适应度值：根据步骤(5)更新的萤火虫位置后，得到更新后各只萤火虫的适应度值；如果获得的单只萤火虫局部最佳适应度值或者全局最佳适应度值优于步骤(2)中得到的上一代萤火虫的单只局部最佳适应度值或者全局最佳适应度值，则进行更新；

(7)输出最优随机共振：当迭代次数达到最大迭代次数时，依据最终输出的全局最佳适应度值所对应的萤火虫位置获得的最终二维Duffing振子随机共振系统参数的寻优结果，对原始信号完成随机共振输出。

3.根据权利要求2所述的随机共振微弱信号检测方法，其特征在于：所述二维Duffing振子随机共振系统模型为：

对于一个受噪声信号n(t)与微弱周期信号u(t)作用的二阶Duffing振子随机共振方程如下式(13)所示：

式(13)中k是阻尼比；x是系统输出；V(x)是势函数；a和b是大于0的非线性系统结构参数,系统势函数V(x)在相面上有三个极值点，分别为稳定点以及不稳定点(0,0)，势垒高度为ΔV＝a²/(4b)；u(t)是二阶Duffing振子SR的输入信号，A是输入信号的幅度；δ是单位脉冲信号；是输入信号的相位角；xDξ(t)为α稳定噪声，D代表α稳定噪声的强度放大系数。

4.根据权利要求3所述的随机共振微弱信号检测方法，其特征在于：所述α稳定噪声数据由α稳定分布序列产生，α稳定分布的特征函数表达如下：

Φ(t)＝exp{jut-γ|t|^α[1+jBsign(t)ω(t,α)]} (14)

式(15)中，α∈(0,2]为特征指数，决定分布的脉冲特性和拖尾特性；γ为尺度参数；u为位置参数，表明分布的中心；B∈[-1,1]为对称参数，B＝0时为完全对称，α稳定分布为特殊高斯模型，B≠0时为非对称，α稳定分布为非高斯模型。

5.根据权利要求2所述的随机共振微弱信号检测方法，其特征在于：所述步骤(2)中二维混沌公式如下式(2)所示：

x_i+1＝tan(αx_i)exp(-βx_i) (2)

式(2)中：α,β是控制参量，取α＝5，β＝3，0＜x₀＜1时式(2)达到了混沌状态；所述得到初始质量较优和分布均匀的萤火虫，此时所处位置为x_i+1；所述混沌优化后的萤火虫适应度值通过下式获得：

F(a,b,k)＝SNR_out(sr(a,b,k)) (3)

式(3)中：sr(a,b,k)是PCVSGSO的随机共振输出结果，SNR_out(sr(a,b,k))是SR的输出信噪比，也为适应度函数；式(5)中:X(f₀)是输出信号的功率谱，N是背景的噪声谱，即ω₀频率周围段内的平均值。

6.根据权利要求5所述的随机共振微弱信号检测方法，其特征在于：所述步骤(3)中荧光素根据下式(6)进行更新：

l_i(t)＝(1-p)l_i(t-1)+γf(x_i(t)) (6)

式(6)中：l_i(t)表示在第t代第i只萤火虫的荧光素含量；p表示为荧光素挥发因子，取p＝0.4；γ表示为荧光素更新率，取γ＝0.6；f(x_i(t))表示第i只萤火虫在第t代时，处于x_i位置的目标函数，也就是适应度值函数。

7.根据权利要求6所述的随机共振微弱信号检测方法，其特征在于：所述步骤(4)中要移向目标萤火虫周围的拥挤度通过下式(7)～(9)计算：

δ(t)＝1-e^-ct (8)

式(7)中：表示在t时刻第i只萤火虫的动态决策域，N_i(t)表示在t时刻第i只萤火虫的邻居；式(8)中：c为阈值变化系数；式(9)中：P_ij(t)为第i只萤火虫向邻居集内第j只萤火虫转移概率，l_i(t)表示在t时刻第i只萤火虫的荧光素含量；

当q_i(t)＜δ(t)时，说明要移向的目标萤火虫周围不拥挤，则萤火虫能够直接飞向要移向的目标萤火虫；当q_i(t)＞δ(t)时，说明要移向的目标萤火虫周围拥挤，则该萤火虫在自己的感知半径内以概率P_ij(t)随机挑选一只萤火虫飞行。

8.根据权利要求7所述的随机共振微弱信号检测方法，其特征在于：所述步骤(5)中下一次萤火虫迭代的位置由下式(10)计算得到：

式(10)中：x_i(t+1)表示在t+1时刻第i只萤火虫的位置，s(t)表示变化步长；该变化步长s(t)采用下式(11)的关系变化：

式(11)中：s_min是s的最小值，取s_min＝10^-4，s_max为s的最大值，取s_max＝1，初始步长s(0)＝0.03，t_max为迭代次数的最大值，取t_max＝100，g(t)为当前迭代次数；

所述步骤(5)中在更新萤火虫位置后，根据下式(12)完成更新位置后萤火虫的动态决策域更新：

式(12)中：为迭代到t+1代时，第i只萤火虫的动态决策域半径；r_s表示单只萤火虫的感知半径；β表示动态决策域的更新率，取值为0.08；N_i(t)表示迭代到第t代时，在其动态决策域内比第i只萤火虫亮的萤火虫数目。

9.根据权利要求1所述的随机共振微弱信号检测方法，其特征在于：所述PCVSGSO参数初始化过程如下：给各只萤火虫分配同样大小的荧光素l₀与感应范围γ₀，形成萤火虫初始的种群，并且将萤火虫的初始移动步长设为s(0)、最大迭代次数T_max、维度D与萤火虫种群规模m，在搜索空间内随机产生初始的萤火虫位置向量：

x_i＝(x_i1,x_i2,...,x_iD)(i＝1,2,...,m) (1)

式(1)中：x_i是萤火虫种群内第i只萤火虫所在的位置。