CN103684350B - 一种粒子滤波方法 - Google Patents

一种粒子滤波方法 Download PDF

Info

Publication number
CN103684350B
CN103684350B CN201310645786.8A CN201310645786A CN103684350B CN 103684350 B CN103684350 B CN 103684350B CN 201310645786 A CN201310645786 A CN 201310645786A CN 103684350 B CN103684350 B CN 103684350B
Authority
CN
China
Prior art keywords
particle
time
particles
importance
sampling
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201310645786.8A
Other languages
English (en)
Other versions
CN103684350A (zh
Inventor
夏元清
蒲钒
耿秀美
邓志红
付梦印
闫莉萍
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN201310645786.8A priority Critical patent/CN103684350B/zh
Publication of CN103684350A publication Critical patent/CN103684350A/zh
Application granted granted Critical
Publication of CN103684350B publication Critical patent/CN103684350B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Complex Calculations (AREA)
  • Image Processing (AREA)

Abstract

本发明提供一种粒子滤波方法,其包括:步骤1,初始化粒子;步骤2,在k时刻获取测量值,然后利用粒子滤波方法由N个粒子滤波过程并行计算均值和方差,然后进行近似处理获得重要性密度函数并抽取采样粒子;步骤3,根据步骤2获得的重要性密度函数,计算每一个采样粒子的重要性权值;步骤4,将步骤3中得到的重要性权值进行归一化处理;步骤5,根据步骤4中归一化处理后得到的权值进行重采样,得到新的粒子序列;步骤6,对步骤5得到的新的粒子序列xik计算后验概率密度,输出滤波结果。本发明的计算过程简单,能在一定程度上改善粒子退化问题,提高了粒子滤波性能。

Description

一种粒子滤波方法
技术领域
本发明属于非线性滤波技术领域,尤其涉及一种粒子滤波方法。
背景技术
滤波是伴随着信号传递产生的一门技术,信号的传递过程不可避免的要受到内、外部干扰的影响,为了获取所需信号,排除干扰,就要对信号进行滤波。对于非线性系统,通过贝叶斯估计得到精确的最优滤波解是很困难的。常用的非线性滤波方法主要有扩展Kalman滤波(EKF)和无迹Kalman滤波(UKF)等,然而这两种非线性滤波方法会受到非线性程度以及噪声类型的限制。于是又有学者提出了应用范围较广的粒子滤波(PF)算法。粒子滤波是一种基于贝叶斯估计原理的序贯MonteCarlo模拟方法,核心是利用一些随机粒子来表示系统随机变量的后验概率密度,以得到物理模型的近似最优数值解,是一种顺序重要性采样法。简单来说,粒子滤波法是指通过寻找一组在状态空间传播的随机样本对概率密度函数进行近似,以样本均值代替积分运算,从而获得状态最小方差分布的过程。但是,粒子滤波随着粒子的不断迭代,会出现粒子退化现象,使滤波性能变差。
发明内容
为解决上述问题,本发明提供一种粒子滤波方法,该方法由粒子滤波产生重要性密度函数,在一定程度上改善了粒子退化问题,提高了粒子滤波性能。
本发明的粒子滤波方法包括以下步骤:
步骤1,初始化粒子且粒子权值为其中,x0为初始时刻t0的粒子集合,为初始时刻t0第i个状态向量,N为产生的粒子数量,p(x0)为初始概率密度函数;
步骤2,在k时刻获取测量值yk,根据k-1时刻粒子的集合xk-1和测量值yk递推出xk
xk-1和测量值yk利用粒子滤波方法由N个粒子滤波过程并行计算xk的均值和方差然后由式(1)进行近似处理获得重要性密度函数:
q ( x k i | x k - 1 i , y k ) = N ( x ‾ k i , P ‾ k i ) - - - ( 1 ) ;
并根据式(2)抽取采样粒子:
x ^ k i = x ‾ k i + P ‾ k 1 / 2 γ i γ i ~ N ( 0 , I ) - - - ( 2 ) ;
步骤3,根据步骤2获得的重要性密度函数,利用式(3)计算每一个采样粒子的重要性权值:
ω ~ k i = ω k - 1 i p ( y k | x k i ) p ( x k i | x k - 1 i ) q ( x k i | x k - 1 i , y k ) - - - ( 3 ) ;
步骤4,将步骤3中得到的重要性权值利用式(4)进行归一化处理:
ω k i = ω ~ k i sum { ω ~ k i } i = 1 N - - - ( 4 ) ;
步骤5,根据步骤4中归一化处理后得到的权值进行重采样,得到新的粒子序列
步骤6,按公式对步骤5得到的新的粒子序列计算后验概率密度,输出滤波结果;
进一步的, x k = 0.5 x k - 1 + 25 x k - 1 1 + x k - 1 2 + 8 cos ( 1.2 ( k - 1 ) ) + w k 是所用非线性模型的状态方程,是观测方程,其中,wk,vk分别为过程噪声和观测噪声,均为零均值高斯白噪声,且wk,vk的方差分别假设为1,初始状态为x0=0.1。
进一步的,所述步骤5中还可以采用精细重采样方法得到新的粒子序列。
本发明的有益效果在于:
本发明提出由粒子滤波产生重要性密度函数的重要性采样方法,与由扩展Kalman滤波、无迹Kalman滤波产生重要性密度函数的重要性采样方法相比较,本发明计算过程简单,精度较高,能在一定程度上改善粒子退化问题,提高了粒子滤波性能。
附图说明
图1为本发明的EKF-PF,UKF-PF以及PF-PF的状态估计示意图;
图2为本发明的EKF-PF,UKF-PF以及PF-PF的RMSE示意图;
图3为本发明的采用精细重采样的EKF-PF,UKF-PF以及PF-PF的状态估计示意图;
图4为本发明的采用精细重采样的EKF-PF,UKF-PF以及PF-PF的RMSE示意图;
图5为本发明的不采用精细重采样的PF-PF,EKF,UKF以及PF的状态估计示意图;
图6为本发明的不采用精细重采样的PF-PF,EKF,UKF以及PF的RMSE示意图。
具体实施方式
重要性采样方法(即获取重要性密度函数的方法)对粒子滤波的效果有较大的影响,若重要性采样方法不合适,会导致严重的粒子退化,影响滤波效果。现有的粒子滤波算法中,产生重要性密度函数的方法主要有EKF和UKF,但是这两种方法计算复杂,精度受限。本发明提供粒子滤波方法,其利用改进的粒子滤波算法,由粒子滤波产生重要性密度函数,得到一种新的滤波方法,称为PF-PF。
该方法步骤如下:
步骤1,初始化粒子,粒子权值其中,x0为初始时刻(t0时刻)粒子的集合。为t0时刻第i个状态向量(称为粒子),N为产生的粒子的数量,p(x0)为初始概率密度函数。
步骤2,重要性采样,在k时刻,获取测量值yk。根据xk-1和测量值yk,利用粒子滤波方法,由N个粒子滤波过程并行计算xk的均值和方差然后重要性密度函数由式(1)近似:
q ( x k i | x k - 1 i , y k ) = N ( x ‾ k i , P ‾ k i ) - - - ( 1 )
并根据式(2)抽取采样粒子:
x ^ k i = x ‾ k i + P ‾ k 1 / 2 γ i γ i ~ N ( 0 , I ) - - - ( 2 ) .
其中PF,EKF,UKF重要性密度函数都是(1)式,只不过估计的方式不同,分别利用EKF,UKF和PF来估计
步骤3,根据步骤2得出的重要性密度函数,计算每一个采样粒子的重要性权值:
ω ~ k i = ω k - 1 i p ( y k | x k i ) p ( x k i | x k - 1 i ) q ( x k i | x k - 1 i , y k )
步骤4,将步骤3中得到的权值归一化:
ω k i = ω ~ k i sum { ω ~ k i } i = 1 N
步骤5,根据得到的权值,进行重采样,得到新的粒子序列
步骤6,按式计算后验概率密度,输出滤波结果,后验概率密度就是最终滤波结果。
其中,xk-1的获得:由于x0已知,则可以根据步骤2至步骤6,由x0以及观测值y1可递推出x1,由x1及y2递推出x2
由于利用粒滤波算法对于粒子均值和方差的估计精度较EKF和UKF更高,因此PF生成的重要性分布与真实的后验概率分布更接近,从而提高采样质量。因此PF-PF算法能够有效防止粒子退化。
此外,步骤5还可以采用精细重采样方法,进一步提高PF-PF的滤波性能。精细重采样方法与传统重采样方法相比,优点是:1、在保留或淘汰粒子上尽可能的保留了权重更大的粒子;2、根据母粒子产生新粒子,不单纯完全复制粒子。这两方面的处理可以有效地改善粒子退化问题。
将由PF产生重要性密度函数,与由EKF、UKF产生重要性密度函数的重要性采样方法进行比较,对采用三种重要性采样方法的粒子滤波算法进行了一系列仿真实验。由EKF产生重要性密度函数的粒子滤波方法(EKF-PF)、由UKF产生重要性密度函数的粒子滤波方法(UKF-PF)、由PF产生重要性密度函数的粒子滤波方法(PF-PF)的状态估计情况如图1、图2所示,验证了采用新的重要性采样方法的粒子滤波的状态估计值更接近真实状态,RMSE更小。
为了进一步提高滤波的精度,设计了由PF、EKF、UKF产生重要性密度函数,同时,重采样阶段分别采用精细重采样方法的一组仿真实验,EKF-PF,UKF-PF,PF-PF的状态估计情况如图3、图4所示,可以看出,PF-PF在三种滤波方法中性能最好,同时从RMSE的大小情况来看,采用更好的重采样方法能够在一定程度上改进滤波性能。
此外,将PF-PF方法与传统的EKF,UKF以及PF进行对比,结果如图5、图6所示,可以明显看出PF-PF滤波性能的优越性,即图中可以看出,PF-PF有更小的RMSE,也就是说明有更好的滤波性能。
仿真的参数设置为:
所用非线性模型的状态方程为 x k = 0.5 x k - 1 + 25 x k - 1 1 + x k - 1 2 + 8 cos ( 1.2 ( k - 1 ) ) + w k .
观测方程为为 y k = x k 2 20 + v k .
其中,wk,vk分别为过程噪声和观测噪声,均为零均值高斯白噪声,且wk,vk的方差分别假设为1,初始状态为x0=0.1。
当然,本发明还可有其他多种实施例,在不背离本发明精神及其实质的情况下,熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形,但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (1)

1.一种粒子滤波方法,其特征在于,包括以下步骤:
步骤1,初始化粒子且粒子权值为其中,x0为初始时刻t0的粒子集合,为初始时刻t0第i个状态向量,N为产生的粒子数量,p(x0)为初始概率密度函数;
步骤2,在k时刻获取测量值yk,根据k-1时刻粒子的集合xk-1和测量值yk
由N个粒子滤波过程并行计算xk的均值和方差然后由式(1)进行近似处理获得重要性密度函数:
q ( x k i | x k - 1 i , y k ) = N ( x ‾ k i , P ‾ k i ) - - - ( 1 ) ;
其中,即是预测出的k时刻粒子第i个状态向量即第i个粒子的均值,是k时刻粒子第i个状态向量的方差,表示均值为方差为
并根据式(2)抽取采样粒子:
x ^ k i = x ‾ k i + P ‾ k 1 / 2 γ i γ i ~ N ( 0 , 1 ) - - - ( 2 ) ;
其中,各粒子均值为0,各个粒子方差均为1,即标准正态分布;
步骤3,根据步骤2获得的重要性密度函数,利用式(3)计算每一个采样粒子的重要性权值:
ω ~ k i = ω k - 1 i p ( y k | x k i ) p ( x k i | x k - 1 i ) q ( x k i | x k - 1 i , y k ) - - - ( 3 ) ;
步骤4,将步骤3中得到的重要性权值利用式(4)进行归一化处理:
ω k i = ω ~ k i s u m { ω ~ k i } i = 1 N - - - ( 4 ) ;
步骤5,根据步骤4中归一化处理后得到的权值进行重采样,得到新的粒子序列
步骤6,按公式对步骤5得到的新的粒子序列计算后验概率密度,输出滤波结果;
其中,x0:k是从采样时刻0到采样时刻k,各个时刻粒子的集合,用来表示第0次采样到第k次采样的状态值;y1:k表示1时刻到k时刻各个时刻状态的测量值;而则表示0时刻到k时刻,每一时刻第i个粒子;而δ表示采样时间间隔,即xk与xk-1之间的时间差;
进一步的,是所用非线性模型的状态方程,是观测方程,其中,wk,vk分别为过程噪声和观测噪声,均为零均值高斯白噪声,且wk,vk的方差分别假设为1,初始状态为x0=0.1。
CN201310645786.8A 2013-12-04 2013-12-04 一种粒子滤波方法 Active CN103684350B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310645786.8A CN103684350B (zh) 2013-12-04 2013-12-04 一种粒子滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310645786.8A CN103684350B (zh) 2013-12-04 2013-12-04 一种粒子滤波方法

Publications (2)

Publication Number Publication Date
CN103684350A CN103684350A (zh) 2014-03-26
CN103684350B true CN103684350B (zh) 2016-07-13

Family

ID=50320843

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310645786.8A Active CN103684350B (zh) 2013-12-04 2013-12-04 一种粒子滤波方法

Country Status (1)

Country Link
CN (1) CN103684350B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105157704A (zh) * 2015-06-03 2015-12-16 北京理工大学 一种基于贝叶斯估计的粒子滤波重力辅助惯导匹配方法
CN105180938B (zh) * 2015-06-03 2018-02-02 北京理工大学 一种基于粒子滤波的重力采样矢量匹配定位方法
CN107919940A (zh) * 2016-10-10 2018-04-17 深圳超级数据链技术有限公司 适用于OvXDM系统的前向后向平滑译码方法、装置及OvXDM系统
CN107919939B (zh) * 2016-10-10 2021-09-03 深圳光启合众科技有限公司 适用于OvXDM系统的双滤波平滑译码方法、装置及OvXDM系统
CN107124159B (zh) * 2017-04-27 2020-06-05 鲁东大学 一种基于自适应kld盒子长度的粒子滤波器的实现方法
CN109492769A (zh) * 2018-10-31 2019-03-19 深圳大学 一种粒子滤波方法、系统和计算机可读存储介质
WO2020087362A1 (zh) * 2018-10-31 2020-05-07 深圳大学 一种粒子滤波方法、系统和计算机可读存储介质
CN112613222B (zh) * 2021-01-04 2023-09-15 重庆邮电大学 基于改进粒子滤波的倾斜探测电离层muf短期预报方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6055119A (en) * 1997-02-21 2000-04-25 Samsung Electronics Co., Ltd. Adaptive signal processing method and circuit for a digital recording/reproducing apparatus
CN101826856A (zh) * 2010-03-11 2010-09-08 哈尔滨工程大学 基于超球面采样无迹卡尔曼滤波的粒子滤波方法
CN103389094A (zh) * 2013-07-15 2013-11-13 哈尔滨工程大学 一种改进的粒子滤波方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6055119A (en) * 1997-02-21 2000-04-25 Samsung Electronics Co., Ltd. Adaptive signal processing method and circuit for a digital recording/reproducing apparatus
CN101826856A (zh) * 2010-03-11 2010-09-08 哈尔滨工程大学 基于超球面采样无迹卡尔曼滤波的粒子滤波方法
CN103389094A (zh) * 2013-07-15 2013-11-13 哈尔滨工程大学 一种改进的粒子滤波方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
一种改进的PSO算法;刘羿彤 等;《第二十六届中国控制会议论文集》;20070726;第479-483页 *

Also Published As

Publication number Publication date
CN103684350A (zh) 2014-03-26

Similar Documents

Publication Publication Date Title
CN103684350B (zh) 一种粒子滤波方法
Tung et al. Detecting chaos in heavy-noise environments
CN102721545B (zh) 一种基于多特征参量的滚动轴承故障诊断方法
Grant et al. Comparison of matrix pencil and prony methods for power system modal analysis of noisy signals
CN103679643B (zh) 一种多条纹噪声定位滤除方法
Zhang et al. An improved filtering method based on EEMD and wavelet-threshold for modal parameter identification of hydraulic structure
US20220020123A1 (en) Bayesian image denoising method based on distribution constraint of noisy images
CN102222911A (zh) 基于ar模型和卡尔曼滤波的电力系统间谐波估计方法
CN110175541B (zh) 一种海平面变化非线性趋势提取的方法
CN101561314A (zh) 随机共振-混沌微弱信号检测方法
CN104677632A (zh) 利用粒子滤波与谱峭度的滚动轴承故障诊断方法
CN105785324A (zh) 基于mgcstft的线性调频信号参数估计方法
CN103944174A (zh) 基于互相关函数滤噪算法的低频振荡在线辨识方法
CN107729845B (zh) 一种基于子空间特征值分解的实测频响函数降噪方法
CN107037486A (zh) 地球天然脉冲电磁场数据处理的时频谱分析方法及系统
Chen et al. Automatic modulation classification of radar signals utilizing X-net
CN117688422A (zh) 基于改进稀疏分量分析的欠定模态参数识别方法、计算机设备及存储介质
CN117309079B (zh) 基于时差法的超声飞渡时间测量方法、装置、设备及介质
Zhou et al. Application of IPSO-MCKD-IVMD-CAF in the compound fault diagnosis of rolling bearing
Liu et al. Bearing failure diagnosis at time-varying speed based on adaptive clustered fractional Gabor transform
CN104111109B (zh) 一种基于不同阶次统计量及支持向量机的机械振动状态识别方法
CN106338651A (zh) 应用于电力系统低频振荡模式识别的粒子滤波分析方法
Xiaomeng et al. A sensor fault diagnosis method research based on wavelet transform and hilbert-huang transform
CN112014811B (zh) 一种雷达载波频率的精细估计方法
CN103166887B (zh) 一种边带调制信号分类的方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant