α−β\alpha-\betaα−β滤波器

文章目录

α−β\alpha-\betaα−β滤波器
- 状态方程和观测方程
- α−β\alpha-\betaα−β滤波器最终形式
- α−β\alpha-\betaα−β滤波器作为一种稳态Kalman滤波器
- - 状态方程、量测方程
  - 推导增益KKK
  - 结论
- α,β\alpha,\betaα,β参数的工程化设置
- 仿真
- - 工况1 初始估计误差很大
  - 工况2 观测噪声很大
  - 工况3 恒定加速度扰动
  - 工况4 增大a、b的结果
- 总结
- Ref

本博客来自一个课程的作业，包含以下内容：

一．α-β滤波器推导
二．作为稳态滤波器与Kalman滤波的比较
三．α-β的参数设置及其最优性
四．仿真演示
五．总结

KF是噪声干扰下、线性系统最优的滤波器，而a-b滤波器是KF的衍生物，在性质、结构上都很类似，但在计算复杂性上大大降低。KF动态地调整系数，以最大效率的滤除观测误差，使对应的状态估计值最佳。而a-b滤波则使用固定的滤波增益，或者认为给定好的变化增益。所以a-b滤波器是次优的。通过调整系数，用户可以控制用于从最近的测量中历史数据的长度，并用来估计数据的真实值。

对一标量x(k)x(k)x(k)，存在对它的一系列观测值{z(k)}\{z(k)\}{z(k)}，则α−β\alpha-\betaα−β滤波器的预测和滤波方程为
x(k+1∣k)=x^(k)+Δt⋅v^(k)x^(k+1)=x(k+1∣k)+α⋅[z(k+1)−x(k+1∣k)]v^(k+1)=v^(k)+β/Δt⋅[z(k+1)−x(k+1∣k)\begin{aligned} x(k+1|k)&=\hat x(k) + \Delta t\cdot\hat v(k)\\ \hat x(k+1)&=x(k+1|k)+\alpha \cdot[z(k+1)-x(k+1|k)]\\ \hat v(k+1)&=\hat v(k)+\beta/\Delta t \cdot [z(k+1)-x(k+1|k) \end{aligned} x(k+1∣k)x^(k+1)v^(k+1)=x^(k)+Δt⋅v^(k)=x(k+1∣k)+α⋅[z(k+1)−x(k+1∣k)]=v^(k)+β/Δt⋅[z(k+1)−x(k+1∣k)其中Δt\Delta tΔt为观测的更新时间。参数α,β\alpha,\betaα,β的设置依赖于对数据机动量噪声估计的大小、观测误差的大小。

状态方程和观测方程

a-b滤波和a-b-g滤波都要解决这样一个问题：对1维数据x(k)x(k)x(k)，已知这个数据的一系列观测值{x(0),x(1),⋯,x(k)}\{x(0),x(1),\cdots,x(k)\}{x(0),x(1),⋯,x(k)}，然后要利用已有数据进行滤波估计，得到更准确和平滑的值。

通常使用场景：假设1维场景的目标以匀速直线运动，而能够测量到距离序列，需要估计它的真实距离。

a-b滤波器不要求已知数学模型，它所考虑的只有状态在时域的连续性，也就是状态量xxx的一阶导和二阶导存在如下关系
x˙=vv˙=a\dot x=v \\ \dot v=a x˙=vv˙=a这个式子如果改成离散型，一阶或二阶形式分别就是a-b滤波，a-b-g滤波的系统方程
x(k+1)=x(k)+Δt⋅v(k)x(k+1)=x(k)+Δt⋅v(k)+Δt22⋅a(k)x(k+1)=x(k)+\Delta t\cdot v(k) \\ x(k+1)=x(k)+\Delta t\cdot v(k)+ \frac{\Delta t^2}{2}\cdot a(k) x(k+1)=x(k)+Δt⋅v(k)x(k+1)=x(k)+Δt⋅v(k)+2Δt2⋅a(k)改写成state-space的离散形式：
x(k+1∣k)=Φ(k+1∣k)x(k)2个状态：位置-速度：Φ(k+1∣k)=[1Δt01],x=[xv]3个状态：位置-速度-加速度：Φ(k+1∣k)=[1ΔtΔt2201Δt001],x=[xva]{\boldsymbol x}(k+1|k)=\boldsymbol \Phi(k+1|k){\boldsymbol x}(k)\\ \text{2个状态：位置-速度：} \boldsymbol \Phi(k+1|k)=\begin{bmatrix} 1 & \Delta t \\ 0 & 1 \end{bmatrix} , \boldsymbol{x}=\begin{bmatrix} x \\ v \end{bmatrix}\\ \text{3个状态：位置-速度-加速度：} \boldsymbol \Phi(k+1|k)=\begin{bmatrix} 1 & \Delta t & \frac{\Delta t^2}{2}\\ 0 & 1 & \Delta t\\ 0 & 0& 1 \end{bmatrix} , \boldsymbol{x}=\begin{bmatrix} x \\ v \\ a \end{bmatrix} x(k+1∣k)=Φ(k+1∣k)x(k)2个状态：位置-速度：Φ(k+1∣k)=[10Δt1],x=[xv]3个状态：位置-速度-加速度：Φ(k+1∣k)=100Δt102Δt2Δt1,x=xva而观测方程往往是只有位置，只有1维，对a-b滤波而言
z(k)=[10][x(k)v(k)]z(k)= \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} x(k)\\ v(k) \end{bmatrix} z(k)=[10][x(k)v(k)]对a-b-g滤波而言，
z(k)=[100][x(k)v(k)a(k)]z(k)= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} x(k)\\ v(k)\\ a(k) \end{bmatrix} z(k)=[100]x(k)v(k)a(k)

如果你的实际模型和这个不一样，很显然，它是被一个2阶Markov过程近似了，而其他的内部关系都没有考虑，都被当做了噪声。也就是说，噪声项w(k)w(k)w(k)直接施加在加速度上，以3阶为例：
x(k+1)=[1ΔtΔt2201Δt001][x(k)v(k)a(k)]+[Δt22Δt1]w(k)\boldsymbol{x}(k+1)= \begin{bmatrix} 1 & \Delta t & \frac{\Delta t^2}{2}\\ 0 & 1 & \Delta t\\ 0 & 0& 1 \end{bmatrix} \begin{bmatrix} x(k) \\ v(k) \\ a(k) \end{bmatrix} + \begin{bmatrix} \frac{\Delta t^2}{2} \\ \Delta t \\ 1 \end{bmatrix}w(k) x(k+1)=100Δt102Δt2Δt1x(k)v(k)a(k)+2Δt2Δt1w(k)

z(k)=[100][x(k)v(k)a(k)]+n(k)z(k)= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} x(k)\\ v(k)\\ a(k) \end{bmatrix} + n(k) z(k)=[100]x(k)v(k)a(k)+n(k)其中w(k),n(k)w(k),n(k)w(k),n(k)均是不相关的高斯白噪声，互相独立，设他们的方差为Q=σw2,R=σn2Q=\sigma_w^2,R=\sigma_n^2Q=σw2,R=σn2. 3阶状态方程的Γ=[Δt2/2,Δt,1]T\boldsymbol \Gamma =[\Delta t^2/2, \Delta t, 1]^{\mathrm T}Γ=[Δt2/2,Δt,1]T，2阶方程的噪声增益为Γ=[Δt2/2,Δt]T\boldsymbol \Gamma =[\Delta t^2/2, \Delta t]^{\mathrm T}Γ=[Δt2/2,Δt]T

α−β\alpha-\betaα−β滤波器最终形式

需要从已有的位置观测序列{z(0),z(1),⋯,z(k),z(k+1)}\{z(0),z(1),\cdots,z(k),z(k+1)\}{z(0),z(1),⋯,z(k),z(k+1)}中估计位置x^(k+1)\hat x(k+1)x^(k+1)、速度v^(k+1)\hat v(k+1)v^(k+1)。类似地a-b-g滤波器需要从已有的位置观测序列{z(0),z(1),⋯,z(k),z(k+1}\{z(0),z(1),\cdots,z(k),z(k+1\}{z(0),z(1),⋯,z(k),z(k+1}中估计位置x^(k+1)\hat x(k+1)x^(k+1)、速度v^(k+1)\hat v(k+1)v^(k+1)、加速度a^(k+1)\hat a(k+1)a^(k+1)，实际上由于种种原因，估计出来的导数项如速度、加速度都是不可用的，但滤波方程中仍然包含这些项，以保证连续性。

类似于KF的滤波过程，分为预测和更新2部分。利用上一轮的滤波数据，其中包含在新一轮观测更新时，新的预测为
x(k+1∣k)=x^(k)+Δt⋅v^(k)v(k+1∣k)=v^(k)x(k+1|k)=\hat x(k) + \Delta t\cdot\hat v(k)\\ v(k+1|k)=\hat v(k) x(k+1∣k)=x^(k)+Δt⋅v^(k)v(k+1∣k)=v^(k)新的测量为z(k+1)z(k+1)z(k+1)，则观测误差为
z~(k+1)=z(k+1)−x(k+1∣k)\tilde z(k+1)=z(k+1)-x(k+1|k) z~(k+1)=z(k+1)−x(k+1∣k)

滤波方程使用残差的 α 倍来校正位置估计，使用残差的 β/T 倍来校正速度估计，即最终结果为
x^(k+1)=x(k+1∣k)+α⋅z~(k+1)v^(k+1)=v(k+1∣k)+β/Δt⋅z~(k+1)\hat x(k+1)=x(k+1|k)+\alpha \cdot\tilde z(k+1)\\ \hat v(k+1)=v(k+1|k)+\beta/\Delta t \cdot \tilde z(k+1) x^(k+1)=x(k+1∣k)+α⋅z~(k+1)v^(k+1)=v(k+1∣k)+β/Δt⋅z~(k+1)其中第二个式子也通常把Δt\Delta tΔt直接引进去。a-b滤波的参数α,β\alpha,\betaα,β需要实际调整。

α−β\alpha-\betaα−β滤波器作为一种稳态Kalman滤波器

状态方程、量测方程

预测状态x^(k+1∣k)=Φ(k+1∣k)x^(k)\hat{\boldsymbol x}(k+1|k)=\boldsymbol \Phi(k+1|k)\hat{\boldsymbol x}(k) x^(k+1∣k)=Φ(k+1∣k)x^(k)
修正x^(k+1∣k+1)=x^(k+1∣k)+K[z(k+1)−Hx^(k+1∣k)]\hat{\boldsymbol x}(k+1|k+1)=\hat{\boldsymbol x}(k+1|k)+\boldsymbol K[z(k+1)-\boldsymbol H\hat{\boldsymbol x}(k+1|k)] x^(k+1∣k+1)=x^(k+1∣k)+K[z(k+1)−Hx^(k+1∣k)]其中观测矩阵h=[1,0]\boldsymbol h=[1,0]h=[1,0]，Kalman增益是一个定常矩阵
K=[αβ/Δt]\boldsymbol K=\begin{bmatrix} \alpha \\ \beta/\Delta t \end{bmatrix} K=[αβ/Δt]这2个方程与Kalman滤波完全一样。

推导增益KKK

根据Kalman滤波协方差矩阵稳态的条件，可以推出其余的3个方程：

状态协方差预测

P(k+1∣k)=ΦP(k∣k)ΦT+Γσw2ΓT\boldsymbol P(k+1|k)=\boldsymbol\Phi \boldsymbol P(k|k) \boldsymbol\Phi^{\mathrm T}+ \boldsymbol\Gamma \sigma^2_w \boldsymbol\Gamma^{\mathrm T} P(k+1∣k)=ΦP(k∣k)ΦT+Γσw2ΓT

目标跟踪增益（方括号里的矩阵求逆退化为1维标量求倒数）

K(k+1)=P(k+1∣k)hT[hP(k+1∣k)hT+σn2]−1\boldsymbol K(k+1) =\boldsymbol P(k+1|k)\boldsymbol h^{\mathrm T} [\boldsymbol h \boldsymbol P(k+1|k)\boldsymbol h^{\mathrm T}+\sigma_n^2]^{-1} K(k+1)=P(k+1∣k)hT[hP(k+1∣k)hT+σn2]−1

协方差更新

P(k+1∣k+1)=[I2−K(k+1)h]P(k+1∣k)\boldsymbol P(k+1|k+1)= [I_2-\boldsymbol K(k+1)\boldsymbol h] \boldsymbol P(k+1|k) P(k+1∣k+1)=[I2−K(k+1)h]P(k+1∣k)
构造这样的关系，即预测出的协方差矩阵已经趋于稳态：
P(k+1∣k)=P(k∣k−1)=P\boldsymbol P(k+1|k)=\boldsymbol P(k|k-1)=\boldsymbol P P(k+1∣k)=P(k∣k−1)=P将稳态形式的Kalman滤波的三个方程代入，得到Riccati方程形式
Φ[P−PH⊤(HPH⊤+R)−1HP]Φ⊤+ΓQΓ⊤=P↓ΦPΦ⊤−ΦPH⊤(HPHT+σn2)−1HPΦ⊤+Γσw2Γ⊤=P\Phi\left[P-P H^{\top}\left(H P H^{\top}+R\right)^{-1} H P\right] \Phi^{\top}+\Gamma Q \Gamma^{\top}=P\\ \downarrow\\ \Phi P \Phi^{\top}-\Phi P H^{\top}\left( H P H^{\mathrm T}+\sigma_n^2\right)^{-1} H P \Phi^{\top}+\Gamma \sigma_w^2 \Gamma^{\top}=P Φ[P−PH⊤(HPH⊤+R)−1HP]Φ⊤+ΓQΓ⊤=P↓ΦPΦ⊤−ΦPH⊤(HPHT+σn2)−1HPΦ⊤+Γσw2Γ⊤=P求解预测的协方差矩阵
P=[P11P12P21P22]\boldsymbol P=\begin{bmatrix} P_{11} & P_{12} \\ P_{21} & P_{22} \end{bmatrix} P=[P11P21P12P22]通过将所有的项代入，可以在mathematica中解出3个方程的解，共有4组可行解，再考虑到协方差矩阵的正定性，只剩下一组可行解。

即作为其他几个变量的函数
P11↦(Δt,σw,σn)P12↦(Δt,σw,σn)P22↦(Δt,σw,σn)P_{11}\mapsto(\Delta t,\sigma_w,\sigma_n)\\ P_{12}\mapsto(\Delta t,\sigma_w,\sigma_n)\\ P_{22}\mapsto(\Delta t,\sigma_w,\sigma_n) P11↦(Δt,σw,σn)P12↦(Δt,σw,σn)P22↦(Δt,σw,σn)

将这个表达式代入，就可以得到滤波增益：
K=[αβ/Δt]=1σn2+P11[P11P21]\boldsymbol K=\begin{bmatrix} \alpha \\ \beta/\Delta t \end{bmatrix} =\frac{1}{\sigma_n^2+P_{11}} \begin{bmatrix} P_{11}\\ P_{21} \end{bmatrix} K=[αβ/Δt]=σn2+P111[P11P21]这样，就可以构造关系：
Δt,σw,σn↦P11,P12↦α,β\Delta t,\sigma_w,\sigma_n \mapsto P_{11},P_{12} \mapsto \alpha,\beta Δt,σw,σn↦P11,P12↦α,β但是，由于前一个表达式过于复杂，不可能把这么长的表达式代入。因此求解结果仍然没有应用价值。

结论

结论：a-b滤波器是Kalman滤波器作为稳态滤波器的衍生物，在性质、结构上与KF都很类似。

KF有明确的状态空间模型，而a-b滤波器将状态简化为二阶噪声驱动的Markov随机过程
KF有多维观测量，而a-b滤波器只有所有状态的第1维。但两者都是将观测直接加白噪声。
KF预测协方差矩阵，而a-b滤波器只关注状态的连续性，不考虑不确定部分
KF是噪声干扰下、线性系统最优的滤波器，而a-b滤波器是稳态的固定系数的，次优的滤波器

α,β\alpha,\betaα,β参数的工程化设置

KF动态地调整系数，以最大效率的滤除观测误差，使对应的状态估计值最佳。而a-b滤波则使用固定的滤波增益，或者认为给定好的变化增益。

文献给出了一种根据物理意义明确的设置方法，定义了跟踪指数tracking index，
Λ≜Δt2⋅σwσn\Lambda\triangleq \frac{\Delta t^2\cdot\sigma_w}{\sigma_n} Λ≜σnΔt2⋅σw其实就是位置的机动（加速度上的噪声w(k)w(k)w(k)）误差除以测量误差（测量方程噪声n(k)n(k)n(k)），作为模型误差与观测误差的比值，根据它，可以直接约束α\alphaα，进而确定β\betaβ。

如果观测不可靠，则跟踪指数小，会更倾向于已有的数据；
如果目标机动量大，则跟踪指数大，则更倾向于观测的数据。
跟踪指数小，滤波结果更平滑，但收敛更慢。但是不至于导致发散。

按照它，最优的参数准则应该满足以下关系：
β21−α=Λ2\frac{\beta^{2}}{1-\alpha}=\Lambda^{2} 1−αβ2=Λ2经过在Kalman滤波器框架下的推导，证明等价于最小方差的最优稳态滤波器，实际上就是Riccati方程的化简解。化简解还有
r=4+Λ−8Λ+Λ24,α=1−r2β=2(2−α)−41−αr=\frac{4+\Lambda-\sqrt{8 \Lambda+\Lambda^{2}}}{4}\ ,\ \alpha=1-r^{2}\\ \beta=2(2-\alpha)-4 \sqrt{1-\alpha} r=44+Λ−8Λ+Λ2 , α=1−r2β=2(2−α)−41−α根据跟踪指数可以获得最优的α\alphaα，进而得到最优的β\betaβ，两个参数之间的关系如下图所示。

此外，Kalata还给出了其他的一些类型如α\alphaα滤波器，α−β−γ\alpha-\beta-\gammaα−β−γ滤波器的参数设置。

由于根号内部不能小于0，所以还有一些不等式条件：
0<α<10<β≤20<4−2α−β\begin{aligned} &0<\alpha<1 \\ &0<\beta \leq 2 \\ &0<4-2 \alpha-\beta \end{aligned} 0<α<10<β≤20<4−2α−β

Kalata还给出了滤波过程的初始化，以及性能指标的预估。虽然这几个滤波器都是全局收敛的。：

仿真

搭建一个1维模拟的模块，模拟匀速直线运动，加速度项用噪声模拟。滤波器根据观测噪声的比值，用最优的参数配置。


观测更新率	dTdTdT	10Hz
观测误差（高斯白噪声，1σ\sigmaσ）	σn\sigma_nσn	1m
目标实际运动
机动性模拟（高斯白噪声，1σ\sigmaσ）	σw\sigma_wσw	0.03m/s2
加速度	a(0)a(0)a(0)	0
速度	v(0)v(0)v(0)	10m/s
初始距离	x(0)x(0)x(0)	980m
滤波器
初始距离估计	x^(0)\hat x(0)x^(0)	1200m
初始位置估计	v^(0)\hat v(0)v^(0)	5m/s
加速度估计		无
跟踪指数	Δt2⋅σwσn\frac{\Delta t^2\cdot\sigma_w}{\sigma_n}σnΔt2⋅σw	3e-4
	α\alphaα	0.0242
	β\betaβ	0.0003

工况1 初始估计误差很大

目标机动量较小，基本沿匀速直线运动，观测噪声一般，但初始估值较差。结果如下

估计太差，而观测器显然不会出现这样的问题，所以此时收敛非常慢。要收敛正常，可增大增益a、b，或者按照跟踪指数的原则，增大Λ\LambdaΛ，得到

α\alphaα 0.1362 β\betaβ 0.0099

此时，收敛更快

工况2 观测噪声很大

目标偏离直线运动较多σw=1m/s2\sigma_w=1m/s^2σw=1m/s2，，观测噪声很大σn=3.16m\sigma_n=3.16mσn=3.16m，初始估值较差，此时滤波结果仍为平滑的

工况3 恒定加速度扰动

目标短暂地加速，加速度a=10m/s2a=10m/s^2a=10m/s2，机动量噪声较大σw=1m/s2\sigma_w=1m/s^2σw=1m/s2，观测噪声一般σn=1m\sigma_n=1mσn=1m，初始估值较差，此时

工况4 增大a、b的结果