卡尔曼滤波

卡尔曼滤波
卡尔曼滤波

经典卡尔曼滤波 动态系统的卡尔曼滤波数学模型包括状态方程和观测方程,对于线性系统,其离散形式为 X k = Φk ,k?1X k?1+ G k?1W k?1 L k = H k X k + V k 其中,X k 为系统在t k 时刻的n ×1维状态向量,Φk ,k?1为系统从t k?1时刻到t k 时刻的n ×n 维状态转移矩阵,W k?1为t k?1时刻的r ×1动态噪声,G k?1为t k?1时刻的n ×r 动态噪声矩阵,L k 为系统在t k 时刻的m ×1维观测向量,H k 为系统在t k 时刻的m ×n 维观测矩阵,V k 为系统在t k 时刻的m ×1维观测噪声。 根据最小二乘原则,可推得卡尔曼滤波递推公式为: 1)状态向量一步预测值为 X k ,k ?1 =Φk ,k?1X k ?1,k ?1 2)状态向量一步预测值方差矩阵为 P k ,k ?1 = Φk ,k?1P k ?1,k ?1 Φk ,k?1T +G k?1Q k?1G k?1T 其中,Q k 为动态噪声方差矩阵。 3)状态向量估计值为 X k ,k = X k ,k ?1 +J k L k ? H k X k ,k ?1 4)状态向量估计值方差矩阵为 P k ,k = I ? J k H k P k ,k ?1 其中,J k 为滤波增益矩阵,具体形式如下 J k = P k ,k ?1 H k T H k P k ,k ?1 H k T + R k ?1 R k 为观测噪声方差矩阵。 基于方差补偿的自适应卡尔曼滤波 假定 W k 和 V k 是正态序列,X 0是正态向量。则定义l 步预测残差是 V k +l = L k +l ? L k +l ,k 其中,L k +l ,L k +l ,k 分别为第k +l 期观测值和它的最佳观测值L k +l ,k = H k +l Φk +l ,k X k ,k ,则V k +l 为正态向量。V k +l 的方差阵S vv 为: S vv = H k +l Φk +l,k P k Φk +l ,k T H k +l T +R k +l + H k +l Φk +l ,k +i l i =1 G k +l ,k +i?1Q k +i?1G k +l ,k +i?1T Φk +l ,k +i T H k +l T

文件一 % THIS PROGRAM IS FOR IMPLEMENTATION OF DISCRETE TIME PROCESS EXTENDED KALMAN FILTER % FOR GAUSSIAN AND LINEAR STOCHASTIC DIFFERENCE EQUATION. % By (R.C.R.C.R),SPLABS,MPL. % (17 JULY 2005). % Help by Aarthi Nadarajan is acknowledged. % (drawback of EKF is when nonlinearity is high, we can extend the % approximation taking additional terms in Taylor's series). clc; close all; clear all; Xint_v = [1; 0; 0; 0; 0]; wk = [1 0 0 0 0]; vk = [1 0 0 0 0]; for ii = 1:1:length(Xint_v) Ap(ii) = Xint_v(ii)*2; W(ii) = 0; H(ii) = ‐sin(Xint_v(ii)); V(ii) = 0; Wk(ii) = 0; end Uk = randn(1,200); Qu = cov(Uk); Vk = randn(1,200); Qv = cov(Vk); C = [1 0 0 0 0]; n = 100; [YY XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V); for it = 1:1:length(XX) MSE(it) = YY(it) ‐ XX(it); end tt = 1:1:length(XX); figure(1); subplot(211); plot(XX); title('ORIGINAL SIGNAL'); subplot(212); plot(YY); title('ESTIMATED SIGNAL'); figure(2); plot(tt,XX,tt,YY); title('Combined plot'); legend('original','estimated'); figure(3); plot(MSE.^2); title('Mean square error'); 子文件::function [YY,XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V); Ap(2,:) = 0; for ii = 1:1:length(Ap)‐1 Ap(ii+1,ii) = 1;

卡尔曼滤波算法总结-标准化文件发布号:(9556-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

2015.12.12 void Kalman_Filter(float Gyro,float Accel) { Angle+=(Gyro - Q_bias) * dt; Pdot[0]=Q_angle - PP[0][1] - PP[1][0]; Pdot[1]= - PP[1][1]; Pdot[2]= - PP[1][1]; Pdot[3]=Q_gyro; PP[0][0] += Pdot[0] * dt; PP[0][1] += Pdot[1] * dt; PP[1][0] += Pdot[2] * dt; PP[1][1] += Pdot[3] * dt; Angle_err = Accel - Angle; PCt_0 = C_0 * PP[0][0]; PCt_1 = C_0 * PP[1][0]; E = R_angle + C_0 * PCt_0; K_0 = PCt_0 / E; K_1 = PCt_1 / E; t_0 = PCt_0; t_1 = C_0 * PP[0][1]; PP[0][0] -= K_0 * t_0; PP[0][1] -= K_0 * t_1; PP[1][0] -= K_1 * t_0; PP[1][1] -= K_1 * t_1; Angle += K_0 * Angle_err; Q_bias += K_1 * Angle_err; Gyro_x = Gyro - Q_bias; }

首先是卡尔曼滤波的5个方程: (|1)(1|1)() X k k AX k k Bu k -=--+(1)先验估计 (|1)(1|1)'P k k AP k k A Q -=--+(2)协方差矩阵的预测 ()(|1)'/(|1)')Kg k P k k H HP k k H R =--+(3)计算卡尔曼增益 (|)(|1)()(()(|1))X k k X k k Kg k Z k HX k k =-+--(4)进行修正 5个式子比较抽象,现在直接用实例来说: 一、卡尔曼滤波第一个式子 对于角度来说,我们认为此时的角度可以近似认为是上一时刻的角度值加上上一时刻陀螺仪测得的角加速度值乘以时间,因为d dt θω=?,角度微分等于时间的微分乘以角速度。但是陀螺仪有个静态漂移(而且还是变化的),静态漂移就是静止了没有角速度然后陀螺仪也会输出一个值,这个值肯定是没有意义的,计算时要把它减去。 由此我们得到了当前角度的预测值Angle Angle=Angle+(Gyro - Q_bias) * dt; 其中等号左边Angle 为此时的角度,等号右边Angle 为上一时刻的角度,Gyro 为陀螺仪测的角速度的值,dt 是两次滤波之间的时间间隔,我们的运行周期是4ms 或者6ms 。 同时 Q_bias 也是一个变化的量。 但是就预测来说认为现在的漂移跟上一时刻是相同的,即 Q_bias=Q_bias 将上面两个式子写成矩阵的形式 1_0 1_0 Angle dt Angle dt Q bias Q bia o s Gyr -= + 得到上式,这个式子对应于卡尔曼滤波的第一个式子 (|1)(1|1)() X k k AX k k Bu k -=--+ (|)(|1) P k k I Kg k H P k k =--(())(5)更新协方差阵

卡尔曼滤波简介及其算法实现代码(C++/C/MATLAB) 卡尔曼滤波器简介 近来发现有些问题很多人都很感兴趣。所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。现在先讨论一下卡尔曼滤波器,如果时间和能力允许,我还希望能够写写其他的算法,例如遗传算法,傅立叶变换,数字滤波,神经网络,图像处理等等。 因为这里不能写复杂的数学公式,所以也只能形象的描述。希望如果哪位是这方面的专家,欢迎讨论更正。 卡尔曼滤波器– Kalman Filter 1.什么是卡尔曼滤波器 (What is the Kalman Filter?) 在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”。跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人! 卡尔曼全名Rudolf Emil Kalman,匈牙利数学家,1930年出生于匈牙利首都布达佩斯。1953,1954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器,正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法)。如果对这编论文有兴趣,可以到这里的地址下载: https://www.360docs.net/doc/eb17613470.html,/~welch/media/pdf/Kalman1960.pdf。 简单来说,卡尔曼滤波器是一个“optimal recursive data processing algorithm(最优化自回归数据处理算法)”。对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。 2.卡尔曼滤波器的介绍 (Introduction to the Kalman Filter) 为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是,他的5 条公式是其核心内容。结合现代的计算机,其实卡尔曼的程序相当的简单,只要你理解了他的那5条公式。 在介绍他的5条公式之前,先让我们来根据下面的例子一步一步的探索。

处理线性滤波以及预测问题的一种新途径R.E.Kalman1引言通讯与控制中的理论与实际问题中有很重要的一类具有统计性质。 这样的问题有:(1)、随机信号的预测;(2)、从随机噪声中分离随机信 号;(3)、在有噪声的情况下探测已知形式的信号(脉冲、正弦波)。 在Wiener开拓性的工作中,他证明[1]从问题(1)和问题(2)可导 出所谓Wiener-Hopf积分方程;他同样给出了解决具有实际重要意义的特 殊情况——定态统计和有理数频谱——之积分方程的方法(频谱因式分解)。 在Wiener的基础性工作之后出现了许多延伸和推广。Zadeh与 Ragazzini给出了有限存储器情况的解[2]。Bode和Shannon[3]同时独 立的给出上述情况的解,并且给出了简化的求解方法。Booton讨论了非定 态统计Wiener-Hopf方程[4]。这些结果现在都写入了标准教科书中[5,6]。 最近Darlington[7]沿着这些主线给出了一种稍微有些不同的方法。对抽 样信号的延伸,参见Franklin[8]和Lees[9]的工作。基于Wiener-Hopf方 程(同样应用于非定态问题,尽管前述方法一般并非如此)特征函数的 方法由Davis[10]开创并被许多其他人应用,例如Shinbrot[11],Blum[12], Pugachev[13],Solodovnikov[14]. 在所有这些工作中,目标都是获取一个线性动力系统的明确说明 (Wiener滤波器),由此可以完成预测、分离或者探测随机信号。 现有求解Wiener问题的方法受制于若干限制,这样就使得它们的实际 用处收到削弱: 1.最佳的滤波器由其脉冲响应具体指定。由这些数据合成滤波器并非易 事。 2.数值确定最佳的脉冲响应往往十分复杂并且不很适合机器计算。这种 情况随着问题复杂度的增加而迅速变得更为糟糕。 11引言2 3.重要的推广(如增长存储器滤波器、非定态预测)需要新的推导过 程,经常给非专业人士带来相当大的困难。 4.这些推导过程的数学部分并不透明。基本假设及其后果趋于模糊。 本文回避上述困难,提出看待这些问题的整个集合的新方式。以下是 本文的亮点: 5.最佳估计和正交投影。Wiener问题是以条件分布与期望的观点处理 的。这样,Wiener理论的基本事实可以迅速获取;结果的范围以及基 本假设可以清楚的显现出来。可以看到所有的统计计算以及结果都基 于一阶和二阶平均;不需要其他的统计数据。这样一来困难(4)便被 排除。这种方法在概率论中为人们所熟知(见Doob[15]第148至155 页以及Lo`eve[16]第455至464页),但在工程上还没有大量的应用。 6.随机过程模型。继前人之后,尤其是Bode与Shannon[3],任意随机 信号可以被表示(直到二阶平均统计性质)为线性动力系统受独立或 不相关随机信号(“白噪声”)激励后的输出。这是工程上应用Wiener 理论的标准手法[2,3,4,5,6,7]。这里用到的方法与传统方法相比只 在线性动力系统的描述方法上不同。我们将强调状态以及状态过渡; 换言之,线性系统将以一阶差分(或微分)方程组来刻画。为了利用 (5)中提到的简化,这种观点是自然的,也是必要的。 7.求解Wiener问题。使用状态——过渡方法,单独一次推导即覆盖多

卡尔曼滤波的基本原理及应用卡尔曼滤波在信号处理与系统控制领域应用广泛,目前,正越来越广泛地应用于计算机应用的各个领域。为了更好地理解卡尔曼滤波的原理与进行滤波算法的设计工作,主要从两方面对卡尔曼滤波进行阐述:基本卡尔曼滤波系统模型、滤波模型的建立以及非线性卡尔曼滤波的线性化。最后,对卡尔曼滤波的应用做了简单介绍。 卡尔曼滤波属于一种软件滤波方法,其基本思想是:以最小均方误差为最佳估计准则,采用信号与噪声的状态空间模型,利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计,求出当前时刻的估计值,算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 最初的卡尔曼滤波算法被称为基本卡尔曼滤波算法,适用于解决随机线性离散系统的状态或参数估计问题。卡尔曼滤波器包括两个主要过程:预估与校正。预估过程主要是利用时间更新方程建立对当前状态的先验估计,及时向前推算当前状态变量和误差协方差估计的值,以便为下一个时间状态构造先验估计值;校正过程负责反馈,利用测量更新方程在预估过程的先验估计值及当前测量变量的基础上建立起对当前状态的改进的后验估计。这样的一个过程,我们称之为预估-校正过程,对应的这种估计算法称为预估-校正算法。以下给出离散卡尔曼滤波的时间更新方程和状态更新方程。 时间更新方程: 状态更新方程: 在上面式中,各量说明如下: A:作用在X k-1上的n×n 状态变换矩阵 B:作用在控制向量U k-1上的n×1 输入控制矩阵 H:m×n 观测模型矩阵,它把真实状态空间映射成观测空间 P k-:为n×n 先验估计误差协方差矩阵 P k:为n×n 后验估计误差协方差矩阵 Q:n×n 过程噪声协方差矩阵 R:m×m 过程噪声协方差矩阵 I:n×n 阶单位矩阵K k:n×m 阶矩阵,称为卡尔曼增益或混合因数 随着卡尔曼滤波理论的发展,一些实用卡尔曼滤波技术被提出来,如自适应滤波,次优滤波以及滤波发散抑制技术等逐渐得到广泛应用。其它的滤波理论也迅速发展,如线性离散系统的分解滤波(信息平方根滤波,序列平方根滤波,UD 分解滤波),鲁棒滤波(H∞波)。 非线性样条自适应滤波:这是一类新的非线性自适应滤波器,它由一个线性组合器后跟挠性无记忆功能的。涉及的自适应处理的非线性函数是基于可在学习

自适应卡尔曼滤波 卡尔曼滤波发散的原因 如果卡尔曼滤波是稳定的,随着滤波的推进,卡尔曼滤波估计的精度应该越来越高,滤波误差方差阵也应趋于稳定值或有界值。但在实际应用中,随着量测值数目的增加,由于估计误差的均值和估计误差协方差可能越来越大,使滤波逐渐失去准确估计的作用,这种现象称为卡尔曼滤波发散。 引起滤波器发散的主要原因有两点: (1)描述系统动力学特性的数学模型和噪声估计模型不准确,不能直接真实地反映物理过程,使得模型与获得的量测值不匹配而导致滤波发散。这种由于模型建立过于粗糙或失真所引起的发散称为滤波发散。 (2)由于卡尔曼滤波是递推过程,随着滤波步数的增加,舍入误差将逐渐积累。如果计算机字长不够长,这种积累误差很有可能使估计误差方差阵失去非负定性甚至失去对称性,使滤波增益矩阵逐渐失去合适的加权作用而导致发散。这种由于计算舍入误差所引起的发散称为计算发散。 针对上述卡尔曼滤波发散的原因,目前已经出现了几种有效抑制滤波发散的方法,常用的有衰减记忆滤波、限定记忆滤波、扩充状态滤波、有限下界滤波、平方根滤波、和自适应滤波等。这些方法本质上都是以牺牲滤波器的最优性为代价来抑制滤波发散,也就是说,多数都是次优滤波方法。 自适应滤波 在很多实际系统中,系统过程噪声方差矩阵Q和量测误差方差阵R事先是不知道的,有时甚至连状态转移矩阵 或量测矩阵H也不能确切建立。如果所建立的模型与实际模型不符可能回引起滤波发散。自适应滤波就是这样一种具有抑制滤波发散作用的滤波方法。在滤波过程中,自适应滤波一方面利用量测值修正预测值,同时也对未知的或不确切的系统模型参数和噪声统计参数进行估计修正。自适应滤波的方法很多,包括贝叶斯法、极大似然法、相关法与协方差匹配法,其中最基本也是最重要的是相关法,而相关法可分为输出相关法和新息相关法。 在这里只讨论系统模型参数已知,而噪声统计参数Q和R未知情况下的自适应滤波。由于Q和R等参数最终是通过增益矩阵K影响滤波值的,因此进行自适应滤波时,也可以不去估计Q和R等参数而直接根据量测数据调整K就可以了。

卡尔曼滤波在单片机上的使用 2011-04-13 13:43:43| 分类:默认分类| 标签:|字号大中小订阅

#ifndef _KALMAN_H_ #define _KALMAN_H_ extern KalmanGain; //卡尔曼增益extern EstimateCovariance; //估计协方差extern MeasureCovariance; //测量协方差extern EstimateValue; //估计值extern void KalmanFilterInit( void ) ; extern KalmanFilter( Measure ); #endif #include "config.h" #include "math.h" KalmanGain; //卡尔曼增益EstimateCovariance; //估计协方差MeasureCovariance; //测量协方差EstimateValue; //估计值 void KalmanFilterInit( void ); extern float KalmanFilter( float Measure ); void KalmanFilterInit( void ) { EstimateValue = 0; EstimateCovariance = 0.1; MeasureCovariance = 0.02; }

KalmanFilter( Measure ) { //计算卡尔曼增益 KalmanGain = EstimateCovariance * sqrt( 1 / ( EstimateCovariance * EstimateCovariance + MeasureCovariance * MeasureCovariance ));//计算本次滤波估计值 EstimateValue = EstimateValue + KalmanGain*(Measure –EstimateValue ); //更新估计协方差 EstimateCovariance = sqrt(1 - KalmanGain) * EstimateCovariance; //更新测量方差 MeasureCovariance = sqrt(1 - KalmanGain) * MeasureCovariance; //返回估计值 return EstimateValue; }

10.6 卡尔曼滤波器简介 本节讨论如何从带噪声的测量数据把有用信号提取出来的问题。通常,信号的频谱处于有限的频率范围内,而噪声的频谱则散布在很广的频率范围内。如前所述,为了消除噪声,可以把 FIR滤波器或IIR滤波器设计成合适的频带滤波器,进行频域滤波。但在许多应用场合,需要进行时域滤波,从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波,但所依据的理论,即针对随机信号的估计理论,是自成体系的。人们对随机信号干扰下的有用信号不能“确知”,只能“估计”。为了“估计”,要事先确定某种准则以评定估计的好坏程度。最小均方误差是一种常用的比较简单的经典准则。典型的线性估计器是离散时间维纳滤波器与卡尔曼滤波器。 对于平稳时间序列的最小均方误差估计的第一个明确解是维纳在1942年2月首先给出的。当时美国的一个战争研究团体发表了一个秘密文件,其中就包括维纳关于滤波问题的研究工作。这项研究是用于防空火力控制系统的。维纳滤波器是基于最小均方误差准则的估计器。为了寻求维纳滤波器的冲激响应,需要求解著名的维纳-霍夫方程。这种滤波理论所追求的是使均方误差最小的系统最佳冲激响应的明确表达式。这与卡尔曼滤波(Kalman filtering)是很不相同的。卡尔曼滤波所追求的则是使均方误差最小的递推算法。 在维纳进行滤波理论研究并导出维纳-霍夫方程的十年以前,在1931年,维纳和霍夫在数学上就已经得到了这个方程的解。 对于维纳-霍夫方程的研究,20世纪五十年代涌现了大量文章,特别是将维纳滤波推广到非平稳过程的文章甚多,但实用结果却很少。这时正处于卡尔曼滤波问世的前夜。 维纳滤波的困难问题,首先在上世纪五十年代中期确定卫星轨道的问题上遇到了。1958年斯韦尔林(Swerling)首先提出了处理这个问题的递推算法,并且立刻被承认和应用。1960年卡尔曼进行了比斯韦尔林更有意义的工作。他严格地把状态变量的概念引入到最小均方误差估计中来,建立了卡尔曼滤波理论。空间时代的到来推动了这种滤波理论的发展。 维纳滤波与卡尔曼滤波所研究的都是基于最小均方误差准则的估计问题。 维纳滤波理论的不足之处是明显的。在运用的过程中,它必须把用到的全部数据存储起来,而且每一时刻都要通过对这些数据的运算才能得到所需要的各种量的估值。按照这种滤波方法设置的专用计算机的存储量与计算量必然很大,很难进行实时处理。虽经许多科技工作者的努力,在解决非平稳过程的滤波问题时,给出能用的方法为数甚少。到五十年代中期,随着空间技术的发展,这种方法越来越不能满足实际应用的需要,面临了新的挑战。尽管如此,维纳滤波理论在滤波理论中的开拓工作是不容置疑的,维纳在方法论上的创见,仍然影响着后人。 五十年代中期,空间技术飞速发展,要求对卫星轨道进行精确的测量。为此,人们将滤波问题以微分方程表示,提出了一系列适应空间技术应用的精练算法。1960年

% 扩展卡尔曼滤波器估计单相电压幅值、相位、频率参数(含直流)function test2_EKF close all; clc; tic; %计时 %模型:y=A0+A1*cos(omega*t+phy1) %离散化:y(k)=A0(k)+A1(k)*cos(omega(k)*k*Ts+phy1(k)) %状态变量:x1(k)=A0(k),x2(k)=omega(k),x3(k)=A1(k)*cos(omega(k)*k*Ts+phy1(k) ),x4(k)=A1(k)*sin(omega(k)*k*Ts+phy1(k)) %下一时刻状态变量为(假设状态不突变):A0(k+1)=A0(k),A1(k+1)=A1(k),omega(k+1)=omega(k),phy1(k+1)=phy1 (k); %则对应状态为:x1(k+1)=x1(k),x2(k+1)=x2(k),x3(k+1)=x3(k)*cos(x2(k)*Ts)- x4(k)*sin(x(2)*Ts),x4(k+1)=x3(k)*sin(x2(k)*Ts)+x4(k)*cos(x(2)*Ts); %状态空间描述:X(k+1)=f(X(k))+W(k);y(k)=H*X(k)+v(k) %f(X(k))=[x1(k);x2(k);x3(k)*cos(x2(k)*Ts)- x4(k)*sin(x(2)*Ts);x3(k)*sin(x2(k)*Ts)+x4(k)*cos(x(2)*Ts)] %偏导(只求了三个):f`(X(k))=[1,0,0;0,1,0;0,-x3(k)*Ts*sin(x2(k)*Ts)-x4(k)*Ts*cos(x2(k)*Ts),cos(x2(k)*Ts);0,x3(k)*Ts*cos(x2(k)*Ts)- x4(k)*Ts*sin(x2(k)*Ts),sin(x2(k)*Ts)]

3. TheKalmanfilterisapowerfulrecursivedatafusiontoolbasedonestimationtheory.Itworksas an estimator forastatevector containingoneor morerandomvariablesthat are assumedtobeGaussiandistributed. Theunderlyingprinciple ofaKalman filterisapredictor-correctorstructure.Inthepredictioncycle,thestatevectordynamicsaredescribedb yastatespacemodel. Thetransition from onestatetothenextoneatatimeisexpressedbyastatetransitionmatrixor,inthenonlinearcase,asanonli nearfunction.Assignedtothestatevectorisacovariancematrix,whichdescribesthestatisticaldispersi onoftherandomstatevectorcomponents.Thecovariancematrixalsorunsthroughthepredictioncycle ,wherethequalityofthestatespacemodelisdescribedbyanadditivenoisecomponent Q . AnimportantpreconditionforusingaKalmanfilteristhatatleastoneoftheelementsinthestatevectoris observableandthatthereisadditionalinformationavailablefortheobservedelement.Thisadditionali nformationcanconsistofe.g.measurementsorknownconstraintsandis combinedinameasurement vector. Theobservation is used forthecorrectioncycle(oralso:filterupdate).Leftmultiplicationofthepredictedstatevectorwiththeo bservationmatrix C yieldsavectorthatissubtractedfromthemeasurementvectortoobtainaresidual.A fterweightingbytheKalmangain,theresidualisaddedtothepredictedstatevectortoobtainafilteresti mate. TheKalmangainincludesthepredictionerrorcovarianceoftheobservedelementsinthestatevector andthemeasurement noise R .From the formulaforthecalculationoftheKalmangainitcan beseen thattheresidualwillbe weightedlessthehigherthemeasurementnoiseis.Inotherwords,iftheaccuracyofthepredictionishig h(lowcovariancevalues)comparedtothemeasurementnoise,theKalmanfilterwill “f ollow ”ratherth estatespacemodelthanthemeasurements.Iftheaccuracyofthepredictionislowbutthequality ofthemeasurementsishigh,theoppositecaseoccurs. AdetailedexplanationoftheKalmanfilteralgorithmanditscompletederivationcanbefoundin[LO90a],[LO90b]. a) LinearKalmanFilter Ifthestatetransitionandtheobservationarelinear,theKalmanfiltercanbedescribedbyasimplesetoff ormulas. Predictioncycle: x ?? =A ?x ?+ ? k +1 =A ?P +?A T +Q with x ?: n-elementstatevector A : n×nstatetransitionmatrix P : covarianceassignedtothestatevector Q : "drivingnoise" Kalmangainandresidual: (B.1.1) K =P ? · C T ??C ?P ? ?1 ?C +R ? k +1 k +1 ? k +1 ? res k +1 =y b (k +1)?C ?x ?k +1 with K : Kalman gain res : residual (B.1.2) y : m-elementmeasurementvector b Updatecycle: C : k +1 k P k T ?

clear all; %% 惯性-GPS组合导航模型参数初始化 we = 360/24/60/60*pi/180; %地球自转角速度,弧度/s psi = 10*pi/180; %psi角度/ 弧度 Tge = 0.12; Tgn = 0.10; Tgz = 0.10; %这三个参数的含义详见参考文献 sigma_ge=1; sigma_gn=1; sigma_gz=1; %% 连续空间系统状态方程 % X_dot(t) = A(t)*X(t) + B(t)*W(t) A=[0 we*sin(psi) -we*cos(psi) 1 0 0 1 0 0; -we*sin(psi) 0 0 0 1 0 0 1 0; we*cos(psi) 0 0 0 0 1 0 0 1; 0 0 0 -1/Tge 0 0 0 0 0; 0 0 0 0 -1/Tgn 0 0 0 0; 0 0 0 0 0 -1/Tgz 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0;]; %状态转移矩阵 B=[0 0 0 sigma_ge*sqrt(2/Tge) 0 0 0 0 0; 0 0 0 0 sigma_gn*sqrt(2/Tgn) 0 0 0 0; 0 0 0 0 0 sigma_gz*sqrt(2/Tgz) 0 0 0;]';%输入控制矩阵%% 转化为离散时间系统状态方程 % X(k+1) = F*X(k) + G*W(k) T = 0.1; [F,G]=c2d(A,B,T);

H=[1 0 0 0 0 0 0 0 0; 0 -sec(psi) 0 0 0 0 0 0 0;];%观测矩阵 %% 卡尔曼滤波器参数初始化 t=0:T:50-T; length=size(t,2); y=zeros(2,length); Q=0.5^2*eye(3); %系统噪声协方差 R=0.25^2*eye(2); %测量噪声协方差 y(1,:)=2*sin(pi*t*0.5); y(2,:)=2*cos(pi*t*0.5); Z=y+sqrt(R)*randn(2,length); %生成的含有噪声的假定观测值,2维X=zeros(9,length); %状态估计值,9维 X(:,1)=[0,0,0,0,0,0,0,0,0]'; %状态估计初始值设定 P=eye(9); %状态估计协方差 %% 卡尔曼滤波算法迭代过程 for n=2:length X(:,n)=F*X(:,n-1); P=F*P*F'+ G*Q*G'; Kg=P*H'/(H*P*H'+R); X(:,n)=X(:,n)+Kg*(Z(:,n)-H*X(:,n)); P=(eye(9,9)-Kg*H)*P; end %% 绘图代码 figure(1) plot(y(1,:)) hold on; plot(y(2,:)) hold off; title('理想的观测量'); figure(2)

1.什么是卡尔曼滤波器 在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”。卡尔曼是一个人的名字。 卡尔曼全名Rudolf Emil Kalman,1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器,正是源于他的博士论文和1960年发表的论文 《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法)。 简单来说,卡尔曼滤波器是一个 “optimal recursive data processing algorithm(最优化自回归数据处理算法)”。对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。 2.卡尔曼滤波器的介绍 为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是,他的5条公式是其核心内容。结合现代的计算机,其实卡尔曼的程序相当的简单,只要你理解了他的那5条公式。 在介绍他的5条公式之前,先让我们来根据下面的例子一步一步的探索。

假设我们要研究的对象是一个房间的温度。根据你的经验判断,这个房间的温度是恒定的,也就是下一分钟的温度等于现在这一分钟的温度(假设我们用一分钟来做时间单位)。假设你对你的经验不是100%的相信,可能会有上下偏差几度。我们把这些偏差看成是高斯白噪声(White Gaussian Noise),也就是这些偏差跟前后时间是没有关系的而且符合高斯分配(Gaussian Distribution)。另外,我们在房间里放一个温度计,但是这个温度计也不准确的,测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。(所谓高斯白噪声中的高斯是指概率分布是正态函数,而白噪声是指它的二阶矩不相关,一阶矩为常数,是指先后信号在时间上的相关性。这是考查一个信号的两个不同方面的问题。 高斯白噪声:如果一个噪声,它的幅度分布服从高斯分布,而它的功率谱密度又是均匀分布的,则称它为高斯白噪声。) 好了,现在对于某一分钟我们有两个有关于该房间的温度值:你根据经验的预测值(系统的预测值)和温度计的值(测量值)。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。 假如我们要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值,来预测k时刻的温度。因为你相信温度是恒定的,所以你会得到k时刻的温度预测值是跟k-1时刻一样的,假设是23度,同时该值的高斯噪声的偏差是5度(5是这样得到的:如果k-1时刻估算出的最优温度值的偏差是3,你对自己预测的不确定度是4度,他们平

Kalman滤波算法的特点: (1)由于Kalman滤波算法将被估计的信号看作在白噪声作用下一个随机线性系统的输出,并且其输入/输出关系是由状态方程和输出方程在时间域内给出的,因此这种滤波方法不仅适用于平稳随机过程的滤波,而且特别适用于非平稳或平稳马尔可夫序列或高斯-马尔可夫序列的滤波,所以其应用范围是十分广泛的。 (2)Kalman滤波算法是一种时间域滤波方法,采用状态空间描述系统。系统的过程噪声和量测噪声并不是需要滤除的对象,它们的统计特征正是估计过程中需要利用的信息,而被估计量和观测量在不同时刻的一、二阶矩却是不必要知道的。 (3)由于Kalman滤波的基本方程是时间域内的递推形式,其计算过程是一个不断地“预测-修正”的过程,在求解时不要求存储大量数据,并且一旦观测到了新的数据,随即可以算的新的滤波值,因此这种滤波方法非常适合于实时处理、计算机实现。 (4)由于滤波器的增益矩阵与观测无关,因此它可预先离线算出,从而可以减少实时在线计算量。在求滤波器增益矩阵时,要求一个矩阵的逆,它的阶数只取决于观测方程的维数,而该维数通常很小,这样,求逆运算是比较方便的。另外,在求解滤波器增益的过程中,随时可以算出滤波器的精度指标P,其对角线上的元素就是滤波误差向量各分量的方差。 Kalman滤波的应用领域 一般地,只要跟时间序列和高斯白噪声有关或者能建立类似的模型的系统,都可以利用Kalman滤波来处理噪声问题,都可以用其来预测、滤波。Kalman滤波主要应用领域有以下几个方面。 (1)导航制导、目标定位和跟踪领域。 (2)通信与信号处理、数字图像处理、语音信号处理。 (3)天气预报、地震预报。 (4)地质勘探、矿物开采。 (5)故障诊断、检测。 (6)证券股票市场预测。 具体事例: (1)Kalman滤波在温度测量中的应用; (2)Kalman滤波在自由落体运动目标跟踪中的应用; (3)Kalman滤波在船舶GPS导航定位系统中的应用; (4)Kalman滤波在石油地震勘探中的应用; (5)Kalman滤波在视频图像目标跟踪中的应用;

卡尔曼滤波器综述 瞿伟军 G10074 1、卡尔曼滤波的起源 1960年,匈牙利数学家卡尔曼发表了一篇关于离散数据线性滤波递推算法的论文,这意味着卡尔曼滤波的诞生。斯坦利.施密特(Stanley Schmidt)首次实现了卡尔曼滤波器,卡尔曼在NASA埃姆斯研究中心访问时,发现他的方法对于解决阿波罗计划的轨道预测很有用,后来阿波罗飞船的导航电脑使用了这种滤波器。关于这种滤波器的论文由Swerling (1958)、Kalman (1960)与 Kalman and Bucy (1961)发表。 2、卡尔曼滤波的发展 卡尔曼滤波是一种有着相当广泛应用的滤波方法,但它既需要假定系统是线性的,又需要认为系统中的各个噪声与状态变量均呈高斯分布,而这两条并不总是确切的假设限制了卡尔曼滤波器在现实生活中的应用。扩展卡尔曼滤波器(EKF)极大地拓宽了卡尔曼滤波的适用范围。EKF的基本思路是,假定卡尔曼滤滤对当前系统状态估计值非常接近于其真实值,于是将非线性函数在当前状态估计值处进行台劳展开并实现线性化。另一种非线性卡尔曼滤波叫线性化卡尔曼滤波。它与EKF的主要区别是前者将非线函数在滤波器对当前系统状态的最优估计值处线性化,而后者因为预先知道非线性系统的实际运行状态大致按照所要求、希望的轨迹变化,所以这些非线性化函数在实际状态处的值可以表达为在希望的轨迹处的台劳展开式,从而完成线性化。 不敏卡尔曼滤波器(UKF)是针对非线性系统的一种改进型卡尔曼滤波器。UKF处理非线性系统的基本思路在于不敏变换,而不敏变换从根本上讲是一种描述高斯随机变量在非线性化变换后的概率分布情况的方法。不敏卡尔曼滤波认为,与其将一个非线性化变换线性化、近似化,还不如将高斯随机变量经非线性变换后的概率分布情况用高斯分布来近似那样简单,因而不敏卡尔曼滤波算法没

卡尔曼滤波在目标跟踪中的应用 摘要:机动卡尔曼算法(VD 算法)在扩展卡尔曼滤波诸算法中原理较为简单,目标跟踪效果也较好。 一. 模型建立 (1) 非机动模型(匀速直线运动) 系统模型 )()()1(k GW k X k X +Φ=+ 其中 ?????? ????? ???=)()()()()(k V k y k V k x k X y x ; ? ? ??????????=Φ10001000010001 T T ; ????? ? ? ???? ???=10200102T T G ? ?? ???=)()()(k W k W k W y x ; 0)]([=k W E ; kj T Q j W k W E δ=)]()([ 测量模型为: )()()(k V k HX k Z +=; 其中 ?? ? ???=01000001H )(k V 为零均值,协方差阵为R 白噪声,与)(k W 不相关。 (2) 机动模型 系统模型 );(*)()1(k W G k X k X m m m m m +Φ=+ 其中

?? ? ? ??? ? ?? ??????????=)()()()()()()(k a k a k V k y k V k x k X m y m y m y m m x m m ;??? ???????????? ?????=Φ100 00 00100000100020100000100200 122 T T T T T T m ;??? ???????????????????=10012040020422T T T T G m 0)]([=k W E m , kj m m m Q j W k W E T δ=)]()([ 观测模型与机动模型的相同,只是H 矩阵为m H 。 ?? ? ???=000100000001m H 二.Kalman 滤波算法 作为一般的Kalman 滤波算法其算法可以描述如下: )1/1(?)1/(?--Φ=-k k X k k X T T G k GQ k k P k k P )1()1/)1()1/(-+Φ--Φ=- 1])1/([)1/()(-+--=R H k k HP H k k P k K T T )]1/()()[()1/(?)/(?--+-=k k HX k Z k K k k X k k X )1/()()1/()/(---=k k HP k K k k P k k P 起始估计值为 ()()()()()()()221/?2/2221/x x x y y y z z z T z z z T ????-??????=????????-???? X 起始估计的估计误差为 (2)(1)(2)(1)2(2/2)(2) (1)(2)(1)2x x x x y y y y v v v T u T v v v T u T -?? ??-?? ?+?? =??-?? -???+???? X 起始估计的估计误差协方差矩阵为

相关文档
最新文档