输出倒逼思考

0%

从贝叶斯滤波到扩展卡尔曼滤波

发表于
分类于 ， Valine：
扩展卡尔曼滤波是标准卡尔曼滤波在非线性情形下的一种扩展形式，它是一种高效率的递归滤波器（自回归滤波器）。扩展卡尔曼滤波的基本思想是利用泰勒级数展开将非线性系统的状态转移函数 f(x) 和（或）观测函数 h(x) 线性化，然后采用卡尔曼滤波框架对信号进行滤波，因此它是一种次优滤波。

目录

目录-从贝叶斯滤波到扩展卡尔曼滤波

0 前言

扩展卡尔曼滤波（Extended Kalman Filter，EKF）是标准卡尔曼滤波在非线性情形下的一种扩展形式，它是一种高效率的递归滤波器（自回归滤波器）。

EKF 的基本思想是利用泰勒级数展开将非线性系统的状态转移函数 $f(x)$ 和（或）观测函数 $h(x)$ 线性化，然后采用卡尔曼滤波框架对信号进行滤波，因此它是一种次优滤波。

1 贝叶斯滤波的三大概率密度函数

在此前的文章《从概率到贝叶斯滤波》和《从贝叶斯滤波到卡尔曼滤波》中，已经讲到贝叶斯滤波的先验概率密度函数、似然概率密度函数和后验概率密度函数：

(1) 先验概率密度函数

(2) 似然概率密度函数

(3) 后验概率密度函数

其中，后验概率密度函数中的归一化常数 $\eta_k$ 为：

2 扩展卡尔曼滤波的假设

扩展卡尔曼滤波以贝叶斯滤波为理论基础，并作了五个前提假设。

2.1 与卡尔曼滤波相同的假设

(1) 假设一：状态量服从正态分布

(2) 假设二：观测量服从正态分布

(3) 假设三：过程噪声服从均值为 0 的正态分布

(4) 假设四：观测噪声服从均值为 0 的正态分布

2.2 与卡尔曼滤波不同的假设

(5) 假设五：状态转移函数和（或）观测函数为非线性函数

在卡尔曼滤波的前提假设中，认为状态方程中的状态转移函数 $f(x)$ 以及观测方程中的观测函数 $h(x)$ 均为线性函数。基于这种线性假设，存在常数或常矩阵 $F$，使得 $f(x)$ 可以写成卡尔曼滤波中的线性形式，存在常数或常矩阵 $H$，使得 $h(x)$ 也可以写成卡尔曼滤波中的线性形式。

不同于标准卡尔曼滤波，扩展卡尔曼滤波处理的是非线性系统，假设系统的状态转移函数和（或）观测函数为非线性函数。

3 扩展卡尔曼滤波的公式推导

3.1 预测步的两个公式

当状态转移函数为线性函数时，扩展卡尔曼滤波的预测步与标准卡尔曼滤波相同；当状态转移函数为非线性函数时，扩展卡尔曼滤波的预测步采用下面的推导过程。

根据假设一，$k-1$ 时刻状态量 $X_{k-1}$ 服从均值为 $\mu_{k-1}^+$，方差为 ${\sigma_{k-1}^+}^2$ 的正态分布：

$X_{k-1}$ 的后验概率密度函数为：

对状态转移函数 $f(x)$ 在 $X_{k-1}$ 的后验估计处 $\hat{x}_{k-1}^+$（即 $\mu_{k-1}^+$）进行一阶泰勒级数展开：

其中，$o(X_{k-1}-\hat{x}_{k-1}^+)$ 为高阶无穷小量，对其进行舍弃，便可得到状态转移函数 $f(x)$ 的近似形式：

令 $F_j=f’(\hat{x}_{k-1}^+), \quad C_F=f(\hat{x}_{k-1}^+)-f’(\hat{x}_{k-1}^+)\hat{x}_{k-1}^+$，则近似有：

结合假设三与上式，$k$ 时刻状态量 $X_k$ 的先验概率密度函数为：

类似推导标准卡尔曼滤波时的做法，使用 Mathematica 软件做符号推导，整理表达式结果可知，先验概率密度函数 $f_{X_k}^-(x)$ 为正态分布函数，均值和方差分别为：

其中，$F_j=f’(\hat{x}_{k-1}^+)$，对于一维扩展卡尔曼滤波，$F_j$ 为一常数项，即状态转移函数 $f(x)$ 的一阶导函数 $f’(x)$ 在上一时刻状态量后验估计处的函数值。

3.2 更新步的三个公式

当观测函数为线性函数时，扩展卡尔曼滤波的更新步与标准卡尔曼滤波相同；当观测函数为非线性函数时，扩展卡尔曼滤波的更新步采用下面的推导过程。

对观测函数 $h(x)$ 在 $X_k$ 的先验估计处 $\hat{x}_k^-$（即 $\mu_k^-$）进行一阶泰勒级数展开：

其中，$o(X_k-\hat{x}_k^-)$ 为高阶无穷小量，对其进行舍弃，便可得到观测函数 $f(x)$ 的近似形式：

令 $H_j=h’(\hat{x}_k^-), \quad C_H=h(\hat{x}_k^-)-h’(\hat{x}_k^-)\hat{x}_k^-$，则近似有：

结合假设四、公式 (3.1)、公式 (3.2) 及 $h(x)$ 的线性化近似形式，可知，k 时刻状态量 $X_k$ 的后验概率密度函数为：

其中，归一化常数 $\eta_k$ 为：

类似推导标准卡尔曼滤波时的做法，使用 Mathematica 软件做符号推导，整理表达式结果可知，后验概率密度函数 $f_{X_k}^+(x)$ 为正态分布函数，均值和方差分别为：

$\mu_k^+$ 即 k 时刻状态量 $X_k$ 的后验估计 $\hat{x}_k^+$。其中，$K$ 被称为卡尔曼增益系数：

其中，$H_j=h’(\hat{x}_k^-)$，对于一维扩展卡尔曼滤波，$H_j$ 为一常数项，即观测函数 $h(x)$ 的一阶导函数 $h’(x)$ 在当前时刻状态量先验估计处的函数值。

4 矩阵形式的扩展卡尔曼滤波

上文内容所描述的是一维的扩展卡尔曼滤波，当状态量和观测量不再是单一的随机变量而是由多个随机变量组成的序列时，扩展卡尔曼滤波中各个量的维数也将随之改变：

  • 状态量 $X$ 由随机变量演变为随机向量，随机向量中的每一个分量为一个状态量随机变量。维数为 $n_X \times 1$
  • 状态转移比例项 $F_j$ 演变为雅可比矩阵，维数为 $n_X \times n_X$。此时，$F_j$ 是由每个状态转移函数对每个状态量分量求偏导后，将上一时刻状态量的后验估计代入得到：
  • 状态量概率密度函数均值 $\mu$ 演变为矩阵，维数为 $n_X \times 1$
  • 状态量概率密度函数方差 $\sigma^2$ 演变为协方差矩阵，用 $\Sigma$ 表示，维数为 $n_X \times n_X$
  • 过程噪声方差 ${\sigma_Q}^2$ 演变为协方差矩阵，用 $\Sigma_Q$ 表示，维数为 $n_X \times n_X$
  • 观测量 $Y$ 由随机变量演变为随机向量，随机向量中的每一个分量为一个观测量随机变量。维数为 $n_Y \times 1$
  • 观测值 $y_k$ 由单一值演变为由单一值组成的值矩阵，维数为 $n_Y \times 1$
  • 观测比例项 $H_j$ 演变为雅可比矩阵，维数为 $n_Y \times n_X$。此时，$H_j$ 是由每个观测函数对每个观测量分量求偏导后，将当前时刻状态量的先验估计代入得到：
  • 观测噪声方差 ${\sigma_R}^2$ 演变为协方差矩阵，用 $\Sigma_R$ 表示，维数为 $n_Y \times n_Y$
  • 卡尔曼增益系数 $K$ 演变为矩阵，维数为 $n_X \times n_Y$

对应的五个公式演变为：

公式 (4.3) 中 $\mu_k^+$ 即 k 时刻状态量 $X_k$ 的后验估计 $\hat{x}_k^+$，$y_k-h(\mu_k^-)$ 常被称为残差（Residual）或新息（Innovation）；公式 (4.4) 中的 $I$ 代表单位矩阵，维数为 $n_X \times n_X$。

当然，对于某个非线性系统，不一定状态转移和观测都是非线性的：

  • 线性的状态转移 + 非线性的观测
    此时，滤波递推公式由卡尔曼滤波的预测步两公式和扩展卡尔曼滤波的更新步三公式组成。

  • 非线性的状态转移 + 线性的观测
    此时，滤波递推公式由扩展卡尔曼滤波的预测步两公式和卡尔曼滤波的更新步三公式组成。

5 应用实例——基于毫米波雷达与扩展卡尔曼滤波的目标跟踪

5.1 系统分析

毫米波雷达（Radar）与激光雷达的的检测原理不同。激光雷达利用光的直线传播原理获得目标在笛卡尔坐标系下的距离信息；毫米波雷达利用电磁波的多普勒效应，获得目标在极坐标系下的距离 $\rho$、方向角 $\varphi$ 和距离变化率（径向速度）$\dot{\rho}$。

假设我们使用毫米波雷达对某作匀速直线运动的目标在笛卡尔坐标系内的横坐标 $p_x$、纵坐标 $p_y$、横向速度 $v_x$、纵向速度 $v_y$ 进行跟踪，则状态量可设为：

观测量可设为：

故，毫米波雷达的测量数据特性可用下图进行表征。

毫米波雷达的测量数据特性

其中，距离变化率 $\dot{\rho}$ 是绝对速度在径向上的投影。利用几何关系与三角函数进行推导，可知：

故，观测函数 $h(X)$ 为：

$h(X)$ 显然为非线性函数，而系统的状态转移（CV 模型）又明显是线性的，因此，此处的扩展卡尔曼滤波形式为：线性的状态转移 + 非线性的观测。

此时，滤波递推公式由卡尔曼滤波的预测步两公式和扩展卡尔曼滤波的更新步三公式组成：

卡尔曼滤波的预测步两公式

扩展卡尔曼滤波的更新步三公式

由于目标作匀速直线运动，故公式 (5.1) 中的控制项 $B*u_k$ 可以忽略。

5.2 雅可比矩阵求解

公式 (5.4) 和公式 (5.5) 中的雅可比矩阵 $H_j$ 根据定义可表示为：

直接求解 $H_j$ 较为繁琐，可通过 matlab 的内置函数 jacobian() 进行求解。创建 matlab 脚本，命名为 CalculateHj.m 并输入如下内容后保存：

1
2
3
4
5
6
7
8
9
10
11
12
13
% define symbol variables
syms px py vx vy

% state random vector
X = [px py vx vy];

% measurement functions
h1 = sqrt(px^2 + py^2);
h2 = atan(py / px);
h3 = (px * vx + py * vy) / sqrt(px^2 + py^2);

% measurement jacobian matrix
Hj = jacobian([h1 h2 h3], X)

运行 CalculateHj.m 脚本，得到 $H_j$ 如下运算结果：

1
2
3
4
5
Hj =

[                                            px/(px^2 + py^2)^(1/2),                                            py/(px^2 + py^2)^(1/2),                      0,                      0]
[                                        -py/(px^2*(py^2/px^2 + 1)),                                            1/(px*(py^2/px^2 + 1)),                      0,                      0]
[ vx/(px^2 + py^2)^(1/2) - (px*(px*vx + py*vy))/(px^2 + py^2)^(3/2), vy/(px^2 + py^2)^(1/2) - (py*(px*vx + py*vy))/(px^2 + py^2)^(3/2), px/(px^2 + py^2)^(1/2), py/(px^2 + py^2)^(1/2)]

可以直接使用上述结果，也可以对其进行化简：

5.3 代码

代码与此前的文章《（十三）手把手教你写卡尔曼滤波器》相似，区别在于状态量和观测量有所不同，且使用雅可比矩阵 $H_j$ 代替了观测矩阵 $H$（对于状态转移非线性的系统，还需使用雅可比矩阵 $F_j$ 代替状态转移矩阵 $F$）。

ExtendedKalmanFilter.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#ifndef EXTENDED_KALMAN_FILTER_
#define EXTENDED_KALMAN_FILTER_

#include "D:/eigen3/Eigen/Dense"

class ExtendedKalmanFilter
{
private:
    bool is_inited;          // flag of initialization
    bool just_begin_filt;    // flag of just begining filt data
    Eigen::VectorXd X;       // state vector
    Eigen::MatrixXd F;       // state transition matrix
    Eigen::MatrixXd P;       // state covariance matrix
    Eigen::MatrixXd Q;       // process noise covariance matrix
    Eigen::MatrixXd Hj;      // measurement jacobian matrix
    Eigen::MatrixXd R;       // measurement noise covariance matrix
    Eigen::MatrixXd K;       // extended kalman gain coefficient
    uint64_t timestamp_last; // timestamp of last frame: us
    float dt;                // delta time: s

public:
    Eigen::VectorXd Z;      // measurement vector
    uint64_t timestamp_now; // timestamp of current frame: us

public:
    ExtendedKalmanFilter();
    ~ExtendedKalmanFilter();

    inline bool &IsInited()
    {
        return is_inited;
    }

    void Init();

    inline void Reset()
    {
        is_inited = false;
    }

    inline void SetF(const Eigen::MatrixXd &f)
    {
        F = f;
    }

    inline void SetF()
    {
        F << 1.0f, 0.0f, dt, 0.0f,
            0.0f, 1.0f, 0.0f, dt,
            0.0f, 0.0f, 1.0f, 0.0f,
            0.0f, 0.0f, 0.0f, 1.0f;
    }

    inline void SetP(const Eigen::MatrixXd &p)
    {
        P = p;
    }

    inline void SetQ(const Eigen::MatrixXd &q)
    {
        Q = q;
    }

    inline void SetHj(const Eigen::MatrixXd &hj)
    {
        Hj = hj;
    }

    inline void SetR(const Eigen::MatrixXd &r)
    {
        R = r;
    }

    /**
     * @brief Predict step.
     *
     */
    void Predict();

    /**
     * @brief Update step.
     *
     */
    void Update();
};

/**
 * @brief To update measurement vector and timestamp.
 *
 * @param Z0 Z(0)
 * @param Z1 Z(1)
 * @param Z2 Z(2)
 * @param timestamp Current timestamp.
 */
void RecvRawData(float &Z0, float &Z1, float &Z2, uint32_t &timestamp)
{
    // Z0 = rho;
    // Z1 = phi;
    // Z2 = rho_dot;
    // timestamp = timestamp_new;
}

#endif

ExtendedKalmanFilter.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#include "math.h"

#include "ExtendedKalmanFilter.h"

/**
 * @brief Construct a new Extended Kalman Filter:: Extended Kalman Filter object
 *
 */
ExtendedKalmanFilter::ExtendedKalmanFilter()
    : is_inited(false),
      just_begin_filt(false),
      X(4),
      F(4, 4),
      P(4, 4),
      Q(4, 4),
      Z(3),
      Hj(3, 4),
      R(3, 3),
      K(4, 3),
      timestamp_now(0),
      timestamp_last(0),
      dt(0.0f)
{
}

/**
 * @brief Destroy the Extended Kalman Filter:: Extended Kalman Filter object
 *
 */
ExtendedKalmanFilter::~ExtendedKalmanFilter() {}

/**
 * @brief Initialize the extended kalman filter.
 *
 */
void ExtendedKalmanFilter::Init()
{
    if (!IsInited())
    {
        X << Z(0) * cos(Z(1)),
            Z(0) * sin(Z(1)),
            Z(2) * cos(Z(1)),
            Z(2) * sin(Z(1));
        P << 1.0, 0.0, 0.0, 0.0,
            0.0, 1.0, 0.0, 0.0,
            0.0, 0.0, 10.0, 0.0,
            0.0, 0.0, 0.0, 10.0;
        Q << 1.0, 0.0, 0.0, 0.0,
            0.0, 1.0, 0.0, 0.0,
            0.0, 0.0, 1.0, 0.0,
            0.0, 0.0, 0.0, 1.0;
        R << 0.09, 0.0, 0.0,
            0.0, 0.009, 0.0,
            0.0, 0.0, 0.09;
        timestamp_last = timestamp_now;
        is_inited = true;
        just_begin_filt = true;
    }
    else
        just_begin_filt = false;
}

/**
 * @brief Predict step.
 *
 */
void ExtendedKalmanFilter::Predict()
{
    if (just_begin_filt)
        return;

    dt = (timestamp_now - timestamp_last) / 1E6;
    SetF();
    timestamp_last = timestamp_now;

    // predict state vector
    X = F * X;

    // predict state covariance matrix
    P = F * P * F.transpose() + Q;
}

/**
 * @brief Update step.
 *
 */
void ExtendedKalmanFilter::Update()
{
    if (just_begin_filt)
        return;

    float c1 = X(0) ^ 2 + X(1) ^ 2;
    float c2 = sqrt(c1);
    float c3 = c1 * c2;
    if (fabs(c1 > 0.0001))
    {
        Hj << X(0) / c2, X(1) / c2, 0, 0,
            -X(1) / c2, X(0) / c2, 0, 0,
            X(1) * (X(2) * X(1) - X(3) * X(0)) / c3, X(0) * (X(3) * X(0) - X(2) * X(1)) / c3, X(0) / c2, X(1) / c2;
    }

    Eigen::VectorXd(3) hx;
    hx << c2, atan2(X(1), X(0)), (X(0) * X(2) + X(1) * X(3)) / c2;

    static Eigen::MatrixXd I = Eigen::MatrixXd::Identity(X.size(), X.size());

    // calculate extended kalman gain
    K = P * Hj.transpose() * (Hj * P * Hj.transpose() + R).inverse();

    // update state vector
    X = X + K * (Z - hx);

    // update state covariance matrix
    P = (I - K * Hj) * P;
}

main.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include "ExtendedKalmanFilter.h"

#include <iostream>

int main()
{
    ExtendedKalmanFilter ekf;

    while (1)
    {
        RecvRawData(ekf.Z(0), ekf.Z(1), ekf.Z(2), ekf.timestamp_now);
        Init();
        Predict();
        Update();
    }

    return 0;
}

点击这里下载完整工程。

6 扩展卡尔曼滤波器的优缺点

优点

扩展卡尔曼滤波与标准卡尔曼滤波有着相似的计算形式，因此运算速度同样很快。

缺点

扩展卡尔曼滤波对非线性的状态转移函数和（或）观测函数使用一阶泰勒级数展开作了线性近似，忽略了二阶及以上的高阶项，因此精度一般（通常称为一阶精度），对于高度非线性问题效果较差。此外，雅可比矩阵的计算较为繁琐，容易出错。

总结

扩展卡尔曼滤波所应用的线性近似是否具有优势主要取决于两个因素：被近似的局部非线性化程度和不确定程度。

参考

  1. b 站忠实的王大头《贝叶斯滤波与卡尔曼滤波》第十三讲：扩展卡尔曼滤波及其代码
  2. 无人驾驶技术入门（十八）| 手把手教你写扩展卡尔曼滤波器
  3. 百度百科 - 泰勒公式
猜你喜欢
Thank you for your donate!