添加过程噪声

本章中,我们给一维卡尔曼滤波模型加入过程噪声。

一维卡尔曼滤波的完整模型

过程噪声

真实世界中,系统动力模型总是有不确定性的。比如我们想测量一个电阻的阻值,我们假设它是不变的,即阻值不随测量过程而改变。但实际上阻值会随着环境温度的改变而轻微改变。再比如用雷达追踪弹道导弹时,导弹动态模型的不确定性会包含一些随机的加减速。对于飞行器之类的目标,模型不确定性更大,因为飞行员随时可能进行机动。

另一方面,当我们用GPS接收机计算一个固定物体的位置时,由于固定物体不会动所以动态模型不确定性为0. 动态模型的不确定性称为 过程噪声。一些文献中也叫模型噪声、驱动噪声、动态噪声或系统噪声。过程噪声也会带来估计误差。

上一个例子里,我们估计了大楼的高度。由于大楼高度客观上不会改变,我们没有考虑过程噪声。

过程噪声方差q 表示。

协方差外插方程 应包含 过程噪声方差

恒定动态模型的协方差外插方程为:

pn+1,n=pn,n+qn

下面为加入了过程噪声后的一维卡尔曼滤波方程:

方程 方程名 文献中可能出现的其他名字
状态更新 Kn=pn,n1pn,n1+rn 卡尔曼增益 Kalman Gain 权重方程 Weight Equation
x^n,n= x^n,n1+Kn(znx^n,n1) 状态更新 State Update 滤波方程 Filtering Equation
pn,n= (1Kn)pn,n1 协方差更新 Covariance Update 修正方程Corrector Equation
状态预测 x^n+1,n=x^n,n
(恒定动态模型)

x^n+1,n=x^n,n+Δtx˙^n,n
x˙^n+1,n=x˙^n,n
(速度恒定动态模型)
状态外插 State Extrapolation 预测器方程 Predictor Equation
转移方程 Transition Equation
预测方程 Prediction Equation
动态模型 Dynamic Model
状态空间模型 State Space Model
pn+1,n=pn,n+qn
(恒定动态模型)

pn+1,nx=pn,nx+Δt2pn,nv
pn+1,nv=pn,nv+qn
(速度恒定动态模型)
协方差外插 Covariance Extrapolation 预测器协方差方程 Predictor Covariance Equation
注1:状态外插方程和协方差外插方程都依赖系统动态模型。
注2:表所描述的卡尔曼滤波方程是针对特定问题的特定形式。其一般形式在后续会以矩阵形式给出。目前我们意在理解卡尔曼滤波的基本原理和概念。

示例 6 - 估计缸中液体的温度

我们想估计一个水缸里的液体的温度。

Estimating the liquid temperature

假设在稳态情况下,液体温度是不变的。但是因为一些微小扰动可能导致液体温度的变化,我们可以用下面的方程描述系统的动态模型:

xn=T+wn

式中:

T 是恒定不变的温度

wn 是一个方差为 q 的随机过程噪声

数值示例

  • 假设真实温度是50摄氏度。
  • 假设系统模型准确,我们设过程噪声的方差为0.0001.
  • 测量误差(标准差)是0.1摄氏度。
  • 每5s测量一次。
  • 每次测量时真实的液体温度分别为:50.005oC, 49.994oC, 49.993oC, 50.001oC, 50.006oC, 49.998oC, 50.021oC, 50.005oC, 50oC, and 49.997oC.
  • 温度的测量值为:49.986oC, 49.963oC, 50.09oC, 50.001oC, 50.018oC, 50.05oC, 49.938oC, 49.858oC, 49.965oC, and 50.114oC.

下表比较了真实液体温度和测量值。

True temperature vs. measurements

第0次迭代

在第一次迭代之前,我们必须先初始化卡尔曼滤波器,并预测下一个状态(即第一个状态)

初始化

我们不知道液体的真实温度是多少,猜测是60oC.

x^0,0=60oC

我们的猜测是不准确的,所以设它的误差 σ 是100. 则初始化的估计方差(σ2)

p0,0=1002=10,000

这个方差非常大(滤波器收敛则很慢)。如果能更可靠地初始化的话,卡尔曼滤波会收敛得相对快一点。

预测

由于我们的模型是恒定的,状态预测等于当前的估计:

x^1,0=60oC

状态预测方差的外插则为:

p1,0=p0,0+q=10000+0.0001=10000.0001

第1次迭代

第1步 - 测量

测量值:

z1= 49.986oC

由于测量误差( σ )是0.1,方差( σ2 )就是0.01:

r1=0.01

第2步 - 更新

计算卡尔曼增益:

K1=p1,0p1,0+r1=10000.000110000.0001+0.01=0.999999

卡尔曼增益几乎就是1. 故我们的预测误差远大于测量误差,预测的权重可以忽略不计,而测量的权重几乎就是1.

估计当前状态:

x^1,1= x^1,0+K1(z1x^1,0)=60+0.999999(49.98660)=49.986oC

更新状态估计方差:

p1,1= (1K1)p1,0=(10.999999)10000.0001=0.01

第3步 - 预测

由于系统模型恒定,液体温度不变:

x^2,1=x^1,1=49.986oC

状态预测方差的外插为:

p2,1=p1,1+q=0.01+0.0001=0.0101

第2次迭代

第1步 - 测量

测量值:

z2= 49.963oC

由于测量误差( σ )是0.1,方差( σ2 )就是0.01:

r2=0.01

第2步 - 更新

计算卡尔曼增益:

K2=p2,1p2,1+r2=0.01010.0101+0.01=0.5

卡尔曼增益是0.5,即预测和测量的权重相等。

估计当前状态:

x^2,2= x^2,1+K2(z2x^2,1)=49.986+0.5(49.96349.986)=49.974oC

更新状态估计方差:

p2,2= (1K2)p2,1=(10.5)0.0101=0.005

第3步 - 预测

由于系统模型恒定,液体温度不变:

x^3,2=x^2,2=49.974oC

状态预测方差的外插为:

p3,2=p2,2+q=0.005+0.0001=0.0051

第3-10次迭代

后续迭代的计算结果汇总在下面的表里:

n zn 当前状态估计 ( Kn , x^n,n , pn,n ) 预测 ( x^n+1,n , pn+1,n )
3 50.09oC
K3=0.00510.0051+0.01=0.3388
x^3,3= 49.974+0.3388(50.0949.974)=50.016oC
p3,3=(10.3388)0.0051=0.0034
x^4,3=x^3,3=50.016oC
p4,3=0.0034+0.0001=0.0035
4 50.001oC
K4=0.00350.0035+0.01=0.2586
x^4,4= 50.016+0.2586(50.00150.016)=50.012oC
p4,4=(10.2586)0.0035=0.0026
x^5,4=x^4,4=50.012oC
p5,4=0.0026+0.0001=0.0027
5 50.018oC
K5=0.00270.0027+0.01=0.2117
x^5,5=50.012+0.2117(50.01850.012)=50.013oC
p5,5=(10.2117)0.0027=0.0021
x^6,5=x^5,5=50.013oC
p6,5=0.0021+0.0001=0.0022
6 50.05oC
K6=0.00220.0022+0.01=0.1815
x^6,6= 50.013+0.1815(50.0550.013)=50.02oC
p6,6=(10.1815)0.0022=0.0018
x^7,6=x^6,6=50.02oC
p7,6=0.0018+0.0001=0.0019
7 49.938oC
K7=0.00190.0019+0.01=0.1607
x^7,7= 50.02+0.1607(49.93850.02)=50.007oC
p7,7=(10.1607)0.0019=0.0016
x^8,7=x^7,7=49.978oC
p8,7=0.0016+0.0001=0.0017
8 49.858oC
K8=0.00170.0017+0.01=0.1458
x^8,8=50.007+0.1458(49.85850.007)=49.985oC
p8,8=(10.1458)0.0017=0.0015
x^9,8=x^8,8=49.985oC
p9,8=0.0015+0.0001=0.0016
9 49.965oC
K9=0.00160.0016+0.01=0.1348
x^9,9= 49.985+0.1348(49.96549.985)=49.982oC
p9,9=(10.1348)0.0016=0.0014
x^10,9=x^9,9=49.982oC
p10,9=0.0014+0.0001=0.0015
10 50.114oC
K10=0.00150.0015+0.01=0.1265
x^10,10= 49.982+0.1265(50.11449.982)=49.999oC
p10,10=(10.1265)0.0015=0.0013
x^11,10=x^10,10=49.999oC
p11,10=0.0013+0.0001=0.0014

结果分析

下图给出了卡尔曼增益的变化。

The Kalman Gain

可见,卡尔曼增益逐渐下降,故卡尔曼滤波收敛了。

下图比较了真值、测量值和估计值。图中置信区间为95%。

置信区间的计算可以在这里找到介绍。

True value, measured values and estimates

可见,估计值逐渐收敛到真值了,只是其对应的95%置信区间有些过大,说明估计不确定性较高。

示例小结

我们用一维卡尔曼滤波测量了液体的温度。尽管系统动态模型含有随机过程噪声,卡尔曼滤波也给出了较为不错的估计。

示例 7 - 估计加热中的液体的温度

类似上一个示例,我们估计缸中的液体的温度。本例里,系统的动态模型不是恒定的 - 液体在以0.1oC每秒的速度加热。

数值示例

卡尔曼滤波参数和上一个示例类似:

  • 假设模型精确,设过程噪声方差 ( q ) 为 0.0001.
  • 测量误差(标准差)为0.1oC.
  • 每5s测量一次。
  • 系统动态模型为恒定。

注:尽管真实的动态模型不是恒定的(液体在加热中),我们暂且认为系统动态模型是恒定的(温度不变)。

  • 测量点处的真实液体温度分别为:50.505oC, 50.994oC, 51.493oC, 52.001oC, 52.506oC, 52.998oC, 53.521oC, 54.005oC, 54.5oC, and 54.997oC.
  • 测量值为:50.486oC, 50.963oC, 51.597oC, 52.001oC, 52.518oC, 53.05oC, 53.438oC, 53.858oC, 54.465oC, and 55.114oC.

下图对比了真实液体温度和测量值。

True temperature vs. measurements

第0次迭代

第0次迭代和上例类似。

第1次迭代前,我们必须初始化卡尔曼滤波,并且预测下一个时刻的状态(即第1个状态)

初始化

我们不知道缸中液体的真是温度,猜测是10oC.

x^0,0=10oC

这个猜测是不准确的,于是设估计误差 ( σ ) 为 100。初始化的估计方差就是误差方差 (σ2)

p0,0=1002=10,000

这个方差很高,如果选择更可靠的初始值卡尔曼滤波收敛会更快。

预测

现在我们基于初始值来预测下一个状态。由于模型是恒定的,预测值就等于当前的估计值:

x^1,0=10oC

状态预测方差的外插为:

p1,0=p0,0+q=10000+0.0001=10000.0001

第1-10次迭代

后续迭代的情况汇总在下面的表里:

n zn 当前状态估计 ( Kn , x^n,n , pn,n ) 预测 ( x^n+1,n , pn+1,n )
1 50.486oC
K1=10000.000110000.0001+0.01=0.999999
x^1,1= 10+0.999999(50.48610)=50.486oC
p1,1=(10.999999)10000.0001=0.01
x^2,1=x^1,1=50.486oC
p2,1=0.01+0.0001=0.0101
2 50.963oC
K2=0.01010.0101+0.01=0.5025
x^2,2= 50.486+0.5025(50.96350.486)=50.726oC
p2,2=(10.5025)0.0101=0.005
x^3,2=x^2,2=50.726oC
p3,2=0.005+0.0001=0.0051
3 51.597oC
K3=0.00510.0051+0.01=0.3388
x^3,3= 50.726+0.3388(51.59750.726)=51.021oC
p3,3=(10.3388)0.0051=0.0034
x^4,3=x^3,3=51.021oC
p4,3=0.0034+0.0001=0.0035
4 52.001oC
K4=0.00350.0035+0.01=0.2586
x^4,4= 51.021+0.2586(52.00151.021)=51.274oC
p4,4=(10.2586)0.0035=0.0026
x^5,4=x^4,4=51.274oC
p5,4=0.0026+0.0001=0.0027
5 52.518oC
K5=0.00270.0027+0.01=0.2117
x^5,5=51.274+0.2117(52.51851.274)=51.538oC
p5,5=(10.2117)0.0027=0.0021
x^6,5=x^5,5=51.538oC
p6,5=0.0021+0.0001=0.0022
6 53.05oC
K6=0.00220.0022+0.01=0.1815
x^6,6= 51.538+0.1815(53.0551.538)=51.812oC
p6,6=(10.1815)0.0022=0.0018
x^7,6=x^6,6=51.812oC
p7,6=0.0018+0.0001=0.0019
7 53.438oC
K7=0.00190.0019+0.01=0.1607
x^7,7= 51.812+0.1607(53.43851.812)=52.0735oC
p7,7=(10.1607)0.0019=0.0016
x^8,7=x^7,7=52.0735oC
p8,7=0.0016+0.0001=0.0017
8 53.858oC
K8=0.00170.0017+0.01=0.1458
x^8,8=52.0735+0.1458(53.85852.0735)=52.334oC
p8,8=(10.1458)0.0017=0.0015
x^9,8=x^8,8=52.334oC
p9,8=0.0015+0.0001=0.0016
9 54.523oC
K9=0.00160.0016+0.01=0.1348
x^9,9= 52.334+0.1348(54.52352.334)=52.621oC
p9,9=(10.1348)0.0016=0.0014
x^10,9=x^9,9=52.621oC
p10,9=0.0014+0.0001=0.0015
10 55.114oC
K10=0.00150.0015+0.01=0.1265
x^10,10= 2.621+0.1265(55.11452.621)=52.936oC
p10,10=(10.1265)0.0015=0.0013
x^11,10=x^10,10=52.936oC
p11,10=0.0013+0.0001=0.0014

结果分析

下图中比较测量值、估计值以及真值。

True value, measured values and estimates

可见,卡尔曼滤波没能给出一个合理的估计,其估计产生了滞后误差。在示例3里,用 αβ 滤波器估计加速中的飞行器时我们已经遇到过一次滞后误差。在示例4里我们通过将 αβ 滤波器更换为能处理加速度估计的 αβγ 滤波器消除了滞后误差。

刚才的卡尔曼滤波中之所以产生滞后误差,有两个原因:

  • 动态模型和实际不符。
  • 我们设定了过低的过程噪声 (q=0.0001),而真实的温度变化要远大于这个值。
注:滞后误差是个稳态误差,故估计值的斜率应该和真值一样。上图只画出了前10次迭代,而10次迭代并不能达到收敛。下图画出了前100次迭代,可以看到稳态的滞后误差。
True value, measured values and estimates

有两个方法能消除滞后误差:

  • 其一,如果我们能预先知道液体温度是线性变化的,就可以设计一个考虑温度线性增长的动态模型。在示例 4中我们就是这样做的,也是一般情况下推荐采用的方法。但是如果温度变化的规律没办法建模,那这种手段就不能改善卡尔曼滤波的性能。
  • 其二,既然是模型失准了,我们可以把模型差异当作一种噪声,可以提高模型的过程噪声 (q)下一个示例会具体讲这种方法。

另一个问题是估计值明明不对,但不确定性却很低(置信区间很窄)。这说明卡尔曼滤波盲目信任预测值,无法给出更准确的估计值。这些现象说明这是一个没设计好的卡尔曼滤波器。

示例小结

本例中,我们用一个面向恒定模型的一维卡尔曼滤波器估计了一个正在加热中的液体的温度。我们观察到卡尔曼滤波的输出产生了滞后误差,这个滞后误差是由错误的动态模型所带来的。

合理的模型设计能够消除滞后误差。

示例8 - 估计加热中的液体的温度

与前一个示例类似,只有一点改变:由于动态模型不是很匹配,我们把过程噪声方差 (q) 从 0.0001 增大到 0.15.

数值示例

第0次迭代

第1次迭代前,我们必须初始化卡尔曼滤波,并且预测下一个时刻的状态(即第1个状态)

初始化

初始化过程和上例类似。

我们不知道缸中液体的真是温度,猜测是10oC.

x^0,0=10oC

这个猜测是不准确的,于是设估计误差 ( σ ) 为 100。初始化的估计方差就是误差方差 (σ2)

p0,0=1002=10,000

这个方差很高,如果选择更可靠的初始值卡尔曼滤波收敛会更快。

Prediction

现在我们基于初始值来预测下一个状态。由于模型是恒定的,预测值就等于当前的估计值:

x^1,0=10oC

状态预测方差的外插为:

p1,0=p0,0+q=10000+0.15=10000.15

第1-10次迭代

后续迭代的情况汇总在下面的表里:

n zn 当前状态估计 ( Kn , x^n,n , pn,n ) 预测 ( x^n+1,n , pn+1,n )
1 50.486oC
K1=10000.1510000.15+0.01=0.999999
x^1,1= 10+0.999999(50.48610)=50.486oC
p1,1=(10.999999)10000.15=0.01
x^2,1=x^1,1=50.486oC
p2,1=0.01+0.15=0.16
2 50.963oC
K2=0.160.16+0.01=0.9412
x^2,2= 50.486+0.9412(50.96350.486)=50.934oC
p2,2=(10.9412)0.16=0.0094
x^3,2=x^2,2=50.934oC
p3,2=0.0094+0.15=0.1594
3 51.597oC
K3=0.15940.1594+0.01=0.941
x^3,3= 50.934+0.941(51.59750.934)=51.556oC
p3,3=(10.941)0.1594=0.0094
x^4,3=x^3,3=51.556oC
p4,3=0.0094+0.15=0.1594
4 52.001oC
K4=0.15940.1594+0.01=0.941
x^4,4= 51.556+0.941(52.00151.556)=51.975oC
p4,4=(10.941)0.1594=0.0094
x^5,4=x^4,4=51.975oC
p5,4=0.0094+0.15=0.1594
5 52.518oC
K5=0.15940.1594+0.01=0.941
x^5,5=51.975+0.941(52.51851.975)=52.486oC
p5,5=(10.941)0.1594=0.0094
x^6,5=x^5,5=52.486oC
p6,5=0.0094+0.15=0.1594
6 53.05oC
K6=0.15940.1594+0.01=0.941
x^6,6= 52.486+0.941(53.0552.486)=53.017oC
p6,6=(10.941)0.1594=0.0094
x^7,6=x^6,6=53.017oC
p7,6=0.0094+0.15=0.1594
7 53.438oC
K7=0.15940.1594+0.01=0.941
x^7,7= 53.017+0.941(53.43853.017)=53.413oC
p7,7=(10.941)0.1594=0.0094
x^8,7=x^7,7=53.413oC
p8,7=0.0094+0.15=0.1594
8 53.858oC
K8=0.15940.1594+0.01=0.941
x^8,8=53.413+0.941(53.85853.413)=53.832oC
p8,8=(10.941)0.1594=0.0094
x^9,8=x^8,8=53.832oC
p9,8=0.0094+0.15=0.1594
9 54.523oC
K9=0.15940.1594+0.01=0.941
x^9,9= 53.832+0.941(54.52353.832)=54.428oC
p9,9=(10.941)0.1594=0.0094
x^10,9=x^9,9=54.428oC
p10,9=0.0094+0.15=0.1594
10 55.114oC
K10=0.15940.1594+0.01=0.941
x^10,10= 54.428+0.941(55.11454.428)=55.074oC
p10,10=(10.941)0.1594=0.0094
x^11,10=x^10,10=55.074oC
p11,10=0.0094+0.15=0.1594

结果分析

下图中比较测量值、估计值以及真值。

True value, measured values and estimates

可以看到,估计值跟上了测量值,滞后误差被消除了。

把过程不确定性加大能够消除稳态误差,但是由于模型仍然是不准的,过程不确定性变大会让卡尔曼滤波更相信测量,因此最终的估计值会几乎等于充满噪声的测量值,整个滤波的目的就达不到了。

我们看看卡尔曼增益的情况。

The Kalman Gain

由于过程不确定性很大,测量值的权重远比估计值的权重高,即卡尔曼增益很高,最终收敛到了0.94.

好消息是本例中测量值都还不错,真值(绿线)基本都落在估计值的95%置信区间内了。

示例小结

一个完美的卡尔曼滤波器要求系统的动态模型非常精确,尽量降低过程噪声。但有时候精确的模型就是无法获得的,比如飞机的飞行员随时都可能突然进行机动,从而大幅改变飞机的轨迹,这是无法预料和建模的,过程噪声就是很大。

上一章 下一章