（交通信息工程及控制专业论文）基于FPGA实现的粒子滤波算法研究.pdf

中文摘要 摘要：粒子滤波是上世纪9 0 年代发展起来的一种基于蒙特卡罗方法和递推贝叶斯 估计的新滤波方法，在处理非线性、非高斯系统的参数估计和状态滤波方而具有 独到的优势。但是其复杂的算法结构、庞大的计算量和缓慢的运算速度限制了其 在实时系统中的应用。本文研究的目的是降低粒子滤波算法的复杂度，提高运算 速度，设计一种运算速度快、性能可靠、占用硬件资源少的粒子滤波器，使其能 应用于目标跟踪等实时系统中。 现场町编程门阵y u ( f i e l dp r o g r a m m a b l eg a t ea r r a y ，f p g a ) 是一种硬件逻辑器 件，执行真正意义卜的并行运算。理论卜，个时钟周期，f p g a 能输一个结果， 粒了滤波在f p g a 系统中能达到最快运算速度。所以本文选取f p g a 作为粒了滤 波硬件实现的f | 标器件。 本文首先介绍了粒| 了滤波的现状和本文研究的意义。接着详细介绍了粒了滤 波基本原理和f p g a 的结构特点。然后对s i r f 算法的步骤、复杂度、并行件进行深 入分析，在此基础上对重采样、权值计算、采样三步进行改进，降低了粒子滤波 算法的复杂度，提高了运算速度。并以被动定位系统中目标跟踪为例进行m a t l a b 仿真，验证改进后算法的正确性。然后基。j ：改进后的s i r f 算法，在f p g a 中设汁出高 斯随机数生成器、重采样、采样、权值计算各个模块，并进行波形图仿真，验证 模块设计的正确性。最后，连接各个模块，在f p g a 中实现粒子滤波算法。 本文在f p g a 中设计并实现了这种改进的粒子滤波算法。波形图仿真和f p g a 资 源使用情况报告表明这种算法器具有速度快、占用资源少的特点，能j 泛适用于 目标跟踪、导航与制导、图像处理、参数估计等实时系统中。 关键词：粒子滤波；f p g a ；目标跟踪；并行算法 分类号：u 2 8 3 1 a b s t r a c t a b s t r a c t ：d e v e l o p i n g i n9 0 sl a s tc e n t u r y ，p a r t i c l ef ii t e r , an e wf i l t e rm e t h o db a s e d o nm o n t ec a r l oa n dr e c u r s i v eb a y e s i a ne s t i m a t i o n ，h a ss p e c i a la d v a n t a g e si nd e a l i n gw i t h t h es t a t ea n dt h ep a r a m e t e re s t i m a t i o ni nt h en o n l i n e a ra n dn o n - g a u s s i a ns y s t e m h o w e v e r ， t h ed i s a d v a n t a g e so fc o m p l e xa l g o r i t h ma r c h i t e c t u r e ，e n o r m o u sc o m p u t a t i o n sa n dl o w s p e e dh a v er e s t r i c t e di t si m p l e m e n t a t i o ni nr e a l - t i m es y s t e m t h i st h e s i sa i m st or e d u c e t h ec o m p l e x i t yo ft h ep a r t i c l ef i l t e ra l g o r i t h mw i t hb e t t e rc o m p u t a t i o ns p e e di no r d e rt o d e s i g nap a r t i c l ef i l t e rw h i c hi sf a s t ，c r e d i b l ea n du s i n gl e s sh a r d w a r er e s o u r c es oa st o b ei m p l e m e n t e di nr e a l t i m es y s t e m ，i e o b j e c tt r a c k i n ge t c f i e l dp r o g r a m m a b l eg a t ea r r a yi sah a r d w a r el o g i cd e v i c e ，w h i c hi m p l e m e n t st h e r e a lp a r a l l e lc o m p u t a t i o n i nt h e o r y , t h ef p g ac a no u t p u to n er e s u l ti no n ec l o c kp e r i o d p a r t i c l ef i l t e rc a nr e a c ht h eh i g h e s tc o m p u t a t i o ns p e e di nf p g as y s t e m s o ，t h i st h e s i s c h o o s e st h ef p g aa st h et a r g e td e v i c et oi m p l e m e n tp a r t i c l ef i l t e r t h i st h e s i s ，a tf i r s t ，i n t r o d u c e st h es t a t u si nq u oo fp a r t i c l ef i l t e ra n dt h em e a n i n go f r e s e a r c h s e c o n d l y ，i t i n t r o d u c e st h eb a s i c st h e o r yo ft h e p a r t i c l e f il t e ra n dt h e c h a r a c t e r i s t i c so ff p g a t h i r d l y ，t h ec o m p l e x i t yo ft h ea l g o r i t h mh a sb e e nr e d u c e da n d t h ec o m p u t a t i o ns p e e dh a sb e e na d v a n c e db yi m p r o v i n gm e t h o d so nt h et h r e es t e p si n s a m p l i n g ，w e i g h t i n gc o m p u t a t i o n ，a n dr e s a m p l i n gb a s e do na n a l y z i n gt h o r o u g h l yt h e p r o c e s s e s ，c o m p l e x i t y ，p a r a l l e lo fs i r fa l g o r i t h m 。a n di t s i m u l a t e si nm a t l a bu s i n g b e a r i n g - o n l ys y s t e mo b j e c tt r a c k i n g f o re x a m p l et oc h e c kt h ec o r r e c t n e s so ft h e a d v a n c e da l g o r i t h m s a n dt h e n ，i td e s i g n st h em o d e l s ，s u c ha st h eg a u s sr a n dn u m b e r g e n e r a t o r ，t h er e s a m p l i n gm o d e l ，s a m p l em o d e la n dw e i g h tc o m p u t a t i o nm o d e li n f p g ai nt h eb a s i so ft h ea d v a n c e ds i r fa l g o r i t h m ，a sw e l la su s e st h ew a v e f o r mb e n c h s i m u l m i o nt ot e s tt h ec o r r e c t n e s so ft h em o d e l s a tl a s t ，t h e s em o d e l sa r ec o n n e c t e dt o i m p l e m e n tt h ep a r t i c l ef i l t e ra l g o r i t h mi nf p g a t h i st h e s i sh a sd e s i g n e da n di m p l e m e n t e dt h ea d v a n c e dp a r t i c l ef i l t e ra l g o r i t h m m e t h o di nf p g a ，a n dt h ew a v e f o r mb e n c hs i m u l a t i o na n dt h er e p o r to fr e s o u r c eu s i n g i n d i c a t et h a tt h i sa d v a n c e da l g o r i t h m ，w i t hi t sc h a r a c t e r i s t i c so fh i g hs p e e d ，u s i n gs m a l l r e s o u r c e ，c a nb eu s e di nr e a l - t i m es y s t e m ，i e o b j e c tt r a c k i n g ，n a v i g a t i o n ，i m a g e p r o c e s s i n g ，p a r a m e t e re s t i m a t i o ne t c k e y w o r d s ：p a r t i c l ef i l t e r ；f p g a ；o b j e c tt r a c k i n g ；p a r a l l e la l g o r i t h m c l a s s n o ：1 3 2 8 3 1 学位论文版权使用授权书 本学位论文作者完全了解北京交通大学有关保留、使用学位论文的规定。特 授权北京交通大学可以将学位论文的全部或部分内容编入有关数据库进行检索， 并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校向困 家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名： 签字日期：d 8 年 导师签名：7 彩锣日 签字日期： 6 苫年石月日 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的研 究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表或 撰写过的研究成果，也不包含为获得北京交通大学或其他教育机构的学位或证书 而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作 了明确的说明并表示了谢意。 学位论文作者签名：签字日期：年月口 5 8 致谢 深深感谢我的导师张三同副教授，正是张老师对我的指导、关怀和严格要求， 使我顺利完成了硕士阶段的学习和研究。在导师渊博的知识、严谨的工作作风、 敏锐的洞察力和高尚的师德的熏陶下，我的学识快速增长，提高了科研技能，具 备了科技人员的优良作风、品格和素质，树立了正确的人生观和批界观。正足张 老师的严格要求，通过两年的潜心学习和研究，我才对粒子滤波算法及f p g a 实 现有了比较深入的研究，在粒子滤波算法f p g a 实现方面提出了自己的一些观点、 见解和实施方案，完成了论文相关的科研t 作和撰写t 作。 张三同老师悉心指导我完成了实验室的科研工作，存学习和生活卜都给予了 我很人的关心和帮助，在此向张老师表示衷心的感谢。 在实验室工作和撰写论文期间，实验室各位同学对我的研究工作给予了热情 帮助，在此向大家表达我的感激之情。 最后我也要特别感谢我的家人和朋友们，是他们给予我精神上、生活上的支 持和关怀，使我安心从事学习和研究工作，彳得以完成学业。 1 引言 滤波，简单米说就是从一些混杂有噪声的信号中滤除噪声，得到有用信号。 有的信号的变化规律是确定的，而有的则是随机的。对于前者，可以根据信号频 率特件设计相应的滤波器，使有用信号无衰减通过，而噪声被抑制，限制其通过， 如低通滤波器、高通滤波器；对于后者，无法用常规滤波器来抑制噪声，但随机 信号有确定的功率谱，可以根据有用信号和噪声的功率谱设计滤波器，如卡尔曼 滤波器( k a l m a nf i l t e r ，k f ) 。r 尔曼滤波器足一种线性、最小方差的最优状态估计， 广泛应用于信号、通信、控制等领域。但卡尔曼滤波算法不能解决非线性、非高 斯问题。为启犁决非线性f i u j 题，最为常用的滤波方法足扩展卡尔曼滤波( e x t e n d e d k a l m a nf i l t e r , e k f ) 算法，但该方法的原理是把非线性问题线性化，将导致模型误 差，只适用于弱非线性系统，对于强非线性系统很容易导致发散。近来，j u l i e r 等 提出一种新的推广卡尔曼滤波算法一非追踪卡尔曼滤波( u n s c e n t e d k a n l m a n f i l t e r , u k f ) ，该算法直接利用非线性模犁，避免引入线性误差，t i 二e k f 具有更 高的滤波精度。但该方法仍需假设后验概率密度函数为高斯分布，只能用于高斯 分布模犁，使用受到限制。 粒了滤波( p a r t i c l ef i l t e r , p f ) 是上世纪9 0 年代发展起来的一种基于蒙特卡罗方 法( m o n t ec a r l om e t h o d s ，m c m ) 和递推贝叶斯估计的新的用于非线性、非高斯 系统的滤波方法。粒了滤波的基本思想是：首先依据系统状态向量的经验条件分 布，在状态空间抽样产生。组随机样本集合，这些样本集合称为“粒子”；然后根 据观测值彳i 断地调整粒子的权重大小和样本位置；最后，通过调整后的粒子的信 息修正最初的经验条件分布，估计出系统状态和参数。该算法是种递推滤波算 法，原则上可以用来估计任意非线性、非高斯随机系统的状态和参数，有效克服了 k f 、e k f 等滤波算法的缺点。 1 1 粒子滤波的研究现状和发展方向 在2 0 世纪5 0 年代，出现了一种被称为“序贯重要性采样( s e q u e n t i a li m p o r t a n c e s a m p l e ，s i s ) ”的蒙特卡罗方法【2 1 ，它通过离散的随机测度逼近概率分布，应用在 物理和工程领域。然而研究人员始终未能解决退化( d e g e n e r a c y ) i 矗j 题和计算量的制 约，相当长一段时间s i s 未能进一步发展。直到1 9 9 3 年，g o r d o n 等【3 1 提出了一种新 的基t s i s i 排线性滤波方法，采用重采样( r e s a m p l i n g ) 克服早期算法的粒子退化问 题，出现了第一个可操作的m o n t ec a r l o 滤波器，从而奠定了粒子滤波算法的基础。 现代计算机技术使m o n t ec a r l o 滤波方法得到迅速发展，这些m o n t ec a r l o 滤波方法 在各种领域分别被称为b 0 0 t s t r a p 【3 1 ，适者生存( s u r v i v a lo f t h ef i t t e s t ) t 4 l ，s e q u e n t i a l m o n t ec a r l o t 5 j f 6 1 等，现在通称为粒子滤波。现代计算机技术的发展以及粒子滤波方 法本身具有的巨大潜力使得其成为当前一个非常活跃的研究领域。粒子滤波被广 泛应用于通信、目标跟踪、导航与制导、图像处理、故障诊断、参数估计、模式 识别等方面。 虽然粒子滤波算法能适用于任意非线性、非高斯系统，但其本身也存在一些 问题，如粒子退化。g o r d o n 等人采用重采样解决这个问题的同时又带米粒子耗尽 问题。许多研究人员对粒子滤波加以改进，提出各种改进算法，如改进重要性密 度函数选取方法的u 一粒子滤波( u n s c e n t e dp a r t i c l ef i l t e r ，u p f ) ；对再采样过程 进行优化的i f 则粒子滤波( r e g u l a r i z e dp a r t i c l ef i l t e r ，r p f ) ；增加粒子多样性的 辅助粒子滤波( a u x i l i a r yp a r t i c l ef i l t e r ，a p f ) ；高斯滤波和粒子滤波结合的高斯 粒子滤波( g a u s s i a np a r t i c l ef i l t e r ，g p f ) 等。 u k f 用高斯分布来逼近状态分布，但不进行线性化，认为状态密度分布可以 由有限个称为s i g m a 点的样本来描述。这些样本点通过非线性模型后可以得到精确 至l j t a y l o r 级数三阶展开式相当的均值和方差。m e r w e 等人【7 】提出用u k f 来产生粒子 滤波的重要性密度函数万【ix 0 ：足一1 ，z l ：kj ，即： ( l 础一i ，刁： ) = ( n ，曩o ) i = l ，2 ，n ( 1 1 ) 其中或n ，砭是通过u k f 算法得到的均值和方差。 u k f 产生的重要性密度函数更加逼近于概率密度函数( p r o b a b i l i t yd e n s i t yf u n c t i o n ， p d f ) ，因而估计精度更高。这种改进的粒子滤波算法被称为u p f 。 通常的重采样是在离散分布中采样，m u s s o 等提r p f 8 1 ，从连续分布式： 多( x k lz 0 ：k ) 以n k h ( x , 一o ) f = 1 ，2 ，( 1 2 ) 其中蚝o ) = ( 去) k ( 云) 是尺度可调的k e m e l 密度，j i i o 为核带宽，为x 的维 数。 中重采样，可以有效缓解重采样过程中粒子多样性丧失问题，在过程噪声较小的 情况下可以获的较好的滤波精度。 p i t t 等人【9 1 提出了一种改进的粒子滤波算法一a p f ，引入辅助变量，选择重要性 2 密度函数满足式： 刀r ( x k ，m i z o ：) 芘p ( z ki 从d ) p ( x ki 三) 嘭一l f = l ，2 ，n ( 1 3 ) 其中m 为辅助变量，表示k 一1 时刻采样粒子的索引，? 为状态转移先验概率 密度的均值。 a p f 辅助采样在k l 时刻进行，并考虑k 时刻的测量信息乙，这样使得在过程噪声 较小时，采样粒子更接近k 时刻的真实状态，从而可获得比标准粒子滤波更高的滤 波精度：但当过程噪声较大时，滤波效果不如标准粒子滤波算法。 k o t e c h a ，d i u r i c 1 0 l 提出了另一种改进粒子滤波算法，将高斯滤波和粒子滤波结 合，称为g p f 。这种算法采用粒子滤波得到高斯密度函数，通过这一高斯密度函数 来近似估计状态的后验概率分布，即： p ( x kiz o ) g p ( z ii 耳) ( _ ；历，艺女) ( 1 4 ) p ( x k ) = n ( x k ；1 1 k ， ) ( 1 5 ) 其中q = ( ，p ( iz o ：) p ( 气ix k ) d x i ) ，一t k 、艺。为随机变量xk 的均值 和方差。式中的均值t 和方差t 将通过粒子滤波得到： 厄= 墨。以o o ( 1 6 ) 艺i = 兰。w ? ( n 一以) ( n 一心) 7 ( 1 7 ) 与其它的粒子滤波相比，g p f 用高斯分布来近似后验概率分布，不会产生粒子退化 问题，从向不需要对粒子进行重采样，使得计算量、复杂度降低。但后验概率分 布不能用高斯分布近似的非线性动态空问模型或非线性、非加性高斯噪声模型， 使用高斯粒子滤波效果不佳。 粒子滤波适用于各种非线性、非高斯系统的参数估计和状态滤波，对其的研 究已取得可喜的进步，应用领域不断扩展。但由于粒子滤波是近年来出现的新算 法，算法本身还不很成熟，仍有大量的问题亟待解决，主要体现在以下几个方面： 1 ) 重要性函数的选取。重要性函数直接影响剑粒子滤波性能的高低。虽然出 现了许多针对实际问题的改进算法，但对通用问题缺乏一般借鉴意义。如何选取 重要性函数成为目前粒子滤波研究的一个热点。 2 ) 粒子滤波的收敛性。从数学基础上看，粒了滤波的收敛性尚未得到解决。 但实际| 、口j 题中要求计算结果收敛到真实值，且必须具有一定的收敛速度。如何找 到一个可靠的途径检测粒子滤波收敛性和收敛速度是至关重要的。 3 3 ) 拓展粒子滤波的应用领域。由于粒子滤波是近年才出现的滤波算法，它的 应用还限于几个领域。原则上讲在非线性、非高斯系统中原有滤波器滤波效果不 佳的情况下都可以尝试采用粒子滤波。不断拓展粒子滤波的应用领域也将成粒子 滤波研究的重要方向。 4 ) 粒子滤波的硬件实现。粒子滤波虽然具有适h j 各种非线性、非高斯系统的 优点，但庞大的计算景制约了其实时系统中的应用。如何简化计算复杂度、提高 运算速度，设计可靠的算法器将是粒子滤波实际应用的关键。 1 2 本文研究的背景和意义 粒子滤波的硬件实现研究是目前国内外正在蓬勃发展的粒子滤波算法研究领 域的一个重要组成部分。在无线通信、目标跟踪、导航与制导等非线性、非高斯 场合中，传统的卡尔曼滤波k f 、扩展卡尔曼滤波e k f 很难精确的实现参数估计和 状态滤波，出现精度不高、可靠性差，甚至发散等问题。而粒子滤波被认为能很 好的适用于这些场合。许多文献对粒子滤波在这些领域的应用进行了理论上的研 究【 f 1 2 l 。但是如何利用硬件设计粒子滤波器，应用剑实际的实时系统巾就很少有 文献提剑。在实际的应用系统巾，通常需要实时的计算出动态参数，估计出现有 状态，如目标跟踪系统，需要根据含有测量噪声的目标方位角估计出目标在j 维 坐标系中的方位和速度；而住无线通信中，对于通过衰落信道传输的数据，我们 需要实时进行序列的信道估计和数据检测。而粒子滤波算法中，为了得到很好的 滤波性能，估计每一个状态通常需要成百上千个粒子，多维的状态向量更是成倍 地增加了计算量。粒子滤波主要有三步基本操作：采样( 从小含观察值的状态空 间产生新的粒子) 、权值计算( 根据观测值计算各个粒子的权值) 、重采样( 抛弃 权值小的粒子，使用权值大的粒子代替) ，这三步构成了粒子滤波的基本算法 s i r f ( s a m p l ei m p o r t a n c er e s a m p l i n gf i l t e o 。粒子滤波三步循环执行，加卜输出步 骤，最后才估计出状态值，这使得算法结构复杂。计算量大、结构复杂成为粒子 滤波应用于实时系统中的最大瓶颈。 为了快速估计出结果，在计算量大的情况下需要很高的计算速度，并行计算 是一种有效的途径。软件意义上的并行在微观上其实是交替执行的，而硬件上的 并行才是真止意义上的并行操作。考虑到粒子滤波巾各个粒子的独立性和操作的 并行性，设计一个高速的并行硬件粒子滤波算法器是可行的。粒子滤波硬件实现 的基本思想是：把粒子滤波计算过程划分为采样、权值计算、重采样三个步骤， 按粒子滤波的各个步骤设计相应的硬件模块，使得各个模块相互独立，流水线并 4 行执行。采样、权值计算两个阶段中，各个粒子是相互独立的，可以并行运算， 而重采样需要考虑所有粒子的似然分布，只有待全部粒子权重计算完成才能开始 执行。如何简化重采样算法，减弱粒子间的相互联系，使得重采样过程能和前面 两步并行执行，也将是粒子滤波硬件实现研究的重要内容。 对于粒子滤波的硬件实现，主要有两个方面的研究内容：一是简化粒子滤波 算法。粒子滤波算法结构复杂，含有大景乘除、指数运算、三角运算，使得计算 量庞大、运算缓慢。如何优化算法，减少乘除、指数运算等非线性运算是研究的 一个重要内容；二是各个模块的硬什设计。如何设计运算速度快、性能可靠、占 用硬件资源少的符个模块也是粒子滤波硬件实现研究的一个重要内容。采用一种 存减少硬件复杂度、提高运算速度的同时，保持滤波性能的有效机制是粒子滤波 硬件实现的关键点1 3 1 。 粒子滤波在无线通信、目标跟踪等非线性、非高斯系统中具有广阔的发展前 景，但计算量大、实时性差成为了其在实际中应用的瓶颈。降低粒子滤波算法的 复杂度，提高运算速度，设计一种运算速度快、性能可靠、占用硬件资源少的粒 子滤波算法器将具有重要的理论和现实意义。 1 3 本文研究的主要内容 随着研究发展，粒予滤波扩展到不同的领域，许多研究人员认为在非线性、 非高斯系统中，粒子滤波是获得最好的滤波性能的滤波算法之一。虽然与e k f 、 u k f 等非线性滤波算法比有很大的优势，但其本身也存在着实时性差等问题。设 计一种运算速度快、性能可靠、占用硬件资源少的粒子滤波算法硬件模块将是本 文研究的主要内容。 随这数字信号处理技术的发展，d s p 芯片( d i g i t a ls i g n a lp r o c e s s o r ) 在各个领域 大显身手，但在运行粒子滤波算法时遇到了很大挑战。通常，为了得到很好的滤 波性能，估计出一个状态值需要成百上千个粒子，当状态维数数增加时，计算量 更是成倍增加，使得运算速度缓慢。在二维纯方位目标跟踪系统中，使用1 0 2 4 个 粒子的粒子滤波d s p 系统的运算速度只达到几h z ，严重影响了其在实时系统中应 用。 与d s p 芯片软件上的并行相比，现场可编程门阵歹u ( f i e l dp r o g r a m m a b l eg a t e a r r a y ，f p g a ) 才是真正意义上的并行运算。f p g a 具有运算速度快、可靠性高、 编程灵活等优点。理论上一个时钟周期，f p g a 输出一个结果。所以理论上粒子滤 5 波在f p g a 系统中最快速度可以达到尼= k ) ，其中n 为粒子数，为系统时 钟周期。所以本文把f p g a 作为粒子滤波硬件实现的目标器件。 本论文的组织结构为： 第一章：引言。简单介绍粒子滤波的发展现状和本课题的研究背景和意义。 第二章：粒子滤波算法原理。首先介绍所用到的状态空问模型，再简要说明 贝叶斯估计理论和蒙特卡罗方法，最后介绍基本的粒了滤波算法原理。 第三章：可编程逻辑器件和f p g a 。首先介绍了可编程器件的发展情况和f p g a 的硬件结构。再主要介绍x i l i n x 公司高阡能v i r t e x - 4 系列f p g a 芯片。 第四章：基于f p g a 实现的粒子滤波算法改进。分析粒子滤波各个步骤的算 法结构，针埘每一个模块进行简化改进。最后进行了m a t l a b 算法仿真。 第五章：基i j ：改进后s i r f 算法的目标跟踪在f p g a 中的实现。以被动定位系 统中二维纯方位目标跟踪为实现对象设计f p g a 模块。首先介绍了被动定位系统 中的h 标跟踪问题，再设计了s i r f 算法实现的总体结构和参数设定，最后，针对 每一滤波步骤设计出f p g a 模块并进行波形图仿真。 第六章：总结和展望。对本文t 作的总结和粒子滤波硬件实现将来发展的展 望。 6 2 粒子滤波算法原理 本章将首先介绍状态空间模型，然后介绍贝叶斯估计和蒙特卡罗方法，在此 基础上引入粒子滤波，对该算法原理进行分析研究。 2 1 状态空间模型 状态空间模型是描述动态系统的模型，隐含时间为白变量，它表达了由于输 入引起系统内部状态的变化，并由此使输出发生的变化。在信号处理、通信、目 标跟踪等动态系统中，许多问题都丌j 以由状态空间模型加以描述。其中，模型的 状态向量包含了与系统相关的信息。如在目标跟踪问题中，状态向量包含了目标 动力学特征参数一其在坐标系中的方位和速度值。状态空问模型包括两个模型： 一是状态方程模型，反映动态系统在输入变量作用下在某时刻所转移到的状态； 二是输出或测量方程模型，它将系统在某时刻的输出和系统的状态及输入变量联 系起来。 状态方程模型： 鼍= 八k - l ，心一1 ) ( 2 1 ) 测量方程模型： 乙= g ( 鼍，吮) ( 2 2 ) 式中，k n 为离散时间k 时刻，鼍r 畋为系统在k 时刻的状态向量， z ：r 以为系统在k 时刻的观测向量，它的维数一般小于状态向量。以尺如、 r 以分别为系统噪声和观测噪声。厂：私水舻专肜、g ：r d 木庐 弘 分别为系统状态转移函数和观测函数，通常为未知的。以、是相互独立的随机 噪声。 引入状态向量是为了对系统内部结构进行数学描述，但在许多情况下，系统 的状态部分是无法直接测量到的，有时甚至全部不能测量到。在实际工作中，通 常能测量到的值只是系统的输入和输出。建立模型的目标就是通过一系列的观测 值( 即输h 值) e , 1 ：k ，z & = 乜，z 2 ，反) 估计出任意时刻的状态向最五。 7 2 2 贝叶斯估计和蒙特卡罗方法 2 2 1 信号处理中的滤波、预测和平滑问题 正如上文指出，建立状态窄间模型进行序贯信号处理的目的是从一系列观测 值z l ：七中估计出状态向量x k 。在处理过程中，有三个概率密度函数扮演了重要 角色，它们是：滤波函数p ( x kz l ：k ) 、预测函数p ( x k 十，lz l ：) ，z 1 、平滑函数 p ( x ki 毛：r ) ，t k 。这三个函数包含了所有关于状态向量义七的信息。当进行 滤波、预测和平滑时，三个概率密度函数可以分别表示为f 1 4 】： 滤波概率密度函数： ，、p ( z ki 垓) 反l z , ：) i 及兹l 兹一。) 烈蕞一。i z w ：纠地一， 聃旧* 户 阿忑赢莨丽五五再 q 3 从式子巾可以看出，) ak - l 时刻的密度函数p ( x k 一。iz l ：k - i ) 计算出k 时刻的密度函 数p ( x 。 z ) 需要二次积分。同样的，预测密度函数、平滑密度函数也需要积分 计算。 预测概率密度函数： p ( 屯+ ，lz i ：女) = ip ( x k + ，i 妖+ ，_ 1 ) p ( 以+ ，1lz l ：k ) d x k + ，- l ，l 1 ( 2 4 ) 平滑概率密度函数： p ( x kiz - ：七，= p c x 七izt：七，兰丛兰量专当舌芝詈鼍挚七+c2 5 ， 很明显，平滑密度函数p ( 吒iz 。：t ) 是从后面时刻的密度函数p ( 稚+ llz 1 ：七) 获得。 信号处理中，概率密度函数严格按照以上函数分布的被称为最优算法。但是 现实中只有很少情况有最优算法，如带高斯噪声的线性模型中使用的卡尔曼滤波 k f 。在非线性或者非高斯系统中通常是使用次优算法来获得近似解。 2 2 2 贝叶斯估计 在统计估计理论的经典方法中，假定所需估计的状态参数x 是确定的未知数， 而当假定状态参数x 为随机变量，估计的是其特定的一个现实时，这被称为贝叶 8 斯( b a y e s i a n ) 方法【1 5 1 。这样做的好处在于：当有一些关于x 的先验知识，就可 以将先验知识运用到估计量中，而这么做正是要求假定x 是具有给定p d f 的随 机变量。这种贝叶斯方法提供了一种与经典方法不同的估计方式，可以利用以前 的信息来估计现在的概率密度函数，获得更精确的估计值。 jp d f、t ：k 2 4 , i ；之一一名41 一，z 】 亨 1 - i 匡舂司一， 卜x 【门】 门= o ，l ，一l 贝叶斯原理 图2 1 贝叶斯办法的数据模型l i ，j f i g u r e2 1t h ed i g i t a lm o d e lo ft h eb a y e s i a nm e t h o d 设x 兰x o ：t 三( x o ，) 为一离散时间序列信号，坼为k 时刻系统的状态向 量，z 兰z o ：t 三( z o ，毛，z k ) 为观测量序列，在信号处理过程巾将通过观测量z o ：k 来 估计出状态向量x o ：k 。在贝叶斯方法中确定后验p d f 为： 舭|z)2鼍铲=瓦p(z一0：kxo：k)p(xo：k) 6 ， 式中分母正好是个与x 无关的归一化因子，它可以保证p ( xiz ) 的积分是 l 。所以在求取后验p d f p ( x i z ) 时，只注意分子即可。 后验p d f 是指得到观测数据后x 的p d f ，与此相对的，p ( x ) 或者式： p ( x ) = ip ( z ，x ) d z ( 2 7 ) 称为先验p d f ，表示得到观测数据之前的p d f 。在贝叶斯估计中先验p d f 的选择 是很关键的，错误的选择将导致差的估计量。 在贝叶斯估计中，观测到数据后，后验p d f p ( x iz ) 概括了我们对这个参数的 了解情况。对状态向量x 的所有现实和观测量z ，使平均均方误差( m e a ns q u a r e e r r o r ，m s e ) 最小的估计量定义为最佳估计量，即所谓的贝叶斯m s e 。该估计量 是后验p d f 的均值，即j = e ( xz ) 。估计量可以显式地表示为： x = e ( xjz ) ：jx p ( x l z ) d x( 2 8 ) 在实际中实现贝叶斯m s e 估计量时，不得不借助数值积分，在矢量情况下相 9 当复杂，后验p d f 变为： p ( 耶) = 警= 丽p ( z l x 丽) p ( x 面) ( 2 9 ) 它需要求x 的p 维积分，求均值j = e ( x1z ) 时需进一步积分。计算是相当 复杂的。 在许多信号处理的应用i 口j 题中，接收的数据是对连续时问信号进行采样而得 到的。随着时间进展，数据源源不断，要么等所有数据采样完成再一次性全部处 理，称为批处理；要么按时间顺序进行处理，称为点处理。 从式( 2 6 ) 中可以看出，在每一个新时刻当观测值到来时，后验p d f 就被更 新一次，直接计算很不方便。因此，首先用贝叶斯方法获得p d f ，再通过迭代更 新它。p ( x o ：tiz o ： ) 可以以序贯形式表示如下： p ( z 0 ：k h ) ：丝姒暨掣 叫| 强一晓1 0 ) 从式子看，这种序贯形式的p d f 接收到新数据进行序列处理而不用对i f l 的数 据苹复计算，相当于点处理，明显比批处理节约时问。同样的，k 时刻状态向量x 的边缘p d f p ( x kiz o ：七) 年- f l k + l 时刻状态向量x 的预测p d f p ( x k + 1iz o ： ) 都可以通 过这种序贯形式计算。贝叶斯信号处理已经成为解决信号检测问题的一个重要上 2 2 3 蒙特卡罗方法 近年来，基于蒙特卡罗方法【8 】的滤波器也引起了人们的广泛关注，其基本思想 是产生服从后验概率密度分布的样本，并对其进行加权给出密度函数的近似。这 类方法大体可分为两类：一种是马尔可夫链的蒙特卡罗方法( m a r k o vc h a i nm o n t e c a r l o ，m c m c ) ，其中马尔可大链给出服从感兴趣的后验概率密度分布的样本，当 有新的观测值到达时，m c m c 过程要重新开始，因而这是一种用于解决批处理问 题的方法；另一种就是序列蒙特卡罗方法( s e q u e n t i a lm o n t ec a r l o ，s m c ) 。从前 文的介绍，我们看到状态向量x 的后验p d f 计算为多维积分，很难直接计算得出。 1 0 这个问题可以通过蒙特卡罗方法解决。 2 0 世纪4 0 年代中期由于科学技术的发展，特别是核武器研制和电子计算机的 发明，科研人员提出了蒙特卡罗方法，在原子能技术中广泛应用，解决了很多典 型数学问题，如重积分、微分方程边值问题、线性方程组求解。蒙特卡罗方法又 称随机抽样法或统计试验方法，它是一种采用统计抽样近似地求解数学、物理及 工程的方法，与一般数值计算方法有本质区别，属于试验数学的一个分支。蒙特卡 罗方法的基础理论是概率统计，其基本手段是随机抽样或随机变量抽样。它首先 建立一个与求解有关的概率模型或随机过程，使它的参数等于所求问题的解，然 后对随机过程抽样实验米计算所求参数的统计特征，最后给出所求解的近似值。 与其它数值计算方法相比，蒙特卡罗方法有如下的优点： 1 ) 蒙特卡罗方法是通过人量简单的重复抽样实验实现的，故计算方法和程序 结构简单，易于在计算机上实现。 2 ) 收敛速度与问题维数无关。换句话说，要达到同一精度，用蒙特卡罗方法选 取的点数与维数无关，计算时问仪与维数成比例。这一特性，决定了埘多维问题的 适用性。 3 ) 适应性强，受问题条件限制的影响小。 蒙特卡罗方法给出了一种求解积分的数值方法，尤其是当无法得到闭式表达 式时，这一方法是十分有用的。假设需要计算的积分为【1 6 】： ，= 昭f ( x ) d x ( 2 1 1 ) ( x ) 在【以6 】上的平均值为_ l ，令x z ，x 2 ，矗为k 6 】中随机抽出的n 个点，则 通过计算( 一) ，f ( x 2 ) ，厂( 矗) 得到f ( x ) 的高度，形成样本均值： 五= 丙1 善n 厂( 薯) ( 2 1 2 ) 我们希望厶= _ l ，因而 ，= rm 陟= 等善他) 当誓采用随机数时，此方法称为蒙特卡罗方法，计算的实质是大数定理。按照上 述方法，后验p d f 的期望 x = e ( x 0 ：kiz 0 ：后) = jx o ：k p ( x o ：k l z o k ) d x 0 ：k ( 2 1 4 ) 可以表示为： 专兰雠2 ) 眨 其巾x o ( n 盔；j = l ，2 一，) 为服从p ( 飞七iz 0 = k ) 的n 个独立同分布的随机样本。当n 足够大时，瓦面j 丽绝对收敛于( 龟量lz 0 ：露) 。 2 3 基本粒子滤波算法 在非线性系统巾，由于实时处理和计算存储量的要求，通常选用序贯贝叶斯 估计来求解非线性估计问题。这种序贯贝叶斯递推滤波器由两个步骤组成：预测 和更新。预测是利用系统模型预测从一个测量时刻到f + 4 个时刻的状态后验概率 密度函数，而更新操作是利用最新的测量值对这个概率密度函数进行修正。 在片j 上述的序贯贝叶斯递推滤波器对后验概率密度函数进行估计时，只能在 某些特殊情况下获得后验概率密度函数的精确解析解。例如，对+ 二线性高斯随机 系统，卡尔曼滤波器能获得后验概率密度函数的精确解析解。但在大多数情况下 获得精确的后验概率密度函数是不可能的，这就需要迸行各种近似次优估计。推 广卡尔曼滤波就是最常用的近似次优估计方法。然而由于系统的可观测性较差， 状态空间模型的非线性程度较高，均导致了推广卡尔曼滤波算法在收敛精度及收 敛时间上往往满足不了要求。粒子滤波方法【3 】通过非参数化的蒙特卡罗模拟方法来 实现贝叶斯估计，用样本形式而不是以函数形式对先验信息和后验信息进行描述。 当样本点数增至无穷大，蒙特卡罗模拟特性与后验概率密度的函数表示等价，从 而滤波精度可以逼近最优估计。在数学处理的方便性上，它不受线性和高斯分布 的限制，原则上适用于能用状态空间模型表示的任何非线性、非高斯系统，因此 受到了高度重视。 2 3 1 序贯重要性采样 粒子滤波正是一种基于蒙特卡罗方法和递推贝叶斯估计的滤波方法，其中 用到的一个核心算法是序贯莺要性采样算法。 序贯重要性采样算法s i s 通过蒙特卡罗方法模拟实现递推贝叶斯滤波，是粒 子滤波的基础。其核心思想是利用一系列随机样本的加权和表示所需的后验p d f ， 1 2 得到状态向量的估计值。 在序贯信号处理过程中，常用状态方程和观测方程表示动态系统，即如式 ( 2 1 ) 、( 2 2 ) 所示。 粒子滤波算法中，使用 础，w ( k i ) 羔l 表示系统后验概率密度函数p ( 黾。i 孙。) 的 样本集合，其中 x 敬，i 1 ，2 ，n ) 是粒了集合，n 为粒子数，各个粒了相应的权值 n 为 以，扛l ，2 ，) ，且满足以) - l ，而x o ：t = ，j = o ，1 ，k 表示时刻0 到 i = l 时刻k 系统所有状态状态向量x 的序列，所以时刻k 的后验概率密度函数町以近 似表示为【1 7 】： n p ( x o ：。l ： ) 以5 ( x o ：。一础) ( 2 1 6 ) ，= l 其中权值可以通过重要性采样法计算得到。 这时很容易估计出我们所需的参数，例如： e ( f ( x o ：。) ) 以。厂( 础) ( 2 1 7 ) f = l 假设概率密度函数是p ( x ) ，如果能直接从中产生粒子x n ，则每一个粒子的 权值是。当无法直接从p ( x ) 中产生粒予时，我们选择从分布万( x ) 中产生粒予 x n ，这里万( x ) 被称为