基于分类进化重采样的粒子滤波方法与流程

本发明涉及粒子滤波
技术领域：
，特别是涉及一种基于分类进化重采样的粒子滤波方法。
背景技术：
：粒子滤波(particlefilter,pf)是一种基于蒙特卡罗方法(montecarlomethods)仿真的递推贝叶斯估计方法，它是在序贯重要性采样(sequentialimportancesampling,sis)的基础上，再进行重采样的一个递推过程。粒子滤波方法可以有效地处理非线性、非高斯噪声环境下的状态估计问题，被广泛用于在目标跟踪、导航定位等领域。自从1993年戈登(gordon)等提出重采样方法(samplingimportanceresembling,sir)以来，出现了很多重采样方法，比较经典的有残差重采样(residualresampling,rr)、多项式重采样(multinomialresampling,mr)和系统重采样(systematicresampling,sr)。虽然这些重采样方法实现的手段不同，但采样的思想相同，都是通过丢弃小权值粒子，对大权值粒子重新分配权值获取新的粒子集合。虽然在一定程度上缓解了样本退化的问题，但随着时间的延长，大权值粒子多次被复制，粒子集的多样性丧失。经过若干次迭代以后，用于估计的粒子集合只是部分有限状态的子代粒子，严重丧失了状态的遍历性能，特别是在系统状态分布具有严重拖尾现象的时候，迭代时间越长、精度越差。针对这一问题，研究人员提出了不少改进方法，总结起来主要分为两类：第一类是基于“优选建议分布函数”的思想，通过滤波技术将最新的观测信息融入到建议分布函数中，增强当前时刻观测信息对粒子多样性的修正作用。例如，王法胜等采用扩展卡尔曼滤波(extendedkalmanfilter,ekf)对最新观测信息进行滤波优化，生成建议分布函数，改善了粒子退化问题，但是ekf的一阶线性近似性却引入了额外的近似误差，同样会产生误差累积效应。julier等的无迹卡尔曼滤波(unscentedkalmanfilter,ukf)估计生成建议分布函数，可以达到2阶以上的近似能力，进一步提升了滤波精度，降低了粒子滤波权值退化问题，但是复杂的ut(unscentedtransformation)变换却导致算法的实时性太差。wucl等采用具有三阶精度的容积卡尔曼滤波生成建议分布函数，由于容积卡尔曼滤波的多数计算均可采用离线形式，提升滤波精度的同时，并没由增加计算负担，虽然这类基于建议分布函优化的方法一定程度上增强了最新观测信息对粒子多样性的保持作用，但是改善的精度严重依赖观测信息的估计方法，没有改善重采样的实现手段，无法避免权值衰退的产生。第二类是基于“优化重采样策略”的思想，通过智能寻优的方法改变重采样手段。grisetti等利用网格分割的思想权重的重新分配，提出了确定性重采样方法，由于网格整体风格的计算量过大，限制了该技术的进一步发展。cheng等通过设置权重阈值，将权值粒子集分为两个采样集合，保留大权值粒子，并结合小权值进行分层采样，每一次采样后都会摒弃大量的小权值粒子，降低了多样性；在此基础上，zhang等提出采用遗传算法对剩余粒子进行优化分析，针对大小权值粒子采取不同的遗传变异修正处理方式实现了遗传算法重采样方法(geneticalgorithmresampling,gar)，提升了粒子滤波在机动目标跟踪中的滤波精度，但是每一次重采样都要对粒子集合的全体粒子操作，增加了算法本身的计算时间。如何在不增加计算时间的条件下避免粒子滤波方法的粒子滤波权值退化、样本贫化，提高粒子滤波方法的精确度和稳定性是本领域亟待解决的问题。技术实现要素：基于此，有必要针对如何在不增加计算时间的条件下避免粒子滤波方法的粒子滤波权值退化、样本贫化，提高粒子滤波方法的精确度和稳定性的问题，提供一种基于分类进化重采样(classificationevolutionresampling，cer)的粒子滤波方法。本发明提供的一种基于分类进化重采样的粒子滤波方法，包括以下步骤：采样步骤，对建议密度函数进行采样，生成当前粒子集；权重赋值步骤，计算所述当前粒子集的每个粒子的权值；判断步骤，判断所述当前粒子集的有效粒子数neff是否小于预设有效粒子阀值nthr，当所述有效粒子数neff小于所述预设有效粒子阀值nthr时，执行重采样步骤；重采样步骤，根据所述当前粒子集的每个粒子的权值大小，将所述当前粒子集分类为小权值粒子种群、保留粒子种群和大权值粒子种群，将所述大权值粒子种群进行裂变获得大权值裂变粒子种群，将所述大权值裂变粒子种群与所述小权值粒子种群进行差分进化，获得最优粒子集；优化步骤，根据所述最优粒子集及相应的每个粒子的权值优化建议密度函数。在其中的一个实施例中，在所述重采样步骤中，将所述当前粒子集分类为小权值粒子种群、保留粒子种群和大权值粒子种群，包括：比较所述当前粒子集的每个粒子的权值与第一权值阀值ωthrmin以及第二权值阀值ωthrmax的大小，其中，所述第一权值阀值ωthrmin小于所述第二权值阀值ωthrmax；当所述当前粒子集的粒子的权值小于所述第一权值阀值ωthrmin时，所述当前粒子集的该粒子归为小权值粒子种群；当所述当前粒子集的粒子的权值大于所述第二权值阀值ωthrmax时，所述当前粒子集的该粒子归为大权值粒子种群；当所述当前粒子集的粒子的权值介于所述第一权值阀值ωthrmin和所述第二权值阀值ωthrmax之间时，所述当前粒子集的该粒子归为保留粒子种群。在其中的一个实施例中，在所述重采样步骤中，所述大权值裂变粒子种群的裂变个数为其中，n为小权值粒子种群的粒子个数，p为大权值粒子种群的粒子个数，k为保留粒子种群的粒子个数。在其中的一个实施例中，在所述重采样步骤中，所述大权值裂变粒子种群为所述大权值裂变粒子种群的粒子满足所述大权值裂变粒子种群的当前粒子的权值为所述大权值裂变粒子种群的新一代粒子的权值为所述大权值裂变粒子种群的粒子分布为其中，表示大权值裂变粒子种群的当前粒子，表示大权值裂变粒子种群的新一代粒子，且满足在其中的一个实施例中，在所述重采样步骤中，将所述大权值裂变粒子种群与所述小权值粒子种群进行差分进化，获得最优粒子集，包括：突变，从所述小权值粒子种群和所述大权值裂变粒子种群分别随机选取一个粒子，按照进行迭代，获得突变粒子种群，其中，是在第k次迭代中需要干扰的目标矢量，和是两个随机粒子，分别来所述小权值粒子种群和所述大权值裂变粒子种群。i1和i3分别为所述小权值粒子种群和所述大权值裂变粒子种群的索引，参量β控制差分变化的缩放比例；交叉，将所述突变粒子种群按照进行交叉，获的交叉粒子种群，其中，d表示维数，j为交叉点集合，其中的元素从集合{1,2,...,d}中随机的选择，表示矢量的第j个元素；选择，比较所述交叉粒子种群与所述大权值粒子种群的拟合效果的精确度，选择精确度高的粒子种群作为新一代粒子种群进行突变、交叉、选择迭代，直至粒子种群达到最大数目，获得最优粒子集。在其中的一个实施例中，评价所述拟合效果的拟合函数为其中，在其中的一个实施例中，在所述采样步骤中，所述建议密度函数为在其中的一个实施例中，在所述权重赋值步骤中，所述当前粒子集的每个粒子的权值为：在其中的一个实施例中，在所述判断步骤中，所述有效粒子数neff为上述基于分类进化重采样的粒子滤波方法，在重采样过程中，首先对大权值粒子种群进行裂变，然后采用差分进化方法同小权值种群进行优化，通过突变、交叉、选择产生新粒子种群，避免了对历史小权值信息的丢弃，避免了粒子滤波权值退化，保证了样本的多样性。进一步地，该滤波方法仅对满足重采样条件的建议密度函数进行滤波估计，在不增加计算负担的情况下保证了最新观测信息对粒子的修正作用，提高了粒子滤波方法的精确度和稳定性。附图说明为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲，还可以根据这些附图获得其他的附图。图1为本发明基于分类进化重采样的粒子滤波方法一实施例的重采样粒子分类示意图；图2为gfa-pf采用本发明的重采样方法在样本数为100、迭代100次时样本分布图；图3为pf采用本发明的重采样方法样本数为100、迭代100次时样本分布图；图4为不同重采样方法在采样数为300时的状态估计曲线；图5为不同重采样方法在采样数为300时的均方根误差rmse曲线；图6为不同重采样方法在采样数为100时的状态估计曲线；图7为不同重采样方法在采样数为100时的均方根误差rmse曲线；图8为不同重采样方法在采样数为20时的状态估计曲线；图9为不同重采样方法在采样数为20时的均方根误差rmse曲线。具体实施方式为了使本发明的目的、技术方案及优点更加清楚明白，以下通过实施例，并结合附图，对本发明的基于分类进化重采样的粒子滤波方法进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。在本发明一实施例的基于分类进化重采样的粒子滤波方法中，以非线性系统粒子滤波为例，建立如式(1)所示的系统过程方程和观测方程：其中，xt∈rn,yt∈rm为状态向量和观测向量，f(.)和h(.)是已知的非线性传递函数，ωt和vt是过程和量测噪声，为不相关的零均值高斯白噪声过程，协方差为q、r。滤波的目的是在给定观测信息p(xt|y1:t)的情况下计算状态的后验概率密度(pdf)。随着观测信息的不断更新，建立如式(2)所示的系统状态的后验概率密度以及如式(3)所示的条件概率密度：p(yt|y1：t-1)＝∫p(yt|xt)p(xt|y1：t-1)dxt(3)递归贝叶斯滤波解决非高斯系统的难点就是积分的求解，粒子滤波基于蒙特卡洛仿真的思想采用如式(4)所示的权重粒子加和的方式近似计算积分：其中，是第m个粒子的位置和权值，δ(.)是dirac函数，n是粒子数目。本发明基于分类进化重采样的粒子滤波方法，包括以下步骤：采样步骤，对建议密度函数进行采样，生成当前粒子集；权重赋值步骤，计算所述当前粒子集的每个粒子的权值；判断步骤，判断所述当前粒子集的有效粒子数是否小于预设有效粒子阀值，当所述有效粒子数小于所述预设有效粒子阀值时，执行重采样步骤；重采样步骤，根据所述当前粒子集的每个粒子的权值大小，将所述当前粒子集分类为小权值粒子种群、保留粒子种群和大权值粒子种群，将所述大权值粒子种群进行裂变获得大权值裂变粒子种群，将所述大权值裂变粒子种群与所述小权值粒子种群进行差分进化，获得最优粒子集；优化步骤，根据所述最优粒子集及相应的每个粒子的权值优化建议密度函数。上述基于分类进化重采样的粒子滤波方法，在重采样过程中，首先对大权值粒子种群进行裂变，然后采用差分进化方法同小权值种群进行优化，通过突变、交叉、选择产生新粒子种群，避免了对历史小权值信息的丢弃，避免了粒子滤波权值退化，保证了样本的多样性。进一步地，该滤波方法仅对满足重采样条件的建议密度函数进行滤波估计，在不增加计算负担的情况下保证了最新观测信息对粒子的修正作用，提高了粒子滤波方法的精确度和稳定性。可选地，在所述采样步骤中，所述建议密度函数为目前常用的两个建议密度函数分别是状态转移概率和后者因为考虑了最新观测信息对粒子采样的修正作用，具有更好的精度，但是需要一定的滤波方法将观测信息融入到系统中，增加了粒子滤波的计算复杂度。在所述权重赋值步骤中，采用如式(5)所示计算当前粒子集的每个粒子的权值：随着迭代时间的增加，粒子权值的方差趋于零，出现“粒子退化”问题，需要进行重采样步骤，但是何时进行重采样以及如何进行重采样是一个关键性的问题。定义有效粒子数neff如式(6)所示：通过判断所述当前粒子集的有效粒子数neff是否小于预设有效粒子阀值nthr来确定是否进行重采样，即通过度量当前粒子集代表真实后验的好坏程度来决定重采样的时间，当有效粒子数neff低于预设有效粒子阀值nthr，即当前粒子集已退化，不能满足粒子滤波精度的情况下，进行重采样。本发明在判断何时进行重采样时，不是对全部的当前粒子集的粒子都进行ukf滤波，而是在判断需要进行重采样步骤以后，在进行重采样的同时，进行建议分布函数的滤波估计优化。判断阈值的设置为neff＜nthr，满足重采样条件后，ukf优化实现步骤保持不变，具体的可以参考arulampalamms,simonm,gordonn.的在线非线性/非高斯贝叶斯粒子滤波器的教程(atutorialonparticlefiltersforon-linenonlinearnon-gaussianbayesiantracking)。请参阅图1所示，本发明提出的基于分类进化重采样的粒子滤波方法的重采样方法主要包括“粒子分类”和“粒子进化”两个步骤。首先根据所述当前粒子集的每个粒子的权值大小，将所述当前粒子集分类为小权值粒子种群、保留粒子种群和大权值粒子种群，再将所述大权值粒子种群进行裂变获得大权值裂变粒子种群，将所述大权值裂变粒子种群与所述小权值粒子种群进行差分进化，获得最优粒子集。在本发明中，为避免对所有当前粒子集的粒子进行重采样而增加运算量，将当前粒子集进行分类处理。根据聚类思想，基于粒子的权值大小将整个当前粒子集划分为小权值粒子种群、保留粒子种群和大权值粒子种群，其中小权值粒子种群需要进行重采样，保留粒子种群在粒子更新中保持不便，大权值粒子种群裂变后与小权值粒子种群进行差分优化，以更好地保持样本多样性。可选地，将所述当前粒子集分类为小权值粒子种群、保留粒子种群和大权值粒子种群，包括：比较所述当前粒子集的每个粒子的权值与第一权值阀值ωthrmin以及第二权值阀值ωthrmax的大小，其中，所述第一权值阀值ωthrmin小于所述第二权值阀值ωthrmax；当所述当前粒子集的粒子的权值小于所述第一权值阀值ωthrmin时，所述当前粒子集的该粒子归为小权值粒子种群；当所述当前粒子集的粒子的权值大于所述第二权值阀值ωthrmax时，所述当前粒子集的该粒子归为大权值粒子种群；当所述当前粒子集的粒子的权值介于所述第一权值阀值ωthrmin和所述第二权值阀值ωthrmax之间时，所述当前粒子集的该粒子归为保留粒子种群。将当前粒子集分类为小权值粒子种群、保留粒子种群和大权值粒子种群的公式如式(7)所示：进一步可选地，大权值粒子种群的裂变个数主要取决于粒子权值，其计算方法如式(8)所示：其中，n为小权值粒子种群的粒子个数，p为大权值粒子种群的粒子个数，k为保留粒子种群的粒子个数。进一步可选地，所述大权值裂变粒子种群表示为所述大权值裂变粒子种群的当前粒子权值如式(9)所示：所述大权值裂变粒子种群的粒子满足如式(10)所示的条件：所述大权值裂变粒子种群的新一代粒子的权值为如式(11)所示：所述大权值裂变粒子种群的粒子分布为高斯分布，如式(12)所示：其中，表示大权值裂变粒子种群的当前粒子，表示大权值裂变粒子种群的新一代粒子，满足如式(13)所示的关系：其中，进一步地，与传统的重采样方法不同，在本发明的重采样步骤中，不存在粒子删除和复制步骤，而是将大权值粒子裂变集合与小权值粒子进行差分进化，以保持样本的多样性，避免粒子滤波权值退化。假设总的粒子种群目标矢量为相应权值集合为表示当前的差分进化(de)种群，由n个d维目标矢量构成。de种群获取的新矢量集合为u，t时刻，第k次迭代的矢量集合目标如式(14)所示：是在第k次迭代中需要干扰的目标矢量，和是两个随机粒子，分别来所述小权值粒子种群和所述大权值裂变粒子种群。i1和i3分别为所述小权值粒子种群和所述大权值裂变粒子种群的索引，参量β控制差分变化的缩放比例。进一步地，为增加新生矢量的多样性，采用交叉实现的方法为如式(15)、式(16)所示：将所述突变粒子种群按照进行交叉，获的交叉粒子种群，其中，d表示维数，j为交叉点集合，其中的元素从集合{1,2,...,d}中随机的选择，表示矢量的第j个元素；选择，比较所述交叉粒子种群与所述大权值粒子种群的拟合效果的精确度，选择精确度高的粒子种群作为新一代粒子种群进行突变、交叉、选择迭代，直至粒子种群达到最大数目，获得最优粒子集。可选地，评价所述拟合效果的拟合函数为其中，以下利用计算机仿真分析，将本发明的基于分类进化重采样的粒子滤波方法(cer)与rr、mr、sr以及gar进行了对比分析实验。具体地，对比分析采用式(19)所示的单变量增长模型进行仿真分析，该模型具有较强的非线性且噪声满足高斯噪声特性。其中，是相互独立的高斯噪声，α、β、γ是常数参量，在实施例中设置α＝1,β＝12,γ＝7，且σω＝σv＝1。系统初始状态真实值为x0＝1，估计值为对比分析实验评价指标为定量评价不同采样方法的性能指标，在对比分析实验中，设定了均方根误差(rmse)以及rmse的均值mermse和方差varrmse，可选地，rmse及其均值mermse、方差varrmse的计算表达式如式(20)以及式(21)所示：其中，m为采样粒子个数，t为仿真周期，xt和分别为系统状态的真实值和估计值。对比分析实验结果分析实验一：粒子样本空间分布多样性分析。在该实验中，用有效粒子数neff度量粒子滤波方法的退化程度，有效粒子数定义如式(6)所示。其中，neff的大小与退化程度成反比关系。在采样粒子数n＝100的情况下，分别对几种不同重采样方法的粒子滤波方法进行了100次的独立运行，其中有效粒子数neff实验结果如表1所示。进一步采用本发明cer的重采样方法是否能够改善基于引力场的粒子滤波算法(gfa-pf)以及pf两种粒子滤波方法的粒子退化问题。表1不同迭代次数的由效粒子数从表1中可以看出，本发明方法的平均有效粒子数保持了最大值，因此，本发明方法能够很好地解决粒子退化问题。请参阅图2所示，gfa-pf采用本发明cer的重采样方法在样本数为100、迭代100次时样本分布图，从图中可以看出，粒子分布均匀，多样性丰富。请继续参阅图3所示，pf采用本发明的重采样方法样本数为100、迭代100次时样本分布图，粒子出现了明显的聚集现象，单一性强；对比图2和图3可知，gfa-pf方法采用本发明的重采样方法在解决粒子退化问题的同时还能保证粒子的多样性，解决粒子贫化问题。而pf方法采用本发明的重采样方法虽然可以在一定程度上避免粒子退化问题，但是却带来了粒子贫化问题。实验二：固定仿真周期t＝50，分析不同采样粒子数情况下的估计效果。为分析不同采样方法的精确性，首先对固定仿真周期、不同采样粒子数情况下的估计精度进行了对比分析实验。请参阅图4至图9所示，分别给出了粒子采样数分别为300、100和20情况下的状态估计结果的状态估计(ee)曲线以及rmse曲线。从中可以看出，随着采样粒子数的减少，传统的重采样方法状态估计精度较低，这是因为传统重采样方法只考虑了对粒子的权值进行重新优化分配，没有考虑观测信息的修正作用，同时舍弃了大量的小权值粒子导致。从图4至图9可以明显看出，传统的重采样方法中rr方法的精度最低，sr的采样精度最好，出现这种情况的原因是由于rr仅考虑了有区别粒子的采样，随着粒子遍历性能的降低，父代粒子越来越少，导致估计精度越来越差。在粒子数n＝300的情况下，gar和本发明方法都保持了较好的跟踪精度，主要原因是gar对小权值粒子进行了优化遗传重组，有效增加了粒子的多样性；但是随着采样粒子数的减少，gar方法的精度明显降低，这是因为gar方法主要缺少最新观测信息的修正作用，产生了估计的累积误差效应。在粒子数n＝20的情况下，本发明方法的估计精度明显优于其他采样方法，并且其跟踪精度和粒子数n＝100基本相近。本发明方法的这一特性，可以明显改善粒子滤波的计算复杂度和滤波精度间的矛盾，提升高精度应用的实时性。实验三：固定采样粒子个数n＝100，分析不同仿真周期情况下的状态估计效果。为分析不同重采样方法的稳定性，针对给定采样粒子数、不同仿真周期情况下的状态估计精度进行了对比分析实验。实验中，设置采样粒子数n＝100，仿真周期t分别为50s、100s和500s，并给出了100步迭代的rmse均值mermse和方差varrmse，具体结果如表2和表3所示。表2n＝100时的mermse重采样方法t＝50st＝100st＝500srr4.88425.76996.6731mr4.68155.26115.4808sr4.37044.78934.0086gar3.06773.53483.9722cer3.04013.48003.4409表3n＝100时的varrmse从表2和表3中可以看出，随着仿真周期的增加，五种方法的估计精度都在降低，其中，rr方法的降低幅度较大，本发明方法和gar方法表现的相对平稳，在50s到500s的变化过程中，gar方法的mermse增加了0.5952，而本发明方法增加了0.4008。由表3可以看出，gar的varrmse波动范围为0.0977，而本发明方法的波动范围为0.065。由此可以看出，本发明方法不仅具有更好的估计精度，同时具有较好的稳定性。实验四：计算复杂度分析。计算的复杂度主要取决于采样、重采样和权值计算三个阶段。在差分进化种群(de)中迭代数目k和粒子数目m决定着本发明方法的复杂度，由于本发明重采样方法进行了粒子分类处理，需要处理的粒子不会超过全部粒子的2/3。具体地，基于matlab2015软件，在intelcore5duo@4ghzpc上进行了计算复杂度对比分析实验。在该对比分析实验中仍然采用如式(19)给定的模型，在相同参数设置情况下进行150次迭代估计。对不同采样粒子数情况下的计算复杂度进行了对比分析，结果如表4所示。表4计算时间消耗比较从表4中可以看出，随着粒子数的增多，五种重采样方法的时间都呈递增趋势，其中，mr、gar和cer的采样时间接近，说明本发明建议分布函数的优化和粒子分类进化重采样步骤在实现避免粒子滤波权值退化、样本贫化的出现，提高粒子滤波方法的精确度和稳定性时，在相同采样粒子个数情况下，计算复杂性并没有增加。本发明针对粒子滤波方法存在的“粒子退化”及“样本贫化”问题，提出了基于分类进化重采样的粒子滤波方法，对传统建议密度函数的优化思想进行了阈值分割，只对退化情况下的建议密度函数进行优化，既有效增强了最新观测信息对粒子的修正作用，也较好地降低了粒子状态优化的时间消耗。同时，在重采样阶段依据粒子的权值大小将采样粒子集分别聚类为小权值粒子种群、保留粒子种群和大权值粒子种群。对大权值粒子后的粒子进行裂变并与小权值粒子进行差分进化，通过突变、交叉、选择产生新粒子种群，避免了对历史小权值粒子信息的丢弃。对比分析实验结果表明，本发明方法有效提升了采样粒子空间分布的多样性，在合理计算复杂性的基础上，有效提升了算法的粒子多样性，在不增加运算负担的前提下提升了滤波的精度，比传统的重采样方法具有更好的精度和稳定性。以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。当前第1页12