【经典教程翻译】卡尔曼与贝叶斯滤波器：直觉理解滤波器背后的原理（下：滤波器的思考框架）

本期继续大神Roger Labbe 的 Kalman and Bayesian Filters in Python上一期 【经典教程翻译】卡尔曼与贝叶斯滤波器：直觉理解滤波器背后的原理（上），这一期介绍滤波器的一些概念，引出广义卡尔曼滤波器 g-h 滤波器的思考方式，并且通过可以实验的例子来建立 g 和 h 因子的作用。

所有文章首发于 MyEncyclopedia公众号，文章链接为

【经典教程翻译】卡尔曼与贝叶斯滤波器：直觉理解滤波器背后的原理（上）

【经典教程翻译】卡尔曼与贝叶斯滤波器：直觉理解滤波器背后的原理（下：滤波器的思考框架）

英文版原书链接在
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python

可动手的 Python Jupyter Notebook环境在

https://nbviewer.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb

MyEncyclopedia公众号以机翻为主，人肉校对为辅将这本书籍介绍给大家。由于校对时间和水平有限，可能会有翻译错误，如果大家想深入学习咀嚼，请前往原版网站查阅。

g-h 滤波器

该算法称为g-h 滤波器，又或者称为 $\alpha -\beta$ 滤波器。$g$ 和 $h$ 对应于我们在示例中使用的两个比例因子。$g$ 是我们用于测量的缩放比例（在我们的示例中为重量），$h$ 是测量值随时间变化的比例（在我们的示例中为磅/天）。$\alpha$ 和 $\beta$ 是另一种不同的名称。

g-h 滤波器是大量滤波器的基础，包括卡尔曼滤波器。换句话说，卡尔曼滤波器是 g-h 滤波器的一种形式，我将在本书后面进行证明。你可能听说过最小二乘滤波器，还有可能没听说过 Benedict-Bordner 的滤波器，都是 g-h 滤波器。每个滤波器都有不同赋值 $g$ 和 $h$ 的方式，但除此之外算法框架是相同的。例如，Benedict-Bordner 滤波器将一个限制在一定范围内的常量给 $g$ 和 $h$ 。卡尔曼等其他滤波器会在每个时间步动态地生成 $g$ 和 $h$。

让我重复以下关键点，因为它们非常重要。如果你不理解这些，你就不会理解本书的其余部分。如果你确实理解了它们，那么本书的其余部分将自然地为你展开，作为对我们将要问的各种关于 $g$ 和 $h$ “假设”问题的数学阐述。尽管数学公式可能看起来截然不同，但算法将完全相同。

• 多个数据点比一个数据点更准确，所以无论多么不准确，都不要丢弃。
• 始终选择两个数据点之间的数以创建更准确的估计。
• 根据当前估计以及我们认为它会改变多少来预测下一次测量和变化率。
• 然后选择新的估计作为预测值和下一次测量值之间的一部分，根据每个测量值的准确度进行缩放。

让我们看一下该算法的直观描述。

In [20]:

book_plots.predict_update_chart()

让我介绍一些更正式的术语。系统是我们要估计的对象。在本章中，系统就是我们试图权衡的任何东西。有些文献称此为受控体（plant）。该术语来自控制系统理论。

系统的状态是我们感兴趣的系统的当前配置或值。上个例子中我们只对重量读数感兴趣。如果我把一个100公斤的重量放在秤上，状态就是100公斤。我们根据与我们相关的内容来定义状态。刻度的颜色与我们无关，因此我们不将这些值包含在状态中。制造商的 QA 工程师可能会在状态中包含颜色，以便她可以跟踪和控制制造过程。

测量值是系统的测量值。测量值可能不准确，因此它的值可能与状态不同。

状态估计是我们的滤波器对状态的估计。例如，对于 100 公斤的重量，由于传感器误差，我们的估计可能是 99.327 公斤。这通常缩写为估计，我在本章中已经这样做了。

换句话说，状态应该被理解为系统的实际值。这个值通常对我们是隐藏的。如果我站在秤上，你就会得到一个测量值。我们称此为可观察的，因为你可以直接观察此测量。相比之下，你永远无法直接观察我的体重，你只能测量它。

这种隐藏和可观察的语言很重要。任何估计问题都包括通过可观察的测量形成对隐藏状态的估计。如果你阅读过文献，那么在定义问题时会使用这些术语，因此你需要熟悉它们。

我们使用过程模型对系统进行数学建模。在本章中，我们的过程模型假设我今天的体重是昨天的体重加上我最后一天的体重增加。过程模型不对传感器进行建模或以其他方式说明。另一个例子是汽车的过程模型。过程模型可能是“距离等于速度乘以时间”。这个模型并不完美，因为汽车的速度可以在非零时间量内变化，轮胎可以在路上打滑，等等。系统错误或过程错误是这个模型中的错误。我们永远不知道这个值；如果我们知道，我们可以改进我们的模型以达到零错误。你可能还会看到系统模型。它们的意思都是一样的。

预测步骤称为系统传播。它使用过程模型_形成新的状态估计。由于过程误差，这个估计是不完美的。假设我们随着时间的推移跟踪数据，我们说我们将状态传播到未来。一些教材称之为进化。

更新步骤称为测量更新。系统传播和测量更新的一次迭代称为时期(epoch)。

现在让我们探索几个不同的问题领域以更好地理解该算法。考虑试图在轨道上跟踪火车的问题。轨道将火车的位置限制在一个非常特定的区域。此外，火车又大又慢。他们需要很长时间才能显着地减速或加速。因此，如果我知道火车在时间 t 位于 23 公里的公里标记处并以 18 公里/小时的速度行驶，我就可以非常有信心地预测它在时间 t + 1 秒的位置。为什么这很重要？假设我们只能以±250 米以下精度测量它的位置。火车以每小时 18 公里的速度行驶，即每秒 5 米。在 t+1 秒，火车将在 23.005 公里处，但测量值可能在 22.755 公里到 23.255 公里之间的任何地方。因此，如果下一次测量显示位置在 23.4，我们知道那一定是不准确的。即使在时间 t 工程师猛踩刹车，火车仍将非常接近 23.005 公里，因为火车不能在 1 秒内减速很多。如果我们要为这个问题设计一个滤波器（我们将在本章中更进一步！），我们会想要设计一个滤波器，将高权重赋予为预测而非测量。

现在考虑跟踪抛出的球的问题。我们知道弹道物体在重力场中在真空中沿抛物线运动。但是扔到地球上的球会受到空气阻力的影响，因此它不会沿完美的抛物线运动。棒球投手在投出曲线球时利用了这一事实。假设我们正在使用计算机视觉跟踪体育场内的球，这是我在工作中所做的事情。计算机视觉跟踪的准确性可能不高，但通过假设球沿抛物线运动来预测球的未来位置也不是非常准确。在这种情况下，我们可能会设计一个滤波器，为测量和预测赋予大致相等的权重。

现在考虑尝试在飓风中追踪氦气派对气球。我们没有合法的模型可以让我们预测气球的行为，除非是在非常短的时间尺度上（例如，我们知道气球不能在 1 秒内飞 10 英里）。在这种情况下，我们将设计一个滤波器，强调测量结果而不是预测结果。

本书的大部分内容都致力于用数学方式表达最后三段中的关注点，然后让我们找到最佳解决方案（在某种数学意义上）。在本章中，我们将只是直观地给 $g$ 和 $h$ 分配不同的值，尽管这种方式不是最理想的。但基本思想是将有些不准确的测量与有些不准确的系统行为模型混合，以获得比任何一个信息源本身都更好的过滤估计。

我们可以将其表达为一种算法：

初始化

1. 初始化滤波器表示的系统状态
2. 初始化我们关于系统状态的概率分布

预测

1. 使用系统规则来预测下一时刻的状态
2. 根据预测的不确定性调整状态的概率分布

更新

1. 获取测量值并赋予其关于状态的概率分布
2. 计算估计值和测量值的残差
3. 新的估计在残差连线中某个地方

我们将在整本书中使用相同的算法，尽管有一些修改。

数学符号

我将开始介绍文献中使用的符号和变量名称。其中一些已经在上面的图表中使用过。测量值通常表示为z，（本书中的表示，一些文献使用y)。 下标 $k$ 表示时间步长，所以 $z_k$ 是这个时间步长的数据。粗体表示向量或矩阵。到目前为止，我们只考虑了一个传感器，因此只有一个传感器测量，但通常我们可能有 n 个传感器和 n 个测量。$\mathbf{x}$ 表示我们的状态，加粗表示它是一个向量。对于我们的秤示例，它代表初始重量和初始增重率，如下所示：

$$\mathbf{x}=\left[\begin{array}{c}x\\ \dot{x}\end{array}\right]$$

这里，我使用牛顿表示法，在 x 上加个点来表示速度。更准确地说，点表示 x 对时间的导数，当然是速度。对于 62 公斤的体重，每天增加 0.3 公斤，我们有

$$\mathbf{x}=\left[\begin{array}{c}62\\ 0.3\end{array}\right]$$

所以，算法很简单。状态初始化 ${\mathbf{x}}_0$，初始估计。然后我们进入一个循环，预测时间或步骤的状态 $k$ 来自时间（或步骤）的值 $k-1$。然后我们得到测量值 $z_k$ ，再选择测量和预测之间的某个中间点，创建估计 ${\mathbf{x}}_k$.

练习：编写通用算法

在上面的示例中，我明确地对此进行了编码，以解决我们在本章中一直在讨论的称重问题。例如，变量命名为“weight_scale”、“gain”等。我这样做是为了让算法更容易理解——你可以很容易地看到我们正确地实现了每一步。但是，那是为一个特定问题编写的代码，算法对于任何问题都是相同的。因此，让我们将代码重写为通用代码以解决任何问题。通用算法使用以下函数签名：

def g_h_filter ( data, x0,  dx, g,  h,  dt ):
    """
    对具有固定 g 和 h 的 1 个状态变量执行 gh 滤波器。
    'data' 包含要过滤的数据。
    'x0' 是我们状态变量的初始值
    'dx' 是我们状态变量的初始变化率
    'g' 是 g-h's 比例因子
    'h' 是 g-h's 比例因子
    'dt' 是长度时间步长
    """

将数据作为 NumPy 数组而不是列表返回。通过传入与之前相同的权重数据对其进行测试，绘制结果，并直观地确定它是否有效。

In [21]:

from kf_book.gh_internal import plot_g_h_results
def g_h_filter(data, x0, dx, g, h, dt):
    pass #  your solution here

# uncomment to run the filter and plot the results
#book_plots.plot_track([0, 11], [160, 172], label='Actual weight')
#data = g_h_filter(data=weights, x0=160., dx=1., g=6./10, h=2./3, dt=1.)
#plot_g_h_results(weights, data)

解决方案与讨论

In [22]:

import matplotlib.pylab as pylab

def g_h_filter(data, x0, dx, g, h, dt=1.):
    x_est = x0
    results = []
    for z in data:
        # prediction step
        x_pred = x_est + (dx*dt)
        dx = dx

        # update step
        residual = z - x_pred
        dx = dx + h * (residual) / dt
        x_est = x_pred + g * residual
        results.append(x_est)
    return np.array(results)

book_plots.plot_track([0, 11], [160, 172], label='Actual weight')
data = g_h_filter(data=weights, x0=160., dx=1., g=6./10, h=2./3, dt=1.)
plot_g_h_results(weights, data)
print(weights)
print(data)

[158.0, 164.2, 160.3, 159.9, 162.1, 164.6, 169.6, 167.4, 166.4, 171.0, 171.2, 172.6]
[159.2  161.8  162.1  160.78  160.985  163.311  168.1  169.696  168.204 169.164  170.892  172.629]

这应该是直截了当的。我只是将体重增加代码中的变量名称替换为变量名称x0，dx 等，没有其他改变。

选择 $g$ 和 $h$

g-h 滤波器不是一个滤波器——它是滤波器一大分类。Eli Brookner 在Tracking and Kalman Filtering Made Easy中列出了 11 个，我相信还有更多。不仅如此，每种滤波器都有许多子类型。每个滤波器的区别在于 $g$ 和 $h$ 的选择。所以我不能在这里给出“一刀切”的建议。一些滤波器设置 $g$ 和 h h