白晓波，邵景峰，和 征，田建刚

(1.西安工程大学 管理学院，陕西 西安 710048；2.陆军边海防学院 工程基础系，陕西 西安 710108)

0 引 言

很多学者对粒子滤波(particle filter，PF)[1,2]进行了深入研究。文献[3]利用萤火虫算法迭代寻优的思想优化粒子分布，但萤火虫算法本身存在着对于初始解分布的依赖性；文献[4]将粒子滤波算法用于嵌入式系统进行目标跟踪；K.Madhan Kumar等[5]提出混合多核偏最小二乘(PLS)粒子滤波方法，基于多核偏最小二乘法，结合粒子滤波算法提高遥感图像中道路提取的准确度。其它如文献[6-8]提出了相关的改进算法，都集中于粒子采样和粒子权值计算的改进；文献[9]主要解决粒子数的自适应选择问题，在保证适当粒子数的情况下，保证滤波精度，提高粒子滤波的实时性；陈世明等[10]提出一种基于引力场的粒子滤波算法以避免传统粒子滤波算法中粒子贫化与退化现象。

对粒子滤波的研究在一定程度上解决了粒子权值退化和粒子贫化问题，但增加了算法运算量，影响算法实时性。融入新的智能优化算法，在解决粒子权值退化和粒子贫化时，有效降低运算量有着重要的理论意义和应用价值。因此，提出了基于布谷鸟算法改进的粒子滤波(CS-PF)，并实验仿真得出重要参数的取值范围。

1 主要算法基本原理

1.1 布谷鸟优化算法

剑桥大学学者YANG Xinshe和DEB Suash通过对布谷鸟寻窝产卵的行为的模拟，提出了一种全局搜索的Cuckoo Search(CS)算法[11]。其基本思想基于3个理想假设。

(1)布谷鸟一次产一个卵，并随机选择鸟窝位置孵化。

(2)在随机选择的一组鸟窝中，最好的鸟窝会被保留到下一代。

(3)可利用的鸟窝数量n为常量，宿主发现布谷鸟蛋的概率为pα∈[0,1]。

基于以上3个理想条件，布谷鸟寻找孵化鸟窝的路径和位置更新公式如下

(1)

1.2 粒子滤波基本原理

粒子滤波基本原理是利用贝叶斯滤波理论解决状态估计问题，计算粒子权重，使用蒙特卡罗序列[12]方法确定状态的后验概率，其过程如下。

已知状态的初始概率密度函数为p(x0|z0)=p(x0)，状态更新方程为

(2)

(3)

(4)

然而，在实际应用时存在抽取有效样本困难的问题，为此，引入了重要性采样方法[12]，主要计算公式如下

(5)

2 布谷鸟算法优化的粒子滤波(CS-PF)

2.1 改进粒子的目标函数

通过增加扰动项来进一步保证粒子的多样性。从式(1)分析，第i个鸟窝在第t+1代的鸟窝位置与第t代位置加上Levy随机搜索路径，由于Levy飞行是一个非高斯的随机过程[13]。理论上，如果将其应用于粒子进化过程，短步长搜索有利于粒子向高似然区移动，偶尔大步长搜索，避免粒子陷入局部最优，且能保证粒子向下迭代时的多样性。但是，仍然存在局部最优的可能性。因此，增加扰动项对式(1)进行优化，如下所示

(6)

2.2 粒子权值更新

将布谷鸟算法融入PF的核心思想是，将粒子滤波中的每个粒子都进行迭代寻优，从而，粒子向后验概率密度值高的区域移动，进而提高估计结果的准确性。但是，各粒子在状态空间的位置也发生了改变，各粒子所表示的分布密度函数p(xk|y1:k-1)被改变，那么粒子滤波基于贝叶斯的理论基础就丢失。所以，在粒子更新时，对粒子权值进行更新。参照文献[3]，利用式(6)进行粒子位置更新时更新粒子权值。公式如下

(7)

2.3 CS-PF算法和步骤

基于2.1和2.2节的改进，提出的CS-PF算法流程图，如图1所示。

图1 CS-PF算法流程

CS-PF算法的详细过程如下。

步骤1 初始化N个粒子

步骤2 为了让粒子尽量向真实区域移动，又要避免最终收敛，这里利用小规模多种群思想，将k时刻N个粒子平均划分为m个子集，每个子集粒子为

FOR-1:m=1,2,…,log2N

FOR-2:j=1,2,…,N/log2N

步骤3 模拟布谷鸟鸟窝寻优

(1)生成随机数r∈[0,1]，若r>pα，执行步骤(2)，否则，执行步骤(3)。

(2)根据式(6)，对第m子集第k时刻第j个粒子的位置更新公式重新表示，如下所示

(8)

参数∂(通常∂=1)和运算符⊕与式(1)相同，L(λ)为Levy随机搜索路径，且L～u=s-λ(1<λ≤3)，s为Levy飞行随机步长，为了便于计算，通常使用下列公式计算Levy随机数

(9)

u,v为服从正态分布的随机数，u～N(0,σ2u),v～N(0,σ2v)， 0<β<2，通常β=1.5，σu、σv的取值如下

(10)

Γ为伽玛函数。返回步骤(1)。

FOR-2end

FOR-1end

步骤4 归一化粒子权值

(11)

步骤5 输出状态

(12)

由于布谷鸟算法具有全局最优特点，为了让粒子群整体向真实值附近移动，且避免最终收敛，我们通过设置随机数r是否小于阈值pα作为迭代的终止条件。因为通过粒子位置迭代更新的目的，并不是让所有粒子都集中于高似然区，只是为了让粒子集在高似然区合理分布，否则就会降低粒子的多样性。但是，pα的取值必须合适，如果迭代次数过多，增加算法时间复杂度，CS-PF的实时性能会受到很大影响；如果迭代次数过少，粒子集发散，不能合理集中于高似然区，就不能提高滤波精度。因此，在下一节的主要内容之一就是通过实验仿真的方法确定阈值pα的取值范围。

3 实验仿真

通过实验仿真的方法，确定重要参数pa的取值范围，并验证算法的有效性。

3.1 实验内容

利用以上实验环境和第2节提出的CS-PF进行以下3项实验，以充分说明算法的有效性。

(1)确定参数pα。在CS-PF算法中，是否更新粒子位置和算法运算量，主要取决于是否r>pα，r为随机数，r∈[0,1]，pα∈[0,1]，从理论上分析pα越小，通过式(8)粒子位置更新算法更新的粒子越多，较多粒子分布于高似然区，从而滤波精度较高，但也会较大的增加算法运算量；若pα越大，被更新位子的粒子越少，也就只有较小的粒子分布于高似然区，粒子保持了多样性，但滤波精度降低，算法运算量也小，实时性较高。因此，需通过实验仿真的方法，在滤波精度和算法实时性之间找到pα合适的取值范围。

(2)检验在高斯，非高斯噪声下CS-PF的性能。分别在量测噪声符合高斯分布、均匀分布、伽马分布、瑞利分布的噪声序列对CS-PF重复25次实验，对本文算法运算性能进行分析，并计算结果的均方根误差均值，计算公式如下

(13)

(3)CS-PF与其它算法对比分析。本文算法与解决非线性问题的经典算法，如文献[3]的基于萤火虫算法智能优化的粒子滤波(FA-PF)、文献[14]的基于自适应差分进化的粒子滤波(ADE-PF)、文献[15]的imp-WOPF和标准PF，在相同的过程方差和量测方差下进行滤波效果对比分析。

以上3项仿真实验基于以下过程状态方程和观测方程

(14)

(15)

式(14)为状态方程，式(15)为观测方程，式中：Q为制造过程噪声方差，R为量测噪声方差，rand为高斯分布的随机数。由于此系统为强非线性系统，且似然函数为双峰状[3]，一般滤波方法很难对其进行有效处理。参考现有文献，式(8)中，参数∂=1，步长因子α=0.1。式(9)中β=1.5。

3.2 实验结论

(1)参数pα取值范围。初始状态x=1，过程噪声Q和量测噪声R的方差均为5，模拟时序k=1,2,…,600，粒子数N=100。为了取得合适的pα，取pα=0.05+0.05i,(i=0,1,2,…,11)，并记录下pα等于不同值时，被更新的粒子平均数和粒子的运动趋势，结果如下。

(a)pα=0.05。即粒子在向下迭代时，理论上，r为均匀分布的随机数，pα越小，被更新的粒子越多，算法运算量越大，通过多次迭代，多数粒子将会向高似然区域移动，故不能保证粒子的多样性。算法滤波结果如图2所示，k=232时的粒子分布情况如图3所示。

图2 pα=0.05滤波效果

图3 pα= 0.05,k=232粒子分布情况

在图2中，平均每个粒子位置更新19.02次。RMSE=2.12。图3中，滤波进入后续阶段，PF算法的粒子大多数集中于少数状态值，从而失去了粒子多样性，不利于滤波的状态估计。而CS-PF在高似然区和低似然区都有一定数量的粒子分布，可以在一定程度上保证粒子的多样性。但是，在当前状态下，算法运算量较高。

(b)pα=0.10。平均每个粒子位置更新9.23次，RMSE=2.29，如图4所示，在k=232时，粒子分布情况，如图5所示。

图4 pα= 0.10滤波效果

图5 pα=0.10,k=232粒子分布情况

(c)pα=0.15。平均每个粒子位置更新6.04次，RMSE=1.67，如图6所示。在k=232时，粒子分布情况，如图7所示。

图6 pα=0.15滤波效果

图7 pα=0.15,k=232粒子分布情况

从以上3步实验可以初步进行以下推测。在pα逐渐增大时，CS-PF滤波精度略有变化，RMSE逐渐增大，但变化较小，而分布于高似然区的粒子数在逐渐减少，低似然区的粒子数在增加，从实验环节说明pα的增大有利于促进粒子的多样性，而其减小对算法滤波精度的降低有限，却较大增加了算法运算量。因此，参数pα的取值是影响算法精度和算法运算量的重要因素，同时，也极大地影响到粒子的多样性。

为了进一步充分说明pα对算法性能和粒子多样的影响，以及选取pα合适的取值范围，采用同样的方法重复以上实验直到pα= 0.60，经统计得出以下结论，见表1。

合适的pα∈[0.20,0.35]。从以上实验得出以下结论，算法运算量、滤波精度都和pα有一定的相关性。若pα<0.2，平均21.59个粒子在高似然区分布，平均每个粒子更新达20.66次，但同时也丧失了粒子多样性，且算法运算量大增，若粒子数较大，就极大影响了算法实时性；反之，若pα>0.5，绝大多数粒子位置都没有更新，平均更新次数在0.65左右，也就是每个粒子更新不到一次，分布于高似然区的粒子数就相对少很多，不能有效提高算法滤波精度，但算法运算量较低。因此，结合表1中高似然区平均粒子数的均值为10，综合考虑，既要保证足够多的粒子分布于高似然区，以保证滤波精度，又能以较低的运算量保证算法实时性，所以最为合适的pα∈[0.20,0.35]之间。

表1 pα取值、粒子位置更新情况与RMSE

(2)利用MatLab R2012a的gamrnd函数、rand函数、raylrnd函数和randn函数，生成伽马分布、均匀分布、瑞利分布和高斯分布的量测噪声序列。取pα=0.25，设初始值为x=1，噪声序列服从均值为0，过程噪声方差Q=5，量测噪声方差R=5，粒子数为N=100，滤波时序k=600，分析噪声在高斯分布、均匀分布、伽马分布、瑞利分布下，对算法的性能影响。针对不同噪声分布分别做25次独立实验，计算不同噪声分布下的RMSE均值。结果见表2。

从表2可以得出以下结论，过程噪声与量测噪声不变，高斯分布、均匀分布、伽马分布和瑞力分布的均方根误差分别为2.1722、2.1963、2.2491和2.2687，说明CS-PF针对强非线性系统，在固定的过程方差和量测方差、高斯分布和非高斯分布的量测噪声下滤波效果几乎一致。

表2 CS-PF在不同分布下的RMSE均值

(3)取pα=0.25，设初始值为x=1，噪声序列服从均值为0，过程噪声方差Q=5，量测噪声方差R=5，粒子数为N=100，滤波时序k=600，噪声分布为高斯分布。对文献[3]的基于萤火虫算法智能优化的粒子滤波(PA-PF)，文献[14]的ADE-PF和文献[15]的imp-WOPF对比研究。针对每种算法进行了25次滤波，并计算每种算法滤波的均方根误差均值，计算方法和式(13)相同，并统计分析每种滤波各自的运行时间均值。

程序运行环境：处理器，Intel(R) Core(TM) i5-4200U CPU @ 1.60 GHz 2.30 GHz；内存，4.00 G；操作系统，64位Windows 7。实验结果见表3。

表3 不同算法RMSE与运算量对比

在表3中，本文算法CS-PF的RMSE=2.1735，滤波精度低于文献[15]的imp-WOPF算法，却优于标准PF，以及FA-PF和ADE-PF，但运行时间高于标准PF，而少于其它几种滤波算法。主要是因为CS-PF中融入了基于布谷鸟算法的粒子寻优过程，但是，该过程的运算量却低于萤火虫算发、自适应差分进化的粒子滤波和权值改进算法。

综合上述实验结果，得出以下结论。

(1)将布谷鸟算法的路径寻优思想融入粒子滤波代替粒子重采样，能够解决粒子权值退化问题，并保持粒子多样性。

(2)相比标准粒子滤波，在运算效率降低12%的情况下，通过合理的粒子更新阈值pα=0.25，能够提高滤波精度18%。

3.3 实际应用

为了进一步验证算法的可靠性，在纺织实际生产环境中对CS-PF算法进行验证。计算机软硬件环境为一台Dell服务器(处理器为2个Intel(R) Xeon(R) CPU E5-2407 2.40 GHz，Windows Server 2008(64位)操作系统，32 GB内存，8 TB硬盘容量，通信带宽：1 Gb/s)，对织布车间300台织机的原始数据进行实时在线监测。通过实际应用表明：在针对非平稳、非线性的织造过程数据时，具有很好的适用性，结果误差小于3%，误差范围小于生产环境5%的要求。结果表明：本文算法能够有效地处理织机声非线性织造过程的非平稳数据，有效地保证织造过程中采集坯布的质量和产量数据的准确性，从而，很好地满足生产管理的需要，以及生产管理、统计分析数据的准确性也大大提高。

4 结束语

基于布谷鸟路径寻优思想融入粒子滤波，以解决粒子滤波中的权值退化和粒子贫化问题。实验结果表明：CS-PF在利用布谷鸟寻优思想代替粒子重采样后，能够保持粒子多样性和滤波精度。与其它改进的粒子滤波算法相比，CS-PF具有以下有点：①粒子位置的迭代更新采用Levy飞行模式，通过阈值控制迭代次数，能够保持粒子多样性，克服陷入局部最优。使得部分粒子分布于高似然区，部分分布于低似然区。②CS-PF调节参数，易于编程实现，容易应用于实际工程问题。③在运算效率降低较少的情况下，运算精度有了较大提高。未来主要的研究应着眼于，结合现有研究基础，在实际工程问题中自适应的控制粒子位置更新阈值，以及粒子数的自适应选择上，以进一步降低运行时间，提高算法的实时性。