C++下扩展卡尔曼类(EKF)的实现
程序员文章站
2022-07-12 09:36:30
...
最初的卡尔曼滤波用于解决离散系统的滤波问题,然而工程中常遇到的滤波问题是连续系统产生的滤波问题。
卡尔曼滤波器具有如下形式的离散的状态方程:
但这只是一个高度简化的方程,更多的时候我们能获得的是关于连续系统的如下形式的方程:
其中
针对上述的非线性方程,扩展卡尔曼滤波的做法如下:
首先,状态一定存在一个轨迹
而
因此滤波器的时间更新为
而状态转移矩阵
测量方程随可能同样为非线性
但它可以直接线性化
于是
接下来只剩下测量更新,按照卡尔曼滤波的那一套来执行就行了。
另外,上述在式(3)只保留了一阶项,若嫌精度低可保留二阶项,相应的,状态转移矩阵也要保留二阶项。
一般的,卡尔曼滤波的Q阵要比R阵更难准确获得,因此建议先确定R阵,之后经验性的给定Q阵,观察滤波效果以决定对Q阵的调整。
在C++下将上述过程编写为类。考虑通用性,对所有关于状态
mat equation(time_stamp& time,mat& X,mat& para);
mat是矩阵类,第一个参数time是时刻,第二个参数X是状态,第三个参数para是方程中需要用到的参数。
以四元数微分方程为例:
则
头文件kalman.h,其中矩阵类mat和时间戳类time_stamp请自行实现:
#ifndef KALMAN_H
#define KALMAN_H
//时间更新发生在测量完成之后,即测量完成时刻为当前滤波时刻
//因此时间更新指从上一个滤波时刻更新到当前滤波时刻
typedef mat (*pEquation_of_X)(time_stamp& time,mat& X,mat& para);
class EKF
{
public:
EKF();
~EKF();
//指向状态微分方程
pEquation_of_X pState_differential_equation;
//指向状态微分方程的雅可比矩阵的计算
pEquation_of_X pJacob_of_state_DE;
//指向测量方程的雅可比矩阵即H阵的计算
pEquation_of_X pJacob_of_measure_equation;
//指向测量方程
pEquation_of_X pMeasure_equation;
mat X;//状态
mat P;//状态协方差
mat Z;//测量量
mat Q;//Q阵
mat R;//R阵
mat H;//H阵
mat innovation;//新息
mat gain;//增益阵
mat para_of_time_update;//时间更新的参数矩阵
mat para_of_measure_update;//测量更新的参数矩阵
time_stamp time;//指示滤波器的时间戳
void time_update(void);//时间更新
void measure_update(void);//测量更新
};
#endif
cpp文件kalman.cpp:
#include <iostream>
#include "kalman.h"
EKF::EKF()
{
pState_differential_equation=NULL;
pJacob_of_state_DE=NULL;
pJacob_of_measure_equation=NULL;
pMeasure_equation=NULL;
}
EKF::~EKF()
{}
void EKF::time_update(void)
{
if(pState_differential_equation!=NULL&&pJacob_of_state_DE!=NULL)
{
double detT=time.current_sec-time.last_sec;
X=X+detT*pState_differential_equation(time,X,para_of_time_update);
mat STM=detT*pJacob_of_state_DE(time,X,para_of_measure_update);
for(int i=0;i<STM.rows();i++)
{
STM[i][i]=STM[i][i]+1;
}
P=STM*P*STM.trans()+Q;
}
else
{
std::cout<<"state differential equation is NULL or Jacob matrix of SDE is NULL!"<<std::endl;
}
}
void EKF::measure_update(void)
{
if(pMeasure_equation!=NULL&&pJacob_of_measure_equation!=NULL)
{
innovation=Z-pMeasure_equation(time,X,para_of_measure_update);
H=pJacob_of_measure_equation(time,X,para_of_measure_update);
gain=P*H.trans()*inv(H*P*H.trans()+R);
X=X+gain*innovation;
P=(diag(X.rows(),1)-gain*H)*P;
}
else
{
std::cout<<"measure equation is NULL or jacob of ME is NULL!"<<std::endl;
}
}
上一篇: 纯JS网页上的动态折线图
下一篇: [记录二]学习扩展卡尔曼滤波器