配套课件-目标跟踪系统中的滤波方法.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,第,1,章 绪 论,1.1 滤波方法在目标跟踪系统中的地位和作用,1.2,状态估计和融合方法的研究进展及现状,1.3,目标跟踪滤波性能评价准则,1.4,本书内容安排,1.1,滤波方法在目标跟踪系统中的地位和作用,目标跟踪是人们运用各种观测和计算手段，实现主体对被关注运动客体的状态建模、估计、跟踪的过程。随着航空、航天、航海事业的不断发展以及现代战争信息化、网络化特征的日益凸显，对海底、海面、陆地、空中和太空中目标跟踪技术的精确性和实时性的要求在不断提高。毋庸置疑，该技术在国防安全领域中发挥着重要作用。此外，在民用领域，目标跟踪技术也得到广泛应用，例如空中交通管制、机器人、视频监控以及存在于制造业中的工件定位等。,一般意义下的目标跟踪技术通常包括三个部分：数据关联、状态估计及融合、航迹管理,1,，如图,1.1,所示。在多目标多传感器测量环境下，数据关联的作用很重要，不正确的数据关联直接导致跟踪精度降低甚至丢失目标。数据关联从所关联数据的类型的角度可做如下划分：量测和量测关联、量测和局部估计关联以及航迹和航迹关联。数据关联的难度体现在传感器存在漏检、虚警以及在测量值比较密集情况下的关联。,Sittler,最早提出数据关联的概念，,Singer,和,BarShalom,对数据关联理论的研究和发展做出过重大贡献。最为重要的数据关联方法包括最近邻方法、概率数据关联、联合概率数据关联、多假设方法以及上述方法的改进算法等。航迹管理包括航迹起始、航迹终止、航迹维持等。,图,1.1,目标跟踪系统,1.2,状态估计和融合方法的研究进展及现状,1.2.1,信息融合技术,信息融合通俗地讲，是将不同来源、不同模式、不同媒质、不同时间、不同地点、不同表示形式的信息进行综合，最后得到对被感知对象的更精确描述。信息融合与数据融合虽属两个不同的概念，但两者很类似，人们所指的信息融合一般是指数据的融合，因此本书对这两个概念不加区分。由于多传感器测量信息具有冗余性和互补性，因此采用有效的融合方法可以得到更加可靠和更加准确的信息。冗余性可以增加系统的健壮性，如果传感器网络中有部分节点遭到破坏，则可以利用其它正常工作的传感器所携带的相同信息进行弥补。互补性可以提高信息的准确性，可以充分利用不同传感器的测量特征，从而获取被观测目标的更准确、更全面的信息。,根据美国国防联合指导实验室,(the US Joint Directors of Laboratories,，,JDL),的定义，信息融合是组合数据的过程，该过程能够进一步优化估计和预测,3,。也就是说，如果能够有效地利用多源信息，则能够在一定程度上获取更加可靠、更加精确的估计和预测信息。为了说明信息融合的特征，我们给出一个典型的信息融合例子，如图,1.2,所示。在目标跟踪过程中，使用雷达,(Radar),、可见光或红外,(EO/IR),和电子支援设备,(ESM),对目标进行观测，然后进行融合。可以看出，不同类型的测量设备在探测性能、运动参数估计和识别能力方面是不完全相同的。也就是说，这些传感器的性能具有互补性，所获得的测量信息也相应地具有互补性，使用有效的融合方法进行融合，能够获得各项指标都较好的总体性能。此外，以各种途径获得的数据之间存在冗余，当一部分数据被破坏或者丢失之后，可以通过其它传感器数据对该部分数据进行恢复，因此提高了系统的健壮性。,图,1.2,信息融合举例,信息融合在军事上涵盖的范围是比较广泛的，从最底层的数据预处理到战场中的指挥决策，全都包含在内。,JDL,实验室对信息融合进行了分层，给出了各层的作用范围框架，如图,1.3,所示。下面对每一层的定义简单加以说明。,图,1.3,信息融合分层模型,3,第,0,层,次目标数据估计：在对有偏的信号或者像素层数据进行数据关联和特征化的基础上，对信号,(,即目标观测向量,),进行估计和预测。第,1,层,目标估计：在有偏的量测和航迹关联的基础上，对实体状态进行估计和预测，包括连续状态估计和离散状态估计。第,2,层,态势估计：估计和预测实体之间的关系，包括军力结构和交叉军力关系、通信和感知影响、自然条件等。第,3,层,影响估计：战争各方的计划或者对其估计,(,或预测,),的行为对态势影响的估计和预测，包括多个参战方行动计划的相互作用等。,第,4,层,过程优化：自适应地获取和处理数据，以便支持战争目的。该处理涉及到计划和控制，不是估计。该层任务主要是根据各层结果进行判断，指定有利于我方的部署，并将任务分配给各种资源，最终达到目前态势和预测态势有利于我方的目的。然而上述分层模型只是一个框架，不够具体。赵宗贵在此基础上结合军事实例对该模型涵盖的内容进行了细化，如图,1.4,所示。信息融合技术涉及的范围尽管很广泛，但是现有的研究成果主要集中在低层,(,也就是第,0,层和第,1,层,),。这主要是因为高层需要具体的实例环境以及相关的军事理论和实践知识，大部分理论研究者缺乏此类专业知识因此，就公开的文献来看，针对高层信息融合的研究开展得相对较少,5,。如果按照,JDL,分层理论，本书的研究范畴属于第,1,层次。,JDL,划分的第,1,层次中一个重要的研究内容就是运动目标跟踪。目标的精确定位和运动估计在现代战场上对决定战场态势具有重要作用。因此，近半个世纪以来，世界各国尤其是发达国家，对信息融合的基础理论研究和工程实践都非常重视。在这期间，出现了众多理论研究成果和专著。,图,1.4,信息融合分层在军事中的实例,4,最具有影响的专著有,Shalom,等人的文献,1,和文献,6,、,Blackman,等人的文献,7,以及,Hall,的文献,8,等。还有许多优秀的相关专著极大地推进了信息融合技术的不断发展,9-10,。我国在上世纪,80,年代中期开始数字雷达处理技术的研究，也出现了多部信息融合领域的著作,11-18,。尤其是,2009,至,2010,年，涌现出多部信息融合和目标跟踪领域的专著,19-25,，这些新著作必将会为信息融合技术在我国科研人员中的普及发挥重要作用。在信息融合领域做出过巨大贡献的有如下几位：,BarShalom Y.,、,Samuel S.Blackman,、,Farina,和李晓榕,(Li X.R.),等人，他们长期活跃在信息融合领域，不断地推动并引导信息融合技术的发展。国内也展开了一系列研究，主要研究机构有上海交通大学、国防科技大学、杭州电子科技大学、海军航空工程学院、西安交通大学、西北工业大学和西安电子科技大学等单位。,信息融合领域的重要刊物包括,IEEE Transactions on Automatic Control,、,IEEE Transactions on Aerospace and Electronic Systems,、,IEEE Transactions on Signal Processing,和,Elsevier,出版集团的,Automatica,等。信息融合国际会议,(International Conference on Information Fusion),是该领域的学术盛会。,1.2.2,目标跟踪技术,使用状态空间估计方法对目标进行跟踪，首先要对目标的运动模型和传感器的测量模型建模。测量模型根据传感器的性质容易确定，然而，目标的运动模型事先难以确定。为了简化问题，人们对常见的理想状态下的目标运动进行建模,26,。简单的目标运动模型包括匀速运动、匀加速运动、零均值一阶马尔可夫模型、均值自适应加速模型、已知转弯角速度的常速度转弯运动和未知转弯角速度的常速度转弯运动等。在基于模型的目标跟踪过程中，对目标运动情况的建模十分重要，模型的好坏直接决定着跟踪的效果。,然而在实际应用中，目标的运动模型往往是未知的，而且在现代战场上，目标在运动过程中会做各种各样的运动模式变换，以摆脱对方的追踪和打击。如果仅仅使用已知的模型进行状态估计，效果就会很差，在这种情况下，人们对机动目标的跟踪产生了浓厚的兴趣。所谓机动目标，是指目标的运动模式不断变换，跟踪系统无法准确知道当前目标的运动模型。李晓榕及其合作者十分关注机动目标跟踪方法的研究和该领域的进展，并且分主题做了相关的综述。,1.2.3,状态估计技术,在目标跟踪过程中，传感器测量信息不仅包含所需信号，同时也包含随机观测噪声和干扰信号。估计是指通过对一系列带有观测噪声和干扰信号的实际观测数据进行处理，从中得到所需要的各种参量的估计值的过程。通常估计问题可以分为两类：一类是系统的结构参数部分或全部未知，有待确定；二是系统中的部分或全部状态变量不能直接测得。这两类问题通常称为参数估计和状态估计。两者的区别在于：参数估计是不随时间变化或只随时间缓慢变化的随机变量；状态估计是随时间变化的随机过程。对于状态估计的国外专著如文献,41-43,，国内也涌现了一系列著作可供参考,44-48,。,根据状态向量和观测向量在时间上存在的不同对应关系，状态估计问题可以分为预测、滤波和平滑。预测是滤波的基础，滤波是平滑的基础，因此下面仅针对主要的滤波算法进行讨论。假设,x,k,|,j,表示根据,j,时刻和,j,时刻以前的测量值对,k,时刻的状态,x,k,作出的估计，则按照,k,和,j,的不同对应关系，状态估计可作如下划分：,(1),当,k,=,j,时的估计过程称为滤波，即依据过去直到现在的观测值来估计现在的状态；,(2),当,k,j,时的估计过程称为预测，即依据过去直到现在的观测值来预测未来的状态；,(3),当,k,0,，,i,=1,，,2,，,，,r,1,，则点,u,为一个生成元。,把全对称集中的生成元,u,=(,u,1,，,u,2,，,，,u,r,，,0,，,，,0),简记为,u,1,，,u,2,，,，,u,r,。若生成元中的各分量互不相同，则生成元的全对称集有个元素。例如，,1R,2,表示二维空间的点集：，该点集的生成元为。我们用,u,1,，,u,2,，,，,u,r,i,表示全对称集,u,1,，,u,2,，,，,u,r,中的第,i,个点。,2.,球面容积规则,根据不变量理论，假定一个三阶球面容积规则的最简单形式为,(2-139),由,u,i,经过置换和符号变换得到的点集是不变量。对于单项式，若为奇数，则式,(2-139),式可以准确积分。对于所有三阶以下的单项式，只需要考虑单项式处于,2,的情况。因此，对于单项式,f,(,y,)=1,和,f,(,y,)=,y,2,1,，需要找出未知参数,u,和,w,，使其满足全对称容积规则,(2-140),(2-141),当,f,(,y,)=1,时，,当,f,(,y,)=,y,2,1,时，,3.,径向规则,接下来介绍如何利用高斯求积分法求解径向积分。高斯积分法是已知计算一维积分最为有效的数值积分方法。m,点的高斯积分可以很好地逼近,2,m,1,阶多项式，即,(2-142),其中，,w,(,x,),是,a,，,b,上已知的非负权函数，未知积分点,x,i,和相应的权值,w,i,具有唯一的解。比较式,(2-136),和式,(2-142),，我们得到区间,0,，,的权函数为,w,(,x,)=,x,n,1exp(,x,2,),。为了把这种积分转变成熟悉的形式，做变量代换，令,t,=,x,2,，则有,(2-143),对于一阶高斯规则，有,f,(,t,)=1,，,t,;,f,(,x,)=1,，,x,2,。但是它对于奇数阶多项式,f,(,x,)=,x,，,x,3,并不适用。但若把径向规则与球面规则相结合求解式,(2-135),，则所有的奇数阶多项式的积分为零。这是因为对于球面规则,(,如式,(2-136),，所有奇数阶多项式在对称空间的积分为零，因此一阶通用高斯,-,拉格朗日规则数值积分需要一个点和相应的权值。因此有,(2-144),4.,球面径向规则,为了将球面规则和径向规则结合起来，并将球面径向规则用于高斯加权积分，下面需要介绍两个有用的定理。,定理,2.2,利用,m,r,点高斯求积分规则对径向积分进行数值化计算,(2-145),利用,m,s,点球面规则对球面积分进行数值化计算,(2-146),于是得到一个,m,s,m,r,点球面径向积分规则,(2-147),证明：,求容积法规则可以设计成适用单项式的子空间，考虑如下被积函数,(2-148),其中,d,i,是正整数。因此，有如下积分,(2-149),此时，假定上式中的被积函数是一个,d,阶单项式，即 ，作球面径向变换得,(2-150),把上式分解成为一个径向积分和一个球面积分,(2-151),运用数值规则，得到,(2-152),证毕。,对于被积函数是阶数不超过,d,的任何单项式，上述定理均适用。,定理,2.3,令函数,w,1,(,x,),和,w,2,(,x,),分别为,w,1,(,x,)=exp(,x,T,x,),和,w,2,(,x,)=,N,(,x,;,，,),，对于满足的平方根矩阵，可以得到,(2-153),证明：考虑式,2-153),左侧。因为,是一个正定矩阵，故做矩阵分解，令，可以得到,(2-154),证毕。,对三阶的球面径向规则，,m,r,=1,和,m,s,=2,n,，共需要,2,n,个容积点。根据定理,2.2,和定理,2.3,，可用三阶球面径向规则计算标准的高斯加权积分，如下式,(2-155),其中,(2-156),使用求容积点,i,，,w,i,数值化计算积分，就得到了,CKF,算法。需要指出的是，上面的容积点定义在笛卡尔坐标系。,2.6.2,容积卡尔曼滤波,应用统计线性回归方法，可得容积卡尔曼滤波算法的状态更新步骤和量测更新步骤如下,12,，,13,。预测步骤为：,Step 1,：采用,Cholesky,或者奇异值方法分解协方差,(2-157),Step 2,：计算容积点,(,i,=1,，,2,，,，,m,),(2-158),其中,m,=2,n,x,。,Step 3,：计算预测状态容积点,(2-159),Step 4,：估计预测状态,(2-160),Step 5,：估计预测协方差,(2-161),更新步骤为：,Step 1,：预测协方差矩阵进行分解,(2-162),Step 2,：计算容积点,(,i,=1,，,2,，,，,m,),(2-163),Step 3,：计算预测量测容积点,(2-164),Step 4,：估计预测量测,(2-165),Step 5,：估计新息协方差矩阵,(2-166),Step 6,：估计互协方差矩阵,(2-167),Step 7,：计算卡尔曼滤波增益,(2-168),Step 8,：估计状态更新,(2-169),tep 9,：估计状态协方差,(2-170),2.7,傅立叶厄米特卡尔曼滤波算法,2.7.1,傅立叶厄米特级数展开,本节介绍傅立叶厄米特级数如何应用于高斯随机变量的非线性变换，以及如何看待它作为统计线性化的一般形式。对于随机变量,x,N,(,m,，,P,),，考虑统计线性化状态模型,y,=,g,(,x,),(2-171),假设其近似化形式为,g,(,x,),b,+,A,x,(2-172),通过最小化表达式,E,g,(,x,),b,Ax,2,获得最优矩阵,A,和向量,b,，其中,x,=,x,m,。我们试图将其一般化，使得它能进行线性化近似。使用一个,p,阶多项式进行近似,(2-173),通过扩展上述表达式并且令其导数为零，就可以确定最优多项式的系数。该方法对于高阶多项式的求解就会变得复杂。所幸的是，依据希尔伯特空间理论可以让多项式的求解变得简单。通过给标量函数,g,和,f,定义内积,(2-174),其中,x,N,(,m,，,P,),。通过定义下列范数确定一个希尔伯特空间函数,(2-175),希尔伯特空间理论表明，存在一个与希尔伯特空间对应的正交多项式基。结果证明这些多项式基函数是尺度多元厄米特多项式，定义为,(2-176),其中,L是一个满足,P,=,LL,T,的矩阵，且有,(2-177),这样，就可以将一个满足条件,g,i,，,g,i,的任意向量函数,g,(,x,),扩展为一个傅立叶厄米特级数，形式为,(2-178),需要指出的是，因为,H,0,，,，,0,(,x,;,m,，,P,)=1,，级数中的零阶项恰好就是函数的期望,E,g,(,x,),。利用基函数的正交性可以得到外积期望之和,(2-179),在不考虑零阶级数的情况下，我们可以得到,g,(,x,),协方差的表达式,(2-180),由希尔伯特空间理论，可知对函数,g,(,x,),关于,|,2,的最优,p,阶多项式近似，也即对式,(2-173),的多项式展开，可以通过,p,阶厄米特多项式的正交投影得到。那么通过对级数,(,式,(2-178),进行,p,阶级数截断就可以得到最优,p,阶多项式近似。,现在考虑使用数值计算的方法来求取级数的系数，然后选取其中的零阶项作为均值，并通过上述对级数截断的方法计算协方差的近似值。然而，用这种方法计算协方差并不是一个好的选择，因为使用同样的数值方法我们可以直接获取协方差的值。实际上，依据下述结果，可以采用一种替代的方法获取傅立叶厄米特级数系数，即,(2-181),上述公式是统计线性化方法导数的一般形式。由此，傅立叶厄米特级数就可以用下面的形式表达,(2-182),且其协方差为,(2-183),结果证明我们甚至不需要计算系统导数的期望值，因为可以使用以下方法求解。,假定对于任意的,m,，,P,，使用下列形式来计算期望值,(2-184),接下来，有,(2-185),为了获得扩维函数,g(,x,)=(,x,，,g,(,x,),的联合均值和协方差，同样使用上述结果，并计算,(2-186),其中 。从而得到下面的具有加性变换的傅立叶厄米特级数近似规则。,对于,x,N(,m,，,P,),，,q,N(0,，,Q,),，随机变量,x,和非线性变换,y,=,g,(,x,)+,q,的联合分布的傅立叶厄米特级数为,(2-187),其中，,(2-188),且有如下定义,(2-189),(2-190),矩阵,G,定义如下,(2-191),当级数在截断时，上述近似形式的均值及互协方差总是很精确的，协方差也可以精确到,p,阶多项式。一阶近似等同于统计线性化近似方法。,用与上述类似的方法也可以获得非加性变换的傅立叶厄米特近似过程。,2.7.2,傅立叶厄米特卡尔曼滤波,通过傅立叶厄米特级数展开，我们得到了傅立叶厄米特卡尔曼滤波算法,14,。该算法是统计线性化滤波的一般化形式。为了实现这种滤波算法，需要下述表达式的近似形式,(2-192),(2-193),下面的式子同样需要近似形式,(2-194),(2-195),带有加性噪声的傅立叶厄米特卡尔曼滤波的状态预测步骤和更新步骤如下,:,预测步骤为,(2-197),更新步骤为,(2-198),(2-199),(2-200),(2-201),(2-202),其中矩阵,H,的元素定义为,(2-203),2.8,中心差分卡尔曼滤波算法,2.8.1 Stirling,插值公式,Stirling,插值可以解决任意函数的多项近似问题。首先考虑单变量函数情况，然后再将其扩展到多维情况。如果函数,f,是解析函数，将其围绕着,x,=,x,做泰勒级数展开,(2-204),一种最常用的近似方法就是对上述级数进行有限项的截断。保留的项数越多，得到的近似也就越精确。类似地，可以通过插值公式来获取多项式的近似形式。这些插值公式通常不要求对函数求导，而是对有限个函数进行计算。因此，该方法可能会更容易获取非线性函数的线性近似形式。下面考虑一种称为,Stirling,的差值公式。,首先定义两个操作符,，,如下,(2-205),其中,h,是选定的区间长度。在,x,=,x,处，采用,Stirling,插值公式进行多项式逼近，得到,(2-206),(2-207),通常,1,p,1,。,考虑到我们要着重介绍一阶和二阶多项式逼近，因此将式,(2-207),简化为,(2-208),其中,(2-209),(2-210),式,(2-208),可以解释为把求导项替换为中心差分项的泰勒近似。为了评估多项式的逼近程度，用展开的泰勒级数代替和，得到,(2-211),可以看出式,(2-211),右边的前,3,项和区间长度h,无关。现在考虑多维的情况。设,x,为一个向量，,x,R,n,，,y,=,f,(,x,),，有多种不同的方法可以将插值公式扩展成为多维的，但是在这之前，我们首先想到多维泰勒级数展开形式。当,x,=,x,时，函数,f,的泰勒级数展开为,(2-212),其中，算子,D,i,xf,表示为,(2-213),上述算子同样可以表示为,(2-214),那么，采用多维插值公式对函数进行二阶多项式截断可得,(2-215),其中，差分算子,D,x,f,和,D,2,xf,可分别表示为,(2-216),(2-217),式中,p,是局部插分算子，且有,(2-218),式中,e,p,表示第,p,个元素为,1,的单位向量。均值运算符,可以做类似扩展。,式,(2-215),仅仅是一个插值公式向多维扩展的例子，为了说明如何推导其它形式的扩展，首先介绍对,x,的线性化变换,z,=,S,1,x,(2-219),函数,f,定义如下,(2-220),函数,f,和,f,有着相同的泰勒级数逼近形式。而采用式,(2-215),多项式插值公式，,f,，,f,显然得不到相同的结果。这是因为,(2-221),2.8.2,中心差分逼近,随机向量,x,的均值和协方差分别定义为,(2-222),(2-223),然后再定义,(2-224),(2-225),由于,f,是非线性函数，不可能计算出其精确的均值，通常情况下采用一阶或二阶多项式逼近形式代替。此处，我们用插值公式,(2-215),来获取,f,的多项式近似形式。首先把协方差矩阵的,Cholesky,分解因子作为变换矩阵,(2-226),该变换有时也是对,x,进行随机解耦，使得,z,的各个元素之间相互无关，即,(2-227),1.,一阶近似,首先对函数使用一阶截断插值公式进行逼近,(2-228),根据定义有,E,z,=0,，则式,(2-228),的期望为,(2-229),如前所述，由于,z,均值为,0,，因此一阶项可以忽略。此外，有,E,z,i,z,T,j,=0,，,i,j,。一阶估计协方差即式,(2-224),推导为,(2-230),考虑到,，式中,s,x,，,p,表示对协方差矩阵,S,x,进行,Cholesky,分解后的因子的第,p,列，则式,(2-230),可以表示为,(2-231),同样，我们可得到一阶估计互协方差矩阵为,(2-232),上式也同样可以写为,(2-233),可以看出，区间长度对于一阶估计均值没有影响，但是却对一阶协方差和互协方差有影响。分析表明，区间长度,h,最优值的选取受到,z,分布的影响。可证明,h,与其峰值存在关系,h,2,=,4,。,2.,二阶近似,采用二阶截断插值公式，可以得到,f,的更加精确的均值和协方差的逼近值。其二阶截断插值公式为,(2-234),由于,z,的均值为零，且其各个元素是互不相关的，可以得到,f,的二阶估计期望为2-235),上式可被写为,(2-236),可以获取协方差的估计值为,(2-237),那么，可以得到二阶协方差估计值为,(2-238),式,(2-238),中的奇数矩阵均值为零，它的第一项已经在一阶近似中研究过，现在就来看它的第二项和第三项。第二项,由如下三项构成,(2-239),(2-240),(2-241),第三项 由如下两项构成,(2-242),(2-243),式,(2-240),和式,(2-243),是相等的，两者相互可以抵消。另外，考虑到随着向量,z,维数的增加，式,(2-241),的计算量也会大大增加，为了降低计算量，将其忽略。再者，也是因为利用二阶插值公式无法计算出所有的四阶矩，就只能用三阶矩逼近函数,f,。因此，二阶协方差近似式为,(2-244),考虑到,2,=1,，且令,h,2,=,4,(,对于高斯分布,4,=3),，上式可以表示为,(2-245),由于,(2-246),即,4,2,2,对任何分布都存在，因此可以选择,h,2,1,。显然，这意味着协方差估计总是半正定的。,二阶互协方差估计,P,xy,的推导为,(2-247),2.8.3,中心差分卡尔曼滤波,考虑非线性系统模型,(2-248),(2-249),假定,v,k,，,w,k,是独立同分布的，并独立于当前和过去的状态，且有,v,k,(,v,k,，,Q,(,k,),，,w,k,(,w,k,，,R,(,k,)状态向量,x,k,的一步预测的均值和协方差定义为,(2-250),(2-251),其中，,Y,k,1,是一个包含有过去量测值的矩阵，即,(2-252),为方便起见，状态估计的更新假定为线性的。令状态估计的估计误差最小化，可得,(2-253),(2-254),其中,(2-255),(2-256),(2-257),相应地，状态向量更新后的协方差矩阵为,(2-258),采用一阶多项式近似方法，得到的状态转移和量测方程如下,(2-259),(2-260),其中,(2-261),(2-262),在上述近似方程的基础上，可以获得下列滤波步骤。,时间更新步骤为,(2-263),(2-264),(2-265),量测更新步骤为,(2-266),(2-267),(2-268),在接下来的部分，将会介绍采用插值公式的方法来获取非线性系统的状态估计值。,1.,一阶,CDKF,下面我们将运用之前所介绍的一阶中心差分逼近方法来获取非线性系统的状态估计值。首先，引入以下四个,Choleskey,分解,(2-269),令表示的第,j,列，其它类似向量也采用这样的方法表示。四个包含差分因子的矩阵定义为,(2-270),(2-271),(2-272),(2-273),考虑一个包含有状态向量和过程噪声的扩维状态向量为,(2-274),由于假定过程噪声独立于状态向量，那么,的协方差为,(2-275),是随机解耦的，引入,z,，使得。那么，不难看出状态估计问题可以映射到如前所述的一般化的向量函数,f,(,z,),上。,使用式,(2-229),进行状态向量的一步更新,(2-276),在式,(2-231),的基础上，通过应用式,(2-270),式,(2-273),定义的矩阵，协方差更新为,(2-277),由于,v,k,、,x,k,相互独立，协方差更新可以写为两个矩阵积之和的形式。,上述计算过程中，协方差矩阵可能会出现非对称和非正定性，从而带来数值计算问题。与卡尔曼滤波类似，通常采用因子分解的方法加以解决。由于式,(2-277),的计算是两个二次项的和，因而不会出现数值计算问题。但是，在量测更新过程中需要用到分解因子，通过对下面复合矩阵的,Choleskey,分解获取分解因子,(2-278),由于该矩阵是矩形矩阵，为了后面的使用，需要将其转换为方阵,Choleskey,因子，这可以通过,Householder,三角化方法获得。量测估计值为,(2-279),相应的复合矩阵为,(2-280),式,(2-280),是式,(2-279),误差协方差的,Cholekey,因子，即,(2-281),对于矩阵，需要采用,Householder,三角化方法将其都转换为方阵。,为了计算状态向量和量测估计之间的互协方差矩阵，使用式,(2-233),中的结果得,(2-282),依据式,(2-253),，卡尔曼增益为,(2283),依据式,(2-254),，状态向量更新为,(2-284),后验协方差可以通过式,(2-258),更新，然而为了避免可能出现的数值计算问题，可以直接使用,Cholekey,分解因子进行更新，则有下式,(2-285),因此，量测更新可以表达为,(2-286),显然协方差矩阵的,Cholekey,分解因子方阵可以通过复合矩阵的三角化分解获得,(2-287),2.,二阶,CDKF,二阶中心差分滤波器可以通过前面二阶近似中的均值和协方差估计得到。首先定义四个包含有中心差分因子的矩阵,(2-288),(2-289),(2-290),(2-291),如同一阶中心差分滤波处理过程一样，我们可以通过式,(2-236),来获取状态估计方程,(2-292),其中，,n,x,表示状态向量的维数，,n,v,表示过程噪声向量的维数。先验协方差的,Cholekey,因子同样可以通过以下复合矩阵获取,(2-293),使用和一阶中心差分滤波相同的方法，量测估计及其协方差分别为,(2-294),(2-295),其中，,n,w,表示量测噪声向量的维数。,互协方差矩阵为,(2-296),卡尔曼增益为,(2-297),状态向量的量测更新为,(2-298),依据式,(2-286),协方差矩阵，可得状态向量后验更新的协方差矩阵为,(2-299),显然其含有,Cholekey,分解,(2-300),2.9,小 结,本章在卡尔曼滤波的基础上，讨论处理非线性系统滤波算法，主要涉及一些常用算法和国际上近年来提出的与卡尔曼滤波相关的非线性滤波算法。讨论的算法有：扩展卡尔曼滤波、不敏卡尔曼滤波、积分卡尔曼滤波、容积卡尔曼滤波、傅立叶厄米特卡尔曼滤波、中心差分滤波。粒子滤波算法也是一类重要的非线性系统滤波方法，放在下一章讨论。,第,3,章 粒子滤波方法,3.1 引言,3.2,贝叶斯滤波,3.3,贝叶斯重要性采样,3.4,序贯重要性重采样粒子滤波算法,3.5,马尔可夫链蒙特卡罗粒子滤波算法,3.6,辅助粒子滤波算法,3.7,正则化粒子滤波算法,3.8,边缘粒子滤波算法,3.9,扩展卡尔曼粒子滤波算法,3.10,高斯和粒子滤波算法,3.11,小结,3.1,引 言,粒子滤波是指通过寻找一组在状态空间中传播的随机样本，对概率密度函数,p,(,x,k,|,y,k,),进行近似，以样本均值代替积分运算，从而获得状态的最小方差估计的一种算法。粒子滤波算法依据系统状态向量的先验分布在状态空间中产生一组随机样本，然后根据观测量不断地调整粒子的权值和位置，通过调整后的粒子信息修正最初的后验概率函数。用数学语言可描述为：针对平稳的时变系统，假定,k,1,时刻系统的后验概率密度为,p,(,x,k,1,|,y,k,1,),，依据一定规则选取,N,个随机样本点，在,k,时刻获得量测信息后,经过状态更新和时间更新过程,N,个粒子的后验概率密度近似为,p,(,x,k,|,y,k,),随着粒子数的增加，粒子的概率密度函数就能逼近状态真实的概率密度函数，对状态向量的估计结果与最优贝叶斯估计结果接近。粒子滤波适用于非线性非高斯系统的状态估计，突破了传统卡尔曼滤波理论框架，精度可以逼近最优估计，是一种有效的非线性滤波技术，广泛应用于数字通信、图像视频处理、计算机视觉、语音信号处理、机器学习等领域。,3.2,贝 叶 斯 滤 波,对于跟踪问题，目标状态序列的状态转移可以用下面的目标状态序列,x,k,k,N,的演变方程来描述,(3-1),其中，,是关于状态,x,k,1,的非线性函数，是独立同分布的过程噪声序列，,n,x,n,v,分别是状态和过程噪声向量的维数。系统量测方程为,(3-2),从贝叶斯估计观点来看，跟踪问题就是计算,k,时刻状态,x,k,的某种置信程度。从,1,到,k,时刻，由于给定量测数据,z,1:,k,的值不同，得出的,x,k,值也不同，因此需要构造概率密度函数,p,(,x,k,|,z,1:k,),。假定初始概率密度函数,p,(,x,0,|,z,0,),p,(,x,0,),，,x,0,表示初始状态向量，,z,0,表示尚且没有获得量测值，它也是先验概率密度函数。因此从形式上看，通过预测和更新两个步骤就可以递推地得到概率密度函数,p,(,x,k,|,z,1:,k,),的值。假定,k,1,时刻的概率密度函数已知，那么通过,ChapmanKolmogorov,等式及使用模型,(3-1),就可以预测出,k,时刻状态的先验概率密度函数,(3-3),在等式,(3-3),中，,p,(,x,k,|,x,k,1,)=,p,(,x,k,|,x,k,1,，,z,1:,k,1,),，且满足等式,(3-1),所描述的一阶马尔可夫过程。状态估计,p,(,x,k,|,x,k,1,),的概率模型由系统等式,(3-1),和统计值,v,k,1,来确定。在,k,时刻，可以得到量测值,z=,，然后通过贝叶斯规则更新先验概率密度函数,(3-4),式中的常量,(3-5),是由模型,(3-2),定义的似然函数,p,(,z,k,|,x,k,),和一步预测的统计值,p,(,x,k,|,z,1:k,1,),共同确定的。在更新式子,(3-4),中，量测值,z,k,被用来修正先验概率，以获得当前状态的后验概率。,式,(3-3),和式,(3-4),是最优贝叶斯估计的一般形式。通过以上步骤递推得到的后验概率只是一般概念下的表达式，通常情况下难以得到其解析表达式。只有在满足特定条件时，才可以得到最优贝叶斯解。,3.3,贝叶斯重要性采样,在贝叶斯重要性采样中，后验概率分布由一组离散的样本集近似得到。根据大数定理，随着样本粒子数,N,的增加，期望,E,g,(,x,0:,k,),可由近似求出。但是，通常很难从后验概率密度函数中直接抽样。常规的解决办法是从容易采样的概率分布,q,(,x,0:,k,|,y,1:,k,),中采样粒子，由此可以得到,(3-6),式中,w,(,x,0:,k,),是未归一化的重要性权值，可以表示为,(3-7),由于,p,(,y,1:,k,),是未知的，我们可以将式,(3-6),表示为,(3-8),式中,E,q,(|,y,1:,k,),是指在概率分布,q,(|,y,1:,k,),上进行计算的期望。通过从概率函数,q,(|,y,1:,k,),中采样，期望可以近似表示为,(3-9),式中的重要性权值表示为,(3-10),等式,(3-9),计算出来的结果是有偏的。但是，通过以下两个假设可以使逐渐收敛，并且接近于真实值。,(1),x,i,),0:,k,是从后验概率分布中采样得到的一组粒子，,E,g,(,x,0:,k,),存在并且是有限的；,(2),在后验概率分布上计算出的,w,k,和,w,k,g,2,(,x,0:,k,),的期望存在而且有限。只有,g,(,x,0:,k,),的方差和重要性权值有限才能验证第二个假设的成立性。随着,N,值的无限增大，后验概率分布函数就会近似于点估计分布，即,(3-11),3.4,序贯重要性重采样粒子滤波算法,3.4.1,序贯重要性采样,贝叶斯重要性采样,(Sequential Importance Sampling,，,SIS),2,3,是一种简单常用的蒙特卡罗积分方法，但是它不能直接用来做递推估计。这主要是因为估计,p,(,x,0:,k,|,y,1:,k,),需要用到所有的观测数据,y,1:,k,，每次更新观测数据,y,k,+1,时，需要重新计算整个状态序列的重要性权值，因此它的计算量随着时间的推移而不断增加。为了解决该问题，人们提出了序贯重要性采样方法。该方法在,k,时刻采样时不改动过去的状态序列,x,0:,k,1,而采用如下递推形式计算重要性权值,(3-12),这里先假设当前状态不依赖于将来的观测，即只进行滤波而不考虑平滑。需要强调的是，在某些情况下，一些建议分布需要用到过去的状态序列。在本书中，不考虑这种情况。假设状态符合马尔可夫过程，在给定状态下，量测值是条件独立的，则可得,(3-13),(3-14),将式,(3-12),、式,(3-13),和式,(3-14),代入式,(3-7),中得到权值递推公式,(3-15),在给定合适的重要性分布函数,q,(,x,k,|,x,0:,k,1,，,y,1:,k,),的条件下，式,(3-15),提供了一个递推计算重要性权值的方法，重要性权值的计算因此得以简化。,3.4.2,序贯重要性采样问题及策略,1.,选取好的重要性密度函数,选择合适的重要性密度函数是重要性采样算法中的关键步骤。选择重要性密度函数的一个原则是使得权系数的方差最小。,Doucet,等提出了在给定,x,0:,k,1,和,y,1:,k,的条件下权系数方差最小的