粒子滤波算法(PF)可以有效处理非线性、非高斯噪声环境下的状态估计问题，被广泛用于目标跟踪[1-2]、故障诊断[3-4]、机器人定位[5-6]等领域．PF算法利用大权值粒子取代小权值粒子来完成重采样，导致该算法存在粒子退化和样本贫化等现象，降低了算法的精确度．针对该问题，国内外许多学者对PF算法进行了改进，主要分为两类：一类是对算法的建议分布函数进行改进[7-10]，该类方法通过改进粒子滤波算法的建议分布函数提高算法精度，但实际上有效的建议分布函数很难获取，现有改进的建议分布函数对算法的优化效果不佳；另一类是对算法的重采样方法进行改进[11-13]，该类方法是目前较为流行的方法，通过对算法的重采样环节进行改进提高算法的精确度，可以取得较好的效果，但同时也增加了算法的时间复杂度．在对PF算法进行改进的现有研究中，虽然算法的精确度得到了提高，但是大部分方法都增加了算法的时间复杂度，导致算法实时性较差，难以将其应用到实时性要求较高的领域．图形处理器(GPU)具有多个计算单元和多种内存，适合处理计算量较大且计算简单的任务，众多研究学者已将GPU的多线程技术应用到雷达信号处理[14-15]、遥感卫星图像处理[16-17]等实时性要求较高的领域．很多研究者将GPU并行计算技术引入粒子滤波算法，提高算法的实时性．文献[18]提出了一种粒子滤波的并行系统重采样算法，并行实现重采样算法的循环迭代，解决粒子滤波算法重采样耗时的问题．文献[19]针对粒子滤波算法重采样不能充分并行的问题，分析了不涉及集体操作的拒绝重采样方案．文献[20]针对粒子滤波算法在跟踪领域计算量较大的问题，提出了基于GPU的粒子滤波并行算法，提高计算速度．文献[21]针对实时系统中粒子滤波的计算耗时问题，提出了一种基于计算机统一设备架构(CUDA)的并行加速粒子滤波算法，通过以线程块为单位并行处理粒子滤波算法，提高算法实时性．以上文献通过GPU并行计算技术实现粒子滤波算法，提高了算法的实时性，但粒子滤波算法在重采样环节因粒子之间频繁交互，导致算法不能充分完成并行化处理，降低了算法的并行效率．针对以上问题，本研究提出了一种基于GPU的骨干粒子群优化粒子滤波算法，该方法首先通过骨干粒子群算法对粒子滤波算法的重采样环节进行优化，提高算法的并行性．然后从GPU的多线程架构和内存访问原则两个方面对算法进行并行优化，使粒子滤波算法得到充分的并行化处理，从而进一步提高算法的实时性．1 算法原理1.1　骨干粒子群算法骨干粒子群算法(BBPSO)是一种改进的粒子群算法．该算法利用全局最优值和局部最优值进行高斯采样来完成粒子位置的更新，无须粒子交互．与经典的粒子群算法相比，BBPSO算法消除了惯性系数、数度阈值和加速系数等控制参数，使得算法结构更为简洁且易于操作．粒子的位置更新方程为Xijt+1=N(uij,σij)；uij(t)=0.5(pij(t)+Gj(t))；σij(t)=|pij(t)-Gj(t)|，式中：Xijt+1为t+1时刻粒子i位置的第j维分量，i=1，2，…，I，j=1，2，…，D，其中，I为种群规模，D为带求解问题的维数；uij(t)为粒子自身最优值与全局最优值的均值；σij(t)为粒子自身最优值与全局最优值的标准差；N(uij,σij)为高斯分布；pi为粒子i自身的最优位置；G为群体最优位置．可见粒子的位置更新是通过局部最优值与全局最优值的算术平均值uij和标准差σij进行高斯采样得到的．局部最优值和群体最优值更新方程为：pit+1=Xi(t+1)    (f(Xi(t+1))f(pi(t))),pi(t)    (其他);Gt+1=pi(t+1)    (f(pi(t+1))f(G(t))),G(t)    (其他),式中f(t)为适应度函数．骨干粒子群算法结构简单，每个粒子种群分别利用局部最优值与全局最优值联合进行粒子更新．粒子种群相互独立运行，适合并行实现．利用该算法取代PF算法重采样环节，提出了并行骨干粒子群优化的粒子滤波算法．1.2　骨干粒子群算法优化粒子滤波算法骨干粒子群优化粒子滤波算法流程如图1所示，图中FX为PF算法的状态更新方程．首先设置初始状态，通过PF算法的状态方程进行状态更新，然后利用测量方程初始化该时刻的状态估计粒子，并将初始化的粒子传递给BBPSO算法进行处理，根据提前设置好的迭代结束条件结束重采样，骨干粒子群算法将状态估计结果传送到PF算法中进行状态的更新．10.13245/j.hust.210303.F001图1骨干粒子群优化粒子滤波算法流程2 BBPSO-PF算法的并行优化设计在CUDA框架下构建并行骨干粒子群优化粒子滤波算法，通过骨干粒子群算法取代粒子滤波重采样环节，算法流程如图2所示．10.13245/j.hust.210303.F002图2并行BBPSO-PF流程图粒子滤波算法完成状态更新后将初始化的粒子分成多个粒子群，然后利用适应度函数计算每个粒子群的适应度值，根据适应度值更新每个粒子群的个体最优值，并根据整体适应度值更新全局最优值，最后通过局部最优值和全局最优值构造高斯函数进行高斯采样更新粒子．从算法流程中可以看出：BBPSO算法中粒子种群相互独立，具有并行结构的特点，利用该算法取代PF算法重采样环节，可以克服PF算法重采样环节因粒子交互而不能充分并行的缺点．GPU并行BBPSO-PF算法以粒子群为单位，如果每个线程处理一个粒子的数据，那么将导致多个线程处于闲置状态，无法同时工作，从而浪费了大量计算资源．针对以上问题，从多线程和内存访问两个方面对并行的BBPSO-PF算法进行优化．2.1　基于多线程优化BBPSO算法在求粒子适应度、个体最优值、粒子更新时都是以粒子群为单位，相互独立运行．如果每个线程负责一个粒子的数据处理，如图3所示，计算每个粒子群的适应度值时，线程之间须要进行累加运算，导致多个线程处于闲置状态，浪费大量计算资源．因此，利用GPU的多线程优势对该算法进行并行优化设计，模型如图4所示，每个粒子群包含多个粒子(X表示粒子)和一个粒子群的适应度值(y表示适应度值)，当须要计算的粒子个数超过可用线程数量时，每个线程负责一个粒子群的数据处理(计算粒子适应度、寻找个体最优值、粒子更新)，可使每个线程都处于工作状态，让算法在粒子群级别得到并行实现，从而提高算法实时性．10.13245/j.hust.210303.F003图3每个线程负责一个粒子10.13245/j.hust.210303.F004图4多线程处理多粒子群2.2　基于内存访问优化GPU内存的访问遵循合并与对齐原则，在GPU端系统的执行单元为一个线程束(32个线程)，当线程束的一个内存访问事务的起始地址是缓存粒度的偶数倍(32 byte的二级缓存或128 byte的一级缓存)，就会出现对齐内存访问，当一个线程束的32个线程访问一段连续的内存时，就出现合并内存访问，线程束的对齐与合并内存访问模型如图5所示．在GPU物理内存上，数据都是按照一维方式存放，所以粒子群的数据结构在GPU的全局内存上存放形式如图6所示，每个线程负责一个粒子群的数据处理，其中包含3个粒子和1个适应度值(通过实验发现将每个粒子群的个数设置为3个时并行效果最佳，所以实验中将每个粒子群的维度设定为3)，可以看出：数据的存取不符合对齐与合并原则，损耗多次内存访问事务．10.13245/j.hust.210303.F005图5内存访问模型10.13245/j.hust.210303.F006图6粒子群数据存放形式针对以上问题利用GPU内存的对齐与合并访问原则，从内存的访问方面进一步对算法进行并行优化．为了减少内存访问事务，重新设计粒子群的数据结构，可将适应度值与粒子分开存放，如图7所示，将适应度值存放在连续的一段内存中，每个线程束只需要一次内存访问事务就可以完成粒子群的适应度值的访问．由于粒子群中的粒子都是随机组合，因此可以进一步改进粒子群数据的存放方式，如图8所示，粒子群中的每个粒子按照固定间隔存放，这样一个线程束访问的是连续的内存空间，满足了内存访问的对齐与合并原则，极大节省了内存访问事务．10.13245/j.hust.210303.F007图7粒子与适应度值分开存放形式10.13245/j.hust.210303.F008图8粒子群的合并与对齐存放形式3 实验仿真与分析在Win7系统下采用VisualStu-dio2013(CUDA7.5)开发环境进行实验，CPU为Intel(R)Core(TM)i3-4150，GPU为NVIDA Ge-Force GTX 950．利用CUDA架构并行实现BBPSO-PF算法和文献[18]、文献[19]的重采样方法，对比算法并行运行时间．采用一个广泛应用的标量模型，其状态模型和观测模型为Xk=0.5Xk-1+2.5Xk-1/1+Xk-12+8cos1.2k+Wk-1;Zk=Xk2/20+Vk，式中：Wk为模型k时刻的状态噪声；Vk为观测噪声．实验参数设置：初始状态为X0=0.1；状态噪声为Wk～N0,0.1；观测噪声为Vk～N0,0.1；仿真周期为T=50．3.1　算法在CPU端执行时间对比利用BBPSO算法对PF算法的重采样环节进行优化，在保证算法精度的前提下提高算法实时性．通过与随机重采样(random resampling，RAR)、多项式重采样(multinomial resampling，MR)、系统重采样(systematic resampling，SR)和残差重采样(residual resampling，RR)方法的状态估计精度进行对比，状态估计误差如表1所示．从表中数据可以看出：BBPSO-PF算法的状态估计误差低于其他重采样方法，表明BBPSO算法优化PF算法提高了算法的精度．在不同粒子数的情况下，对BBPSO-PF算法和文献[18]、文献[19]的重采样方法进行串行时间的测试，结果如图9所示，可以看出：BBPSO-PF算法的实时性介于文献[18]和文献[19]的算法，原因在于BBPSO算法与PF算法的结合增加了算法的复杂度．10.13245/j.hust.210303.T001表1平均状态估计误差粒子数RARMRSRRRBBPSO1 5000.2720.2780.2840.3150.2033 0000.2670.2740.2740.2920.1964 5000.2500.2650.2670.2660.1926 0000.2500.2540.2610.2130.17210.13245/j.hust.210303.F009图9算法串行执行时间3.2　基于多线程优化的实验效果在多种粒子数的情况下对提出的多线程并行BBPSO-PF算法和文献[18]、文献[19]进行对比，实验结果如图10所示．可以看出：基于多线程的并行BBPSO-PF算法实时性优于文献[18]和文献[19]的重采样方法，可见利用骨干粒子群算法优化粒子滤波重采样环节，并通过多线程优化方案对并行的BBPSO-PF算法进行优化的改进思路使本文算法的实时性得到了有效提高，具体实验数据如表2所示，文献[18]、文献[19]中的优化方案与本文多线程优化方案之间的加速比1优化(文献[18]/多线程优化)和加速比2优化(文献[19]/多线程优化)都超过了1倍，说明提出的多线程优化方案可以有效缩短算法的运行时间，从而验证了多线程优化方案对提高算法实时性具有重要的作用．10.13245/j.hust.210303.F010图10多线程优化后的并行执行时间10.13245/j.hust.210303.T002表2多线程优化时间对比粒子数文献[18]文献[19]多线程优化加速比1优化加速比2优化1 9203473893041.141 41.279 63 8404746024181.134 01.440 24 8005166964301.200 01.618 65 7605407864471.208 11.758 46 7206029194531.328 92.028 77 6806199645631.099 51.712 39 6006961 0956211.120 81.763 3ms3.3　基于内存访问优化的实验效果根据BBPSO-PF算法中粒子群的具体数据结构设计高效内存访问方案．实验效果如图11所示，可以看出：在多线程优化方案的基础上进行内存访问优化，可以进一步降低算法的运行时间，提高并行BBPSO-PF算法的实时性．具体实验数据如表3所示，通过对比多线程优化方案与内存访问优化方案的加速比(多线程优化/内存访问优化)可以看出：内存访问优化可以在多线程优化的基础上进一步提高算法的实时性，充分验证了内存访问优化的必要性．10.13245/j.hust.210303.F011图11内存优化后的并行执行时间10.13245/j.hust.210303.T003表3内存访问优化时间对比粒子数多线程优化内存访问优化加速比优化1 9203042891.051 93 8404183741.117 64 8004303901.102 65 7604474061.101 06 7204534091.107 67 6805634591.226 69 6006215151.205 8ms4 结论针对粒子滤波算法在重采样环节因粒子交互而不能充分并行的缺点，提出基于GPU的BBPSO-PF并行算法．针对粒子滤波算法重采样环节耗时且不易并行的问题，充分利用骨干粒子群算法具有并行结构的特点优化粒子滤波算法重采样环节，提高算法的并行度．通过对BBPSO-PF算法的并行性进行分析，利用GPU的多线程结构对并行的BBPSO-PF算法进行优化，使粒子群之间得到并行化处理，从而提高算法的实时性．在多线程优化的基础上通过GPU合并与对齐的内存访问原则，设计粒子群的数据结构，使多线程以最少的内存访问事务存取数据，从内存存取方面进一步提高并行BBPSO-PF算法的实时性．为了验证并行优化后的算法实时性，并行实现本文算法，统计算法运行时间，并与部分现有文献进行对比，实验结果表明：在算法实时性方面，本文方法相较于对比实验分别提高了1.47和2.24倍。