$α - β - γ$ 滤波器

本章介绍 $a l p h a - β$ 滤波器和 $α - β - γ$ 滤波器，这类滤波器常用来对时间序列数据进行平滑。 $a l p h a - β$ 滤波器和 $α - β - γ$ 滤波器在原理上和卡尔曼滤波高度相关。

示例 1 – 给金条称重

现在介绍第一个简单示例。本例对一个静态系统的状态进行估计。所谓静态系统，是指在合理时间范围内系统状态不会自发改变的系统。例如一座塔便是一个静态系统，高度便是其状态之一，它不随时间改变而变化。

本例中，我们估计一根金条的重量。假定我们用来称金条的秤是无偏的，即称重结果没有系统性偏差，但是有随机噪声。

金条就是我们所关心的系统，金条的重量就是该系统的状态。该系统的动态模型是恒定的，因为我们假定金条的重量（在短时间内）不会发生变化。

为了估计出该系统的状态（金条重量），我们可以对其进行多次称重，然后取多次测量结果的平均值。

获取本书

在时刻 $n$ ，估计值 ${\hat{x}}_{n, n}$ 便是所有之前测量的平均值：

{\hat{x}}_{n, n} = \frac{1}{n} (z_{1} + z_{2} + \dots + z_{n - 1} + z_{n}) = \frac{1}{n} \sum_{i = 1}^{n} (z_{i})

注解：

$x$	是金条重量的真值
$z_{n}$	是 $n$ 时刻对金条重量的测量值
${\hat{x}}_{n, n}$	是在 $n$ 时刻，使用了 $n$ 时刻的测量值 $z_{n}$ ，对 $x$ 的估计值
${\hat{x}}_{n + 1, n}$	是在 $n$ 时刻对未来状态（ $n + 1$ 时刻）的预测，记为 ${\hat{x}}_{n + 1, n}$ ，或者说外插
${\hat{x}}_{n - 1, n - 1}$	是在 $n - 1$ 时刻，使用了 $n - 1$ 时刻的测量值 $z_{n - 1}$ ，对 $x$ 的估计值
${\hat{x}}_{n, n - 1}$	是一个先验估计 - 在 $n - 1$ 时刻对 $n$ 时刻的系统状态所进行的预测（译注：对第n个时刻而言， ${\hat{x}}_{n, n - 1}$ 是先验估计， ${\hat{x}}_{n + 1, n}$ 是预测）

注：在本教程中，变量上的尖号符号（或者叫hat）代表这是一个对该变量的估计值。

由于金条的重量不随时间改变而改变，系统动态模型在本例中是静态的（恒定），因此有 ${\hat{x}}_{n + 1, n} = {\hat{x}}_{n, n}$ .

上面求平均的表达式虽然在数学层面是正确的，但是它不具备可实现性。这是因为根据平均值的定义，为了估计 ${\hat{x}}_{n, n}$ ，我们需要存储下所有的历史测量值，这对内存开销巨大。并且每次获得了新的测量值后都需要完全重新从第一次测量开始计算，这对CPU算力也消耗巨大。

现实一点的考虑是，最好只需存储上一时刻的估计值 ${\hat{x}}_{n - 1, n - 1}$ ，并在新的测量完成后更新它即可。下图描述了这个思路：

根据当前的测量和先验估计，估计当前的状态。
根据当前的状态估计以及系统动态模型，进行下一时刻的预测。

把上述求平均的表达式在数学上等效变换一下，可以得到：

	注解
${\hat{x}}_{n, n} = \frac{1}{n} \sum_{i = 1}^{n} (z_{i})$	求平均公式： $n$ 个测量的和再除以 $n$
$= \frac{1}{n} (\sum_{i = 1}^{n - 1} (z_{i}) + z_{n})$	前 $n - 1$ 个测量的和再加上最近一次的测量值整体除以 $n$
$= \frac{1}{n} \sum_{i = 1}^{n - 1} (z_{i}) + \frac{1}{n} z_{n}$	$\frac{1}{n}$ 乘进去展开
$= \frac{1}{n} \frac{n - 1}{n - 1} \sum_{i = 1}^{n - 1} (z_{i}) + \frac{1}{n} z_{n}$	给求和项同时乘以并除以 $n - 1$
$= \frac{n - 1}{n} \frac{1}{n - 1} \sum_{i = 1}^{n - 1} (z_{i}) + \frac{1}{n} z_{n}$	调整顺序橘黄色的项就是上一时刻的估计值
$= \frac{n - 1}{n} {\hat{x}}_{n - 1, n - 1} + \frac{1}{n} z_{n}$	把求和项用上一时刻的估计值替换掉
$= {\hat{x}}_{n - 1, n - 1} - \frac{1}{n} {\hat{x}}_{n - 1, n - 1} + \frac{1}{n} z_{n}$	把 ${\hat{x}}_{n - 1, n - 1}$ 乘进 $\frac{n - 1}{n}$ 的分子，并拆项
$= {\hat{x}}_{n - 1, n - 1} + \frac{1}{n} (z_{n} - {\hat{x}}_{n - 1, n - 1})$	提出 $\frac{1}{n}$

${\hat{x}}_{n - 1, n - 1}$ 就是在 $n - 1$ 时刻使用 $n - 1$ 时刻的测量值，对 $x$ 的状态估计。

接下来使用 $n - 1$ 时刻的估计值 ${\hat{x}}_{n - 1, n - 1}$ 计算 ${\hat{x}}_{n, n - 1}$ （对 $n$ 时刻 $x$ 的预测）。即把 ${\hat{x}}_{n - 1, n - 1}$ 外插至 $n$ 时刻。

由于系统模型是静态的，当前时刻对 $x$ 的预测就等于上一时刻对 $x$ 的估计： ${\hat{x}}_{n, n - 1} = {\hat{x}}_{n - 1, n - 1}$ .

基于上述推导，对当前时刻状态 ${\hat{x}}_{n, n}$ 的估计可以写成：

{\hat{x}}_{n, n} = {\hat{x}}_{n, n - 1} + \frac{1}{n} (z_{n} - {\hat{x}}_{n, n - 1})

上式即为卡尔曼滤波的五个方程之一。称为状态更新方程。其意为：

系数 $\frac{1}{n}$ 是本例特定的。后面会具体谈到这个系数的重要性，但此刻可以先指出，在卡尔曼滤波的语境中，这个系数被称作卡尔曼增益，符号为 $K_{n}$ . 其具有下标 $n$ 意味着卡尔曼增益随着每次迭代都会改变。

$K_{n}$ 的提出是Rudolf Kalman重要的贡献之一。

在进展到卡尔曼滤波之前，我们先不用 $K_{n}$ ，而是用希腊字母 $α_{n}$ 来表示这个系数。

所以状态更新方程可以写作：

{\hat{x}}_{n, n} = {\hat{x}}_{n, n - 1} + α_{n} (z_{n} - {\hat{x}}_{n, n - 1})

$(z_{n} - {\hat{x}}_{n, n - 1})$ 这一项被称为“测量残差”，也叫更新量。更新量包含新的信息。

本例中，随着 $n$ 的增加， $\frac{1}{n}$ 会下降。在一开始，因为没有足够的信息，第一次估计完全是基于第一次的测量值的（ $\frac{1}{n} |_{n = 1} = 1$ ）。随着迭代进行，每次后续测量的权重都在下降，并且会逐渐变得可以忽略不计。

继续讲示例。在进行第一次测量之前，我们可以根据金条上的钢印来猜测（或粗略估计）金条的重量，这是初始估计，是算法的第一个估计值。

卡尔曼滤波需要一个初始估计作为初始值，这个值可以非常粗略。

估计算法

下图描述了本例中所使用的算法。

现在，可以开始具体测量和估计的过程了。

数值示例

第0次迭代

初始化

金条重量的初始估计是1000g，这个估计仅在滤波器初始化时使用一次，后续迭代不再需要这个值。

{\hat{x}}_{0, 0} = 1000 g

预测

金条的重量不会改变，系统模型是静态的，下一个时刻的预测就等于此时的初始化估计值。

{\hat{x}}_{1, 0} = {\hat{x}}_{0, 0} = 1000 g

第1次迭代

第1步

用秤称重：

z_{1} = 996 g

第2步

计算增益。本例中 $α_{n} = \frac{1}{n}$ ，故：

α_{1} = \frac{1}{1} = 1

用状态更新方程计算当前的估计值：

{\hat{x}}_{1, 1} = {\hat{x}}_{1, 0} + α_{1} (z_{1} - {\hat{x}}_{1, 0}) = 1000 + 1 (996 - 1000) = 996 g

注：初始估计可以是任何值，由于

α_{1} = 1

，初始值在第1次迭代的时候将被抵消。

第3步

系统模型是静态的，金条的重量不应该改变，下一时刻的预测应该等于此时的估计：

{\hat{x}}_{2, 1} = {\hat{x}}_{1, 1} = 996 g

第2次迭代

在一个采样周期过后，上一时刻的预测值成为了这一时刻的先验估计。

{\hat{x}}_{2, 1} = 996 g

第1步

做第二次称重：

z_{2} = 994 g

第2步

计算增益：

α_{2} = \frac{1}{2}

计算估计值：

{\hat{x}}_{2, 2} = {\hat{x}}_{2, 1} + α_{2} (z_{2} - {\hat{x}}_{2, 1}) = 996 + \frac{1}{2} (994 - 996) = 995 g

第3步

{\hat{x}}_{3, 2} = {\hat{x}}_{2, 2} = 995 g

第3次迭代

z_{3} = 1021 g α_{3} = \frac{1}{3}

{\hat{x}}_{3, 3} = 995 + \frac{1}{3} (1021 - 995) = 1003.67 g

{\hat{x}}_{4, 3} = 1003.67 g

第4次迭代

z_{4} = 1000 g α_{4} = \frac{1}{4}

{\hat{x}}_{4, 4} = 1003.67 + \frac{1}{4} (1000 - 1003.67) = 1002.75 g

{\hat{x}}_{5, 4} = 1002.75 g

第5次迭代

z_{5} = 1002 g α_{5} = \frac{1}{5}

{\hat{x}}_{5, 5} = 1002.75 + \frac{1}{5} (1002 - 1002.75) = 1002.6 g

{\hat{x}}_{6, 5} = 1002.6 g

第6次迭代

z_{6} = 1010 g α_{6} = \frac{1}{6}

{\hat{x}}_{6, 6} = 1002.6 + \frac{1}{6} (1010 - 1002.6) = 1003.83

{\hat{x}}_{7, 6} = 1003.83 g

第7次迭代

z_{7} = 983 g α_{7} = \frac{1}{7}

{\hat{x}}_{7, 7} = 1003.83 + \frac{1}{7} (983 - 1003.83) = 1000.86 g

{\hat{x}}_{8, 7} = 1000.86 g

第8次迭代

z_{8} = 971 g α_{8} = \frac{1}{8}

{\hat{x}}_{8, 8} = 1000.86 + \frac{1}{8} (971 - 1000.86) = 997.125 g

{\hat{x}}_{9, 8} = 997.125 g

第9次迭代

z_{9} = 993 g α_{9} = \frac{1}{9}

{\hat{x}}_{9, 9} = 997.125 + \frac{1}{9} (993 - 997.125) = 996.67 g

{\hat{x}}_{10, 9} = 996.67 g

第10次迭代

z_{10} = 1023 g α_{10} = \frac{1}{10}

{\hat{x}}_{10, 10} = 996.67 + \frac{1}{10} (1023 - 996.67) = 999.3 g

{\hat{x}}_{11, 10} = 999.3 g

至此先告一段落。增益随着每次测量而减小，故后面的测量对估计值的贡献总小于前面的测量。我们已经很接近真实的金条重量了（1000g）。如果做更多次的称重，我们会和真值更加接近。

下表汇总并比较了上面的测量值、估计值以及真值。

$n$	1	2	3	4	5	6	7	8	9	10
$α_{n}$	$1$	$\frac{1}{2}$	$\frac{1}{3}$	$\frac{1}{4}$	$\frac{1}{5}$	$\frac{1}{6}$	$\frac{1}{7}$	$\frac{1}{8}$	$\frac{1}{9}$	$\frac{1}{10}$
$z_{n}$	996	994	1021	1000	1002	1010	983	971	993	1023
${\hat{x}}_{n, n}$	996	995	1003.67	1002.75	1002.6	1003.83	1000.86	997.125	996.67	999.3
${\hat{x}}_{n + 1, n}$	996	995	1003.67	1002.75	1002.6	1003.83	1000.86	997.125	996.67	999.3

结果分析

下图中列出了测量值、估计值以及真值。

Measurements vs. True value vs. Estimates

这个估计算法对测量值有平滑的效果，并且能够收敛到真值。

示例小结

本例中，我们设计了一个针对静态系统的简单估计算法。我们还推导出了五个卡尔曼滤波方程之一的状态更新方程。下一章中我们会重新审视这个状态更新方程。

示例 2 - 跟踪直线匀速运动的飞行器

是时候考虑一下状态随时间变化的动态系统了。本例中，我们尝试用 $α - β$ 滤波器对一个直线匀速飞行中的飞行器进行跟踪。

假设一个只有一个维度的世界中，有这样一个飞行器在往远离雷达的方向飞行（或者靠近雷达）。因为是一维空间，飞行器到雷达的角度是恒定的，其高度也是恒定的。

$x_{n}$ 表示 $n$ 时刻飞行器的距离。飞行器速度可以近似用距离差分法得到 - 计算距离随时间的变化率。

获取本书

因此，速度是距离的导数：

\dot{x} = v = \frac{d x}{d t}

雷达向目标的方向以固定频率发射跟踪波束，两次跟踪测量之间的时间间隔为 $Δ t$ .

则匀速运动的动力学模型可以由下面的运动方程给出：

x_{n + 1} = x_{n} + Δ t {\dot{x}}_{n}

{\dot{x}}_{n + 1} = {\dot{x}}_{n}

根据这些方程，下一个采样周期时的飞行器距离等于当前飞行器距离加上目标速度乘以采样间隔时间。由于我们假设飞行器的速度不变，下一时刻的速度等于当前时刻的速度。

上述方程称为状态外插方程 (也叫转移方程或预测方程)。

在上个示例中我们已经用过状态外插方程了，只不过上个例子中这个方程是个恒等式，即下一时刻的状态等于当前时刻的状态。

状态外插方程依赖系统动态模型，因此不同的示例中这个方程也不同。

这个方程有一个以矩阵形式给出的更加一般的形式，后续会讲到。

上述方程的形式是本例特有的。

注：我们已经学到了两个卡尔曼滤波方程：

状态更新方程
状态外插方程

现在我们来把状态更新方程改一改。

$α - β$ 滤波器

假设雷达的跟踪间隔 ( $Δ t$ ) 为5秒，假设 $n - 1$ 时刻飞行器的距离为30,000m，其速度为40m/s.

使用上述状态外插方程，我们能够预测 $n$ 时刻的目标位置为：

{\hat{x}}_{n, n - 1} = {\hat{x}}_{n - 1, n - 1} + Δ t {\hat{\dot{x}}}_{n - 1, n - 1} = 30000 + 5 \times 40 = 30200 m

预测目标 $n$ 时刻的速度为：

{\hat{\dot{x}}}_{n, n - 1} = {\hat{\dot{x}}}_{n - 1, n - 1} = 40 m / s

然而 $n$ 时刻雷达测量的目标距离 ( $z_{n}$ ) 为30,110m而非30,200m。预测和实际测量的距离之间相差了90m. 这个差有可能是两个原因导致的：

雷达测量不够精准
飞行器速度变化了。新的速度是 $\frac{30, 110 - 30, 000}{5} = 22 m / s$

哪个原因是正确的呢？

我们把速度的状态更新方程写下来：

{\hat{\dot{x}}}_{n, n} = {\hat{\dot{x}}}_{n, n - 1} + β (\frac{z_{n} - {\hat{x}}_{n, n - 1}}{Δ t})

系数 $β$ 的值和雷达的测量精度等级有关。假设雷达的 $1 σ$ 精度是20m，那么90m的误差大概率是飞行器速度改变了，我们应该把 $β$ 的值调高一些。如果 $β$ 设为0.9，那么此时估计的速度就应该是：

{\hat{\dot{x}}}_{n, n} = {\hat{\dot{x}}}_{n, n - 1} + β (\frac{z_{n} - {\hat{x}}_{n, n - 1}}{Δ t}) = 40 + 0.9 (\frac{30110 - 30200}{5}) = 23.8 m / s

而另一方面，假如雷达的 $1 σ$ 精度是150m，那90m的误差大概率是雷达测得不准，我们应该把 $β$ 的值降低一些。如果 $β$ 的值降为0.1，那么估计的速度就应该是：

{\hat{\dot{x}}}_{n, n} = {\hat{\dot{x}}}_{n, n - 1} + β (\frac{z_{n} - {\hat{x}}_{n, n - 1}}{Δ t}) = 40 + 0.1 (\frac{30110 - 30200}{5}) = 38.2 m / s

如果是飞行器真实速度从40m/s降到了22m/s，可以看到在10个测量周期之后（把上面的公式以 $β$ = 0.1 带入10次），飞行器速度的估计也逐渐降到了22m/s. 而如果是因为雷达测量不准，则后续测量的位置将会大致均匀散布在真实位置前后，整体上计算出来的平均速度会保持在40m/s左右不变。

飞行器位置的状态更新方程与上一个例子里的方程类似：

{\hat{x}}_{n, n} = {\hat{x}}_{n, n - 1} + α (z_{n} - {\hat{x}}_{n, n - 1})

和上例不同的是，上例里 $α$ 系数每周期都在重新计算 ( $α_{n} = \frac{1}{n}$ )，而本例中 $α$ 则是恒定的。

$α$ 系数的大小和雷达精度有关。对高精度雷达，应该选用高的 $α$ ，以给测量值分配更高的权重。如果 $α = 1$ ，则估计的飞行器距离会等于测量值。

{\hat{x}}_{n, n} = {\hat{x}}_{n, n - 1} + 1 (z_{n} - {\hat{x}}_{n, n - 1}) = z_{n}

如果 $α = 0$ ，则测量值完全起不到任何作用：

{\hat{x}}_{n, n} = {\hat{x}}_{n, n - 1} + 0 (z_{n} - {\hat{x}}_{n, n - 1}) = {\hat{x}}_{n, n - 1}

于是我们推导出了雷达跟踪问题的状态更新方程。这又叫 $α - β$ 跟踪更新方程或 $α - β$ 跟踪滤波方程。

位置的状态更新方程为：

{\hat{x}}_{n, n} = {\hat{x}}_{n, n - 1} + α (z_{n} - {\hat{x}}_{n, n - 1})

速度的状态更新方程为：

{\hat{\dot{x}}}_{n, n} = {\hat{\dot{x}}}_{n, n - 1} + β (\frac{z_{n} - {\hat{x}}_{n, n - 1}}{Δ t})

注：在其他的书中，

α - β

滤波器又叫 g-h 滤波器，字母g代替了希腊字母

α

，字母h代替了希腊字母

β

注：本例中，我们是根据雷达距离测量来间接估计的飞行器速度的 (

\dot{x} = \frac{Δ x}{Δ t}

)，实际上现代雷达可以直接用多普勒效应测量径向速度。但目前我们的目的是解释卡尔曼滤波而非雷达工作原理，所以为了简单起见，我们的示例中将全部用距离观测来间接测量速度。

估计算法

下图描述了本例所使用的估计算法。

与上一个示例不同，增益（ $α$ 和 $β$ ）的值在本例里是给定的。在卡尔曼滤波里， $α$ 和 $β$ 会被卡尔曼增益代替，并且每个采样周期会重新计算，后面会讲到。

现在来看一个数值示例。