卡尔曼滤波两例题含matlab程序汇总
扩展卡尔曼滤波算法的matlab程序

clear allv=150; %%目标速度v_sensor=0;%%传感器速度t=1; %%扫描周期xradarpositon=0; %%传感器坐标yradarpositon=0; %%ppred=zeros(4,4);Pzz=zeros(2,2);Pxx=zeros(4,2);xpred=zeros(4,1);ypred=zeros(2,1);sumx=0;sumy=0;sumxukf=0;sumyukf=0;sumxekf=0;sumyekf=0; %%%统计的初值L=4;alpha=1;kalpha=0;belta=2;ramda=3-L;azimutherror=0.015; %%方位均方误差rangeerror=100; %%距离均方误差processnoise=1; %%过程噪声均方差tao=[t^3/3 t^2/2 0 0;t^2/2 t 0 0;0 0 t^3/3 t^2/2;0 0 t^2/2 t]; %% the input matrix of process G=[t^2/2 0t 00 t^2/20 t ];a=35*pi/180;a_v=5/100;a_sensor=45*pi/180;x(1)=8000; %%初始位置y(1)=12000;for i=1:200x(i+1)=x(i)+v*cos(a)*t;y(i+1)=y(i)+v*sin(a)*t;endfor i=1:200xradarpositon=0;yradarpositon=0;Zmeasure(1,i)=atan((y(i)-yradarpositon)/(x(i)-xradarpositon))+random('Normal',0,azimutherror,1,1); Zmeasure(2,i)=sqrt((y(i)-yradarpositon)^2+(x(i)-xradarpositon)^2)+random('Normal',0,rangeerror,1,1);xx(i)=Zmeasure(2,i)*cos(Zmeasure(1,i));%%观测值yy(i)=Zmeasure(2,i)*sin(Zmeasure(1,i));measureerror=[azimutherror^2 0;0 rangeerror^2];processerror=tao*processnoise;vNoise = size(processerror,1);wNoise = size(measureerror,1);A=[1 t 0 0;0 1 0 0;0 0 1 t;0 0 0 1];Anoise=size(A,1);for j=1:2*L+1Wm(j)=1/(2*(L+ramda));Wc(j)=1/(2*(L+ramda));endWm(1)=ramda/(L+ramda);Wc(1)=ramda/(L+ramda);%+1-alpha^2+belta; %%%权值if i==1xerror=rangeerror^2*cos(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*sin(Zmeasure(1,i))^2; yerror=rangeerror^2*sin(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*cos(Zmeasure(1,i))^2; xyerror=(rangeerror^2-Zmeasure(2,i)^2*azimutherror^2)*sin(Zmeasure(1,i))*cos(Zmeasure(1,i));P=[xerror xerror/t xyerror xyerror/t;xerror/t 2*xerror/(t^2) xyerror/t 2*xyerror/(t^2);xyerror xyerror/t yerror yerror/t;xyerror/t 2*xyerror/(t^2) yerror/t 2*yerror/(t^2)];xestimate=[Zmeasure(2,i)*cos(Zmeasure(1,i)) 0 Zmeasure(2,i)*sin(Zmeasure(1,i)) 0 ]'; endcho=(chol(P*(L+ramda)))';%for j=1:LxgamaP1(:,j)=xestimate+cho(:,j);xgamaP2(:,j)=xestimate-cho(:,j);endXsigma=[xestimate xgamaP1 xgamaP2];F=A;Xsigmapre=F*Xsigma;xpred=zeros(Anoise,1);for j=1:2*L+1xpred=xpred+Wm(j)*Xsigmapre(:,j);endNoise1=Anoise;ppred=zeros(Noise1,Noise1);for j=1:2*L+1ppred=ppred+Wc(j)*(Xsigmapre(:,j)-xpred)*(Xsigmapre(:,j)-xpred)';endppred=ppred+processerror;chor=(chol((L+ramda)*ppred))';for j=1:LXaugsigmaP1(:,j)=xpred+chor(:,j);XaugsigmaP2(:,j)=xpred-chor(:,j);endXaugsigma=[xpred XaugsigmaP1 XaugsigmaP2 ];for j=1:2*L+1Ysigmapre(1,j)=atan(Xaugsigma(3,j)/Xaugsigma(1,j)) ;Ysigmapre(2,j)=sqrt((Xaugsigma(1,j))^2+(Xaugsigma(3,j))^2);endypred=zeros(2,1);for j=1:2*L+1ypred=ypred+Wm(j)*Ysigmapre(:,j);endPzz=zeros(2,2);for j=1:2*L+1Pzz=Pzz+Wc(j)*(Ysigmapre(:,j)-ypred)*(Ysigmapre(:,j)-ypred)';endPzz=Pzz+measureerror;Pxy=zeros(Anoise,2);for j=1:2*L+1Pxy=Pxy+Wc(j)*(Xaugsigma(:,j)-xpred)*(Ysigmapre(:,j)-ypred)';endK=Pxy*inv(Pzz);xestimate=xpred+K*(Zmeasure(:,i)-ypred);P=ppred-K*Pzz*K';xukf(i)=xestimate(1,1);yukf(i)=xestimate(3,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EKF PRO%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if i==1ekf_p=[xerror xerror/t xyerror xyerror/t;xerror/t 2*xerror/(t^2) xyerror/t 2*xyerror/(t^2);xyerror xyerror/t yerror yerror/t;xyerror/t 2*xyerror/(t^2) yerror/t 2*yerror/(t^2)];ekf_xestimate=[Zmeasure(2,i)*cos(Zmeasure(1,i)) 0 Zmeasure(2,i)*sin(Zmeasure(1,i)) 0 ]';ekf_xpred=ekf_xestimate;end;F=A;ekf_xpred=F*ekf_xestimate;ekf_ppred=F*ekf_p*F'+processerror;H=[-ekf_xpred(3)/(ekf_xpred(3)^2+ekf_xpred(1)^2) 0 ekf_xpred(1)/(ekf_xpred(3)^2+ekf_xpred(1)^2)0;ekf_xpred(1)/sqrt(ekf_xpred(3)^2+ekf_xpred(1)^2) 0 ekf_xpred(3)/sqrt(ekf_xpred(3)^2+ekf_xpred(1)^2) 0];ekf_z(1,1)=atan(ekf_xpred(3)/ekf_xpred(1)) ;ekf_z(2,1)=sqrt((ekf_xpred(1))^2+(ekf_xpred(3))^2);PHHP=H*ekf_ppred*H'+measureerror;ekf_K=ekf_ppred*H'*inv(PHHP);ekf_p=(eye(L)-ekf_K*H)*ekf_ppred;ekf_xestimate=ekf_xpred+ekf_K*(Zmeasure(:,i)-ekf_z);traceekf(i)=trace(ekf_p);xekf(i)=ekf_xestimate(1,1);yekf(i)=ekf_xestimate(3,1);errorx(i)=xx(i)+xradarpositon-x(i);errory(i)=yy(i)+yradarpositon-y(i);ukferrorx(i)=xestimate(1)+xradarpositon-x(i);ukferrory(i)=xestimate(3)+yradarpositon-y(i);ekferrorx(i)=ekf_xestimate(1)+xradarpositon-x(i); ekferrory(i)=ekf_xestimate(3)+yradarpositon-y(i);aa(i)=xx(i)+xradarpositon-x(i);;bb(i)=yy(i)+yradarpositon-y(i);sumx=sumx+(errorx(i)^2);sumy=sumy+(errory(i)^2);sumxukf=sumxukf+(ukferrorx(i)^2);sumyukf=sumyukf+(ukferrory(i)^2);sumxekf=sumxekf+(ekferrorx(i)^2);sumyekf=sumyekf+(ekferrory(i)^2);mseerrorx(i)=sqrt(sumx/(i-1));%噪声的统计均方误差mseerrory(i)=sqrt(sumy/(i-1));mseerrorxukf(i)=sqrt(sumxukf/(i-1));%UKF的统计均方误差mseerroryukf(i)=sqrt(sumyukf/(i-1));mseerrorxekf(i)=sqrt(sumxekf/(i-1));%EKF的统计均方误差mseerroryekf(i)=sqrt(sumyekf/(i-1));endfigure(1);plot(mseerrorxukf,'r');hold on;plot(mseerrorxekf,'g');hold on;plot(mseerrorx,'.');hold on;ylabel('MSE of X axis','fontsize',15);xlabel('sample number','fontsize',15);legend('UKF','EKF','measurement error');figure(2)plot(mseerroryukf,'r');hold on;plot(mseerroryekf,'g');hold on;plot(mseerrory,'.');hold on;ylabel('MSE of Y axis','fontsize',15); xlabel('sample number','fontsize',15); legend('UKF','EKF','measurement error');figure(3)plot(x,y);hold on;plot(xekf,yekf,'g');hold on;plot(xukf,yukf,'r');hold on;plot(xx,yy,'m');ylabel(' X ','fontsize',15);xlabel('Y','fontsize',15);legend('TRUE','UKF','EKF','measurements');。
扩展卡尔曼滤波(EKF)理论讲解与实例(matlab、python和C++代码)

扩展卡尔曼滤波(EKF)理论讲解与实例(matlab、python和C++代码)扩展卡尔曼滤波(EKF)理论讲解与实例(matlab、python和C++代码)⽂章⽬录我们上篇提到的 (参见我的另⼀篇⽂章: )是⽤于线性系统,预测(运动)模型和观测模型是在假设⾼斯和线性情况下进⾏的。
简单的卡尔曼滤波必须应⽤在符合⾼斯分布的系统中,但是现实中并不是所有的系统都符合这样 。
另外⾼斯分布在⾮线性系统中的传递结果将不再是⾼斯分布。
那如何解决这个问题呢?扩展卡尔曼滤波就是⼲这个事的。
理论讲解扩展卡尔曼滤波(Extended Kalman Filter,EKF)通过局部线性来解决⾮线性的问题。
将⾮线性的预测⽅程和观测⽅程进⾏求导,以切线代替的⽅式来线性化。
其实就是在均值处进⾏⼀阶泰勒展开。
数学中,泰勒公式是⼀个⽤函数在某点的信息描述其附近取值的公式( ⼀句话描述:就是⽤多项式函数去逼近光滑函数 )。
如果函数⾜够平滑的话,在已知函数在某⼀点的各阶导数值的情况之下,泰勒公式可以⽤这些导数值做系数构建⼀个多项式来近似函数在这⼀点的邻域中的值。
泰勒公式还给出了这个多项式和实际的函数值之间的偏差。
表⽰ 在第 阶导数的表达式,带⼊⼀个值计算后得到的结果(注意,它是个值)是⼀个系数(⼀个值),每⼀项都不同,第⼀项 ,第⼆项 …… 依此类推是⼀个以为⾃变量的表达式 。
是泰勒公式的余项,是 的⾼阶⽆穷⼩KF 和EKF 模型对⽐⾸先,让卡尔曼先和扩展卡尔曼滤波做⼀个对⽐。
在对⽐过程中可以看出,扩展卡尔曼是⼀个简单的⾮线性近似滤波算法,指运动或观测⽅程不是线性的情况,在预测模型部分,扩展卡尔曼的预测模型和量测模型已经是⾮线性了。
为了简化计算,EKF 通过⼀阶泰勒分解线性化运动、观测⽅程。
KF 与EKF 具有相同的算法结构,都是以⾼斯形式描述后验概率密度的,通过计算贝叶斯递推公式得到的。
最⼤的不同之处在于,计算⽅差时,EKF 的状态转移矩阵(上⼀时刻的状态信息)和观测矩阵(⼀步预测)都是状态信息的雅克⽐矩阵( 偏导数组成的矩阵)。
自适应扩展卡尔曼滤波matlab

自适应扩展卡尔曼滤波matlab自适应扩展卡尔曼滤波(Adaptive Extended Kalman Filter,AEKF)是一种用于非线性系统状态估计的滤波算法。
本文将介绍AEKF算法的原理、步骤和实现方法,并结合MATLAB 编写代码进行演示。
一、扩展卡尔曼滤波原理扩展卡尔曼滤波(Extended Kalman Filter,EKF)是一种用于非线性系统状态估计的滤波算法。
它通过使用线性化系统模型的方式将非线性系统转换为线性系统,在每个时间步骤中用线性卡尔曼滤波器进行状态估计。
然而,EKF仅限于具有凸多边形测量特性的问题,并且对线性化过程误差敏感。
为了解决这些问题,AEKF通过自适应更新协方差矩阵的方式提高了滤波器的性能。
AEKF通过测量残差的方差更新协方差矩阵,从而提高了滤波器对非线性系统的适应能力。
AEKF算法的步骤如下:1. 初始化状态向量和协方差矩阵。
2. 根据系统的非线性动力学方程和测量方程计算预测状态向量和协方差矩阵。
3. 计算测量残差,即测量值与预测值之间的差值。
4. 计算测量残差的方差。
5. 判断测量残差的方差是否超过预设阈值,如果超过,则更新协方差矩阵。
6. 利用更新后的协方差矩阵计算最优滤波增益。
7. 更新状态向量和协方差矩阵。
8. 返回第2步,进行下一次预测。
二、AEKF算法的MATLAB实现下面,我们将使用MATLAB编写AEKF算法的代码,并通过一个实例进行演示。
首先,定义非线性系统的动力学方程和测量方程。
在本例中,我们使用一个双摆系统作为非线性系统模型。
```matlabfunction x_next = nonlinear_dynamics(x_current, u)% Nonlinear system dynamicstheta1 = x_current(1);theta2 = x_current(2);d_theta1 = x_current(3);d_theta2 = x_current(4);g = 9.8; % Gravitational accelerationd_theta1_next = d_theta1 + dt * (-3*g*sin(theta1) - sin(theta1-theta2) ...+ 2*sin(theta1-theta2)*(d_theta2^2 + d_theta1^2*cos(theta1-theta2))) .../ (3 - cos(2*(theta1-theta2)));d_theta2_next = d_theta2 + dt * (2*sin(theta1-theta2)*(2*d_theta2^2 ...+ d_theta1^2*cos(theta1-theta2) + g*cos(theta1) +g*cos(theta1-theta2))) .../ (3 - cos(2*(theta1-theta2)));theta1_next = theta1 + dt * d_theta1_next;theta2_next = theta2 + dt * d_theta2_next;x_next = [theta1_next; theta2_next; d_theta1_next;d_theta2_next];endfunction y = measurement_model(x)% Measurement model, measure the angles of the double pendulumtheta1 = x(1);theta2 = x(2);y = [theta1; theta2];end```然后,定义AEKF算法的实现。
matlab扩展卡尔曼滤波实例

matlab扩展卡尔曼滤波实例在MATLAB中实现卡尔曼滤波的过程。
第一步:准备数据要使用卡尔曼滤波算法,首先需要准备尽可能准确的测量数据。
我们假设我们正在跟踪一个匀速运动的物体,我们使用一个简单的模型来生成数据。
以下是一个MATLAB代码示例:matlab生成匀速运动数据t = 0:0.1:10; 时间范围v = 2; 速度x_true = 5 + v*t; 真实位置x_meas = x_true + 0.5*randn(size(t)); 测量位置(加入噪声)上述代码生成了10秒钟内物体的真实位置(x_true)和加入高斯噪声的测量位置(x_meas)。
第二步:初始化卡尔曼滤波器在开始使用卡尔曼滤波器之前,需要初始化滤波器的状态估计和协方差矩阵。
以下是一个示例代码:matlab初始化卡尔曼滤波器参数x_est = 0; 状态估计P_est = 1; 协方差矩阵Q = 1; 过程噪声方差R = 0.5; 测量噪声方差在上述代码中,我们初始化了状态估计(x_est)、协方差矩阵(P_est)、过程噪声方差(Q)和测量噪声方差(R)。
第三步:卡尔曼滤波迭代卡尔曼滤波的核心是迭代过程,其中包含两个关键步骤:预测和更新。
预测步骤是使用系统模型来预测下一时刻的状态和协方差矩阵。
更新步骤将测量值与预测结果进行比较,以改进状态估计和协方差矩阵。
以下是一个MATLAB代码示例:matlab卡尔曼滤波迭代for k = 1:length(t)预测步骤x_pred = x_est;P_pred = P_est + Q;更新步骤K = P_pred / (P_pred + R);x_est = x_pred + K*(x_meas(k) - x_pred);P_est = (1 - K)*P_pred;保存状态估计结果x_est_hist(k) = x_est;end在上述代码中,我们首先进行预测步骤,计算预测状态(x_pred)和预测协方差矩阵(P_pred)。
卡尔曼滤波 matlab算法

卡尔曼滤波 matlab算法卡尔曼滤波是一种用于状态估计的数学方法,它通过将系统的动态模型和测量数据进行融合,可以有效地估计出系统的状态。
在Matlab中,实现卡尔曼滤波算法可以通过以下步骤进行:1. 确定系统的动态模型,首先需要将系统的动态模型表示为状态空间方程,包括状态转移矩阵、控制输入矩阵和过程噪声的协方差矩阵。
2. 初始化卡尔曼滤波器,在Matlab中,可以使用“kf = kalmanfilter(StateTransitionModel, MeasurementModel, ProcessNoise, MeasurementNoise, InitialState, 'State', InitialCovariance)”来初始化一个卡尔曼滤波器对象。
其中StateTransitionModel和MeasurementModel分别是状态转移模型和测量模型,ProcessNoise和MeasurementNoise是过程噪声和测量噪声的协方差矩阵,InitialState是初始状态向量,InitialCovariance是初始状态协方差矩阵。
3. 进行预测和更新,在每个时间步,通过调用“predict”和“correct”方法,可以对状态进行预测和更新,得到最优的状态估计值。
4. 实时应用,将测量数据输入到卡尔曼滤波器中,实时获取系统的状态估计值。
需要注意的是,在实际应用中,还需要考虑卡尔曼滤波器的参数调节、性能评估以及对不确定性的处理等问题。
此外,Matlab提供了丰富的工具箱和函数,可以帮助用户更便捷地实现和应用卡尔曼滤波算法。
总的来说,实现卡尔曼滤波算法需要对系统建模和Matlab编程有一定的了解和技能。
希望以上内容能够对你有所帮助。
卡尔曼滤波计算举例全

卡尔曼滤波计算举例⏹计算举例⏹卡尔曼滤波器特性假设有一个标量系统,信号与观测模型为[1][][]x k ax k n k +=+[][][]z k x k w k =+其中a 为常数,n [k ]和w [k ]是不相关的零均值白噪声,方差分别为和。
系统的起始变量x [0]为随机变量,其均值为零,方差为。
2nσ2σ[0]x P (1)求估计x [k ]的卡尔曼滤波算法;(2)当时的卡尔曼滤波增益和滤波误差方差。
220.9,1,10,[0]10nx a P =σ=σ==1. 计算举例根据卡尔曼算法,预测方程为:ˆˆ[/1][1/1]xk k ax k k -=--预测误差方差为:22[/1][1/1]x x nP k k a P k k -=--+σ卡尔曼增益为:()1222222[][/1][/1][1/1][1/1]x x x nx n K k P k k P k k a P k k a P k k -=--+σ--+σ=--+σ+σˆˆˆ[/][/1][]([][/1])ˆˆ[1/1][]([][1/1])ˆ(1[])[1/1][][]xk k x k k K k z k x k k axk k K k z k ax k k a K k xk k K k z k =-+--=--+---=---+滤波方程:()()2222222222222[/](1[])[/1][1/1]1[1/1][1/1][1/1][1/1]x x x nx n x n x nx nP k k K k P k k a P k k a P k k a P k k a P k k a P k k =--⎛⎫--+σ=---+σ ⎪--+σ+σ⎝⎭σ--+σ=--+σ+σ滤波误差方差起始:ˆ[0/0]0x=[0/0][0]x x P P =k [/1]x P k k -[/]x P k k []K k 012345689104.76443.27012.67342.27652.21422.18362.16832.16089.104.85923.64883.16542.94752.84402.79352.76870.47360.32700.26730.24040.22770.22140.21840.2168ˆ[0/0]0x=[0/0]10x P =220.9110na =σ=σ=2. 卡尔曼滤波器的特性从以上计算公式和计算结果可以看出卡尔曼滤波器的一些特性:(1)滤波误差方差的上限取决于测量噪声的方差,即()2222222[1/1][/][1/1]x nx x na P k k P k k a P k k σ--+σ=≤σ--+σ+σ2[/]x P k k ≤σ这是因为(2)预测误差方差总是大于等于扰动噪声的方差,即2[/1]x nP k k -≥σ这是因为222[/1][1/1]x x n nP k k a P k k -=--+σ≥σ(3)卡尔曼增益满足,随着k 的增加趋于一个稳定值。
容积卡尔曼滤波 matlab

容积卡尔曼滤波matlab摘要:1.容积卡尔曼滤波简介2.容积卡尔曼滤波算法原理3.容积卡尔曼滤波算法在MATLAB 中的实现4.容积卡尔曼滤波算法的应用案例5.结论正文:一、容积卡尔曼滤波简介容积卡尔曼滤波(Cubature Kalman Filter,简称CKF)是一种基于卡尔曼滤波理论的非线性滤波算法。
它通过将非线性系统的状态空间模型和观测模型进行离散化,采用立方插值方法对系统状态进行预测和更新,从而实现对非线性系统的状态估计。
相较于传统的卡尔曼滤波,容积卡尔曼滤波具有更好的性能和鲁棒性,被广泛应用于导航定位、目标跟踪、机器人控制等领域。
二、容积卡尔曼滤波算法原理容积卡尔曼滤波算法主要包括两个部分:预测阶段和更新阶段。
1.预测阶段在预测阶段,首先对系统的状态向量进行初始化,然后通过系统动态模型和观测模型,对系统的状态向量进行预测。
具体来说,根据系统的状态转移矩阵、控制矩阵、观测矩阵和过程噪声协方差矩阵,计算预测状态向量的均值和协方差矩阵。
2.更新阶段在更新阶段,根据预测的观测值和观测协方差矩阵,计算观测均值和协方差矩阵。
然后,利用卡尔曼增益公式,结合预测状态向量和观测均值,更新系统的状态向量和协方差矩阵。
三、容积卡尔曼滤波算法在MATLAB 中的实现在MATLAB 中,可以通过以下步骤实现容积卡尔曼滤波算法:1.导入所需库:`import numpy as np;`2.初始化状态向量和协方差矩阵:`x = np.zeros((2,1)); p =np.zeros((2,2));`3.设置系统参数:`F = np.array([[1, 0.1], [0, 1]]); Q = np.array([[0.1, 0], [0, 0.1]]); H = np.array([[1, 0], [0, 1]]);`4.预测阶段:`F_pred = F; Q_pred = Q; x_pred = F_pred @ x; S_pred = F_pred @ P @ F_pred.T + Q_pred;`5.更新阶段:`y=H@x;S_update=H@*****+R;`6.计算卡尔曼增益:`K=*****@np.linalg.inv(S_update);`7.更新状态向量和协方差矩阵:`x = x + K @ (y - H @ x); P = (np.eye(2) - K @ H) @ P;`四、容积卡尔曼滤波算法的应用案例容积卡尔曼滤波算法在各种领域都有广泛应用,例如:1.导航定位:利用GPS、惯性导航等传感器的数据,实现对飞行器、船舶等移动设备的精确定位。
卡尔曼滤波器例题

卡尔曼滤波器是一种递归滤波器,它使用测量数据和系统状态方程来估计一个动态系统的状态。
下面是一个简单的卡尔曼滤波器例题:假设我们有一个一维系统,其状态方程为:X(k) = X(k-1) + 1其中X(k) 表示系统在时刻k 的状态,X(k-1) 表示系统在时刻k-1 的状态。
测量方程为:Z(k) = X(k) + 2 * R * randn()其中Z(k) 表示在时刻k 的测量值,R 表示测量噪声的方差,randn() 表示随机数生成函数,用于生成符合标准正态分布的随机噪声。
现在,我们使用卡尔曼滤波器来估计系统在时刻k 的状态X(k)。
在卡尔曼滤波器中,我们需要定义以下几个变量:P(k|k-1):表示在时刻k-1 时对时刻k 的状态估计的不确定性。
K(k):表示卡尔曼增益,用于将测量值和估计值结合起来。
P(k):表示在时刻k 时对状态估计的不确定性。
X(k):表示在时刻k 的状态估计值。
在初始状态下,我们设定X(0|0)=0 和P(0|0)=1。
然后,我们使用以下递归公式来更新这些变量:P(k|k-1) = P(k-1|k-1) + Q,其中Q 表示过程噪声的方差。
K(k) = P(k|k-1) / (P(k|k-1) + R),其中R 表示测量噪声的方差。
X(k) = X(k|k-1) + K(k) * (Z(k) - X(k|k-1)),其中Z(k) 是测量值。
P(k) = (1 - K(k)) * P(k|k-1)。
下面是一个使用Python 实现的卡尔曼滤波器示例代码:import numpy as npfrom scipy.linalg import inv# 定义系统参数和测量参数Q = 0.01 # 过程噪声方差R = 0.25 # 测量噪声方差X_true = np.arange(1, 10) # 真实状态序列Z_measured = X_true + 2 * R * np.random.randn(len(X_true)) # 测量序列# 初始化卡尔曼滤波器参数X_est = np.zeros_like(X_true) # 初始状态估计值P_est = np.ones_like(X_true) # 初始状态估计的不确定性X_est[0] = X_true[0] # 初始状态估计值为真实值P_est[0] = 1 # 初始不确定性设为1K = np.zeros_like(X_true) # 卡尔曼增益数组for k in range(1, len(X_true)):# 预测步骤:根据上一时刻的状态和过程方程预测当前时刻的状态和不确定性P_est[k] = P_est[k-1] + Q # P预测= P上一时刻+ QX_est[k] = X_est[k-1] + 1 # X预测= X上一时刻+ 1# 更新步骤:根据当前时刻的测量值和卡尔曼增益更新状态估计和不确定性K[k] = P_est[k] / (P_est[k] + R) # 卡尔曼增益计算公式X_est[k] = X_est[k] + K[k] * (Z_measured[k] - X_est[k]) # 状态估计更新公式P_est[k] = (1 - K[k]) * P_est[k] # 不确定性更新公式。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
设高度的测量误差是均值为0、方差为1的高斯白噪声随机序列,该物体的初始高度0h 和速度0V 也是高斯分布的随机变量,且0000019001000,var 10/02Eh h m P EV m s V ⎡⎤⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦。
试求该物体高度和速度随时间变化的最优估计。
(2/80.9s m g =) 解:1. 令()()()h k X k v k ⎡⎤=⎢⎥⎣⎦t=1 R (k )=1 Q(k)=0 根据离散时间卡尔曼滤波公式,则有: (1)(1,)()()X k k k X k U k φ+=++ (1)(1)(1)(1)Y k H k X k V k +=++++(1,)k k φ+= 11t -⎡⎤⎢⎥⎣⎦ ()U k = 20.5gt gt ⎡⎤-⎢⎥⎣⎦(1)H k +=[]10 滤波初值:^1900(0|0)(0)10X EX ⎡⎤==⎢⎥⎣⎦0100(0|0)var[(0)]2P X P ⎡⎤===⎢⎥⎣⎦一步预测:^^(1|)(1,)(|)()X k k k k X k k U k φ+=++ (1|)(1,)(|)(1,)TP k k k k P k k k k φφ+=++滤波增益:1(1)(1|)(1)[(1)(1|)(1)(1)]TTK k P k k H k H k P k k H k R k -+=+++++++ 滤波计算:^^^(1|1)(1|)(1)[(1)(1)(1|)]X k k X k k K k Y k H k X k k ++=++++-++ (1|1)[(1)(1)](1|)P k k I K k H k P k k ++=-+++ 2. 实验结果高度随时间变化估计速度随时间变化的最优估计高度协方差速度协方差从以上的结果,可以得到高度和速度的估计值,再通过所得到的高度协方差和速度协方差,可见用卡尔曼滤波法,虽然刚开始的初始高度协方差很大为100,但通过2步之后减小到不超过1,逐渐接近于0, 同样的速度协方差刚开始的时候也比较大,为2,但是通过5步之后迅速减小,到10步之后接近于0。
3. 有关参数的影响(例如初始条件、噪声统计特性对滤波结果的影响等);1)初始条件改变时,改变初始高度值,和速度值 00230030/Eh m EV m s ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦由实验结果分析可得度滤波值和速度滤波值在开始几步接近初始值,协方差值基本不变。
2)当初始协方差值改变时,改为0001500var 010h P V ⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎣⎦实验结果分析高度和速度滤波值基本不变,速度协方差和高度协方差开始要接近速度协方差和高度协方差的初始值。
但是经过几步之后,都趋于0。
二.同样考虑自由落体运动的物体,用雷达(和物体落地点在同一水平面)进行测量,如图所示。
如果⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡200050005var ,/1200519950000000V h d P s m m m EV Eh Ed ,且雷达测距和测角的测量噪声是高斯白噪声随机序列,均值为零、方差阵⎥⎦⎤⎢⎣⎡=01.00004.0R ,试根据下列测量数据确定物体的高度和速度随时间变化的估计值。
时间[s]*1000 斜距[km] 俯仰角[rad]*1000 0.00050000000000 2.82741643781891 0.00075850435876 0.00100000000000 2.82519811729771 0.00083282260478 0.00150000000000 2.82066686966236 0.000678082416390.00200000000000 2.81487233105901 0.00085279036802 0.00250000000000 2.80671786536244 0.00072900768452 0.00300000000000 2.79725268974089 0.00080072481819 0.00350000000000 2.78664273475039 0.00075095576213 0.00400000000000 2.77320365026313 0.00065762725379 0.00450000000000 2.75919535464551 0.00081186148545 0.00500000000000 2.74331288628195 0.00079783727034 0.00550000000000 2.72538888482812 0.00073060712986 0.00600000000000 2.70664967712312 0.00063242006530 0.00650000000000 2.68632403406473 0.00063656524495 0.00700000000000 2.66386533852220 0.00080659845639 0.00750000000000 2.64093529707333 0.00067704740069 0.00800000000000 2.61621111727357 0.00076573767706 0.00850000000000 2.59038109850785 0.00054955759081 0.00900000000000 2.56298794272843 0.00058487913971 0.00950000000000 2.53498317950797 0.00055602747368 0.01000000000000 2.50647589372246 0.00033550412588 0.01050000000000 2.47571075016386 0.00056012688452 0.01100000000000 2.44560676000982 0.00056694491978 0.01150000000000 2.41403690772088 0.00059380631025 0.01200000000000 2.38252228611696 0.00053681916544 0.01250000000000 2.35016501182332 0.00065871960781 0.01300000000000 2.31790939837137 0.00068598344328 0.01350000000000 2.28597616656453 0.00060922471348 0.01400000000000 2.25418431681401 0.00057086018918 0.01450000000000 2.22259320219535 0.00041308535708 0.01500000000000 2.19237398969466 0.00047302026281 0.01550000000000 2.16290177997271 0.00030949309972 0.01600000000000 2.13441725793706 0.00040552624986 0.01650000000000 2.10811064690727 0.00037545033142 0.01700000000000 2.08322179823195 0.00017282319262 0.01750000000000 2.06148109026767 0.00020758327980 0.01800000000000 2.04219885094031 0.00037186464579 0.01850000000000 2.02610235314357 0.00018082163465 0.01900000000000 2.01290326863579 0.00023323830160 0.01950000000000 2.00463157388395 -0.00004536186964 0.02000000000000 2.00058143251913 0.00003246284068解: 1.令()()()d X k h k v k ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦t=0.5 0.04()0.01R k ⎡⎤=⎢⎥⎣⎦ Q(k)=0 ()()()l k Y k k θ⎡⎤=⎢⎥⎣⎦根据离散时间扩展卡尔曼滤波公式,则有:状态方程:(1)(1,)()()()X k k k X k U k W k φ+=+++ 20()0.5U k gt gt ⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦(1,)k k φ+=111t ⎡⎤⎢⎥-⎢⎥⎢⎥⎣⎦测量方程:(1)()()arctan Y k V k h k d +=+⎢⎥⎢⎥⎣⎦辅助方程:^(1)(1|)[(1),1](1)|(1)x k x k k h X k k H k X k +=+∂+++=∂+^22^^22^22^^22(1|)0(1|)(1|)0(1|)(1|)(1|)h k k dh k k d h k k d d h k k h k k d h k k d ⎡⎤+⎢⎥⎢⎥++++⎢⎥=⎢⎥⎢⎥-+⎢⎥⎢⎥++++⎣⎦一步预测:^^(1|)(1,)(|)()X k k k k X k k U k φ+=++ (1|)(1,)(|)(1,)TP k k k k P k k k k φφ+=++滤波增益:1(1)(1|)(1)[(1)(1|)(1)(1)]TTK k P k k H k H k P k k H k R k -+=+++++++ 滤波计算:^^^(1|1)(1|)(1)[(1)((1|),1)]X k k X k k K k Y k h X k k k ++=++++-++(1|1)[(1)(1)](1|)P k k I K k H k P k k ++=-+++滤波初值: ^1995(0|0)20051X ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦ 5(0|0)52P ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦2. 实验结果高度随时间变化估计速度随时间变化的估计高度协方差:速度协方差估计实验结果分析:根据图,可得高度和速度的估计值,通过扩展卡尔曼滤波法,高度协方差和速度协方差,刚开始的值比较大,但是迅速减小,在几步之后逐渐趋近于0。
3. 有关参数的影响(例如初始条件、噪声统计特性对滤波结果的影响等);1)初始条件发生变化,改变高度和速度的初始值为^1995 (0|0)230010 X⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦分析:高度和速度在刚开始的时候波动比较大,不过经过10步之后,逐渐趋于平稳,高度协方差收敛变快,速度协方差基本不变2)初始噪声改变,1 ()1 R k⎡⎤=⎢⎥⎣⎦分析:高度和速度滤波值基本不变,速度协方差也基本不变,高度协方差刚开始的时候有波动,10步之后趋于稳定。